利用python绘制EC数据全国降水图

zhaoyu

学生

推荐：
1、 Python在气象与海洋中的实践技术应用

2、 Python在WRF模型自动化运行及前后处理中的实践技术应用

3、 python机器学习与深度学习在气象领域中的应用

4、【夏令营】针对课题组人员AI培养计划:“开启AI科研之路”

5、 Python语言在地球科学领域中的应用

Python的绘图功能非常强大，在大气和海洋常常用来绘制一些有关地理方面的图。本片主要介绍python绘制EC数据（grib格式）的的全国降水分布图。

总体流程：首先利用pygrib读取grib数据中的累计降水量数据，然后利用matplotlib进行绘图，其中利用cartopy进行地图shp文件的处理，并绘制南海小地图。

代码如下：

#导入需要的包，缺什么就pip安装什么,主要依赖matplotlib、pygrib、cartopy三个包。
import matplotlib
matplotlib.use('Agg')
import pygrib as pg
import matplotlib.pyplot as plt
import cartopy.io.shapereader as shpreader
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
from pylab import * 
import matplotlib.colors as col
import matplotlib.cm as cm
from cartopy.mpl.gridliner import \
  LONGITUDE_FORMATTER, LATITUDE_FORMATTER
#读取文件，grib格式，全国的降水预报
ec_file='H:/2021/2021-08-16/D2S08160000081700001'
#读取变量
data=pg.open(ec_file)
rain_data = data.select(name='Total precipitation')[0]
rain= rain_data.values*1000 #转化单位
lat,lon=rain_data.latlons() #获取经纬度
#开始画图
proj = ccrs.PlateCarree()  # 创建投影
fig = plt.figure(figsize=(18, 18),dpi=800) #设置图片大小和分辨率
ax = fig.subplots(1, 1, subplot_kw={'projection':proj})
# 中国经纬度范围
region = [70, 136, 15, 55]
ax.set_extent(region, crs=proj)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),linewidth=0.6, zorder=1)
ax.add_feature(cfeature.RIVERS.with_scale('50m'),zorder=1)
ax.add_feature(cfeature.LAKES.with_scale('50m'),zorder=1)
ax.set_title('EC:Total Precipitation(mm)',loc='left',fontsize=12)
time='Time:'+' 00:00:00(UTC)+24hours'
ax.set_title(time, loc='right',fontsize=12)
gl = ax.gridlines(ylocs=np.arange(region[2],region[3] +10, 10), xlocs=np.arange(region[0], region[1] + 10, 10),
draw_labels=True, linestyle='--', alpha=0.7)
gl.xlabels_top = False
gl.ylabels_right = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
#cpool = ['#FFFFFF', '#98FB98', '#32CD32',  '#00BFFF','#0000FF', '#FF00FF', '#800000']
#cmap3 = col.ListedColormap(cpool, 'indexed')
# 注意：cpool中有15个颜色，但是只使用了10个，其余的会被忽略
# 具体可查阅matplotlib官方文档中对ListedColormap的说明
#cm.register_cmap(name='dbzcmap', cmap=cmap3)
levels = [0,0.1,10,25,50,100,250]
#contours = plt.contourf(lon, lat, rain,levels=levels,extend='both',cmap='Accent',transform=ccrs.PlateCarree())
contours = plt.contourf(lon, lat,rain,levels=levels,extend='both',colors=('#FFFFFF', '#98FB98', '#32CD32',  '#00BFFF','#0000FF', '#FF00FF', '#800000'))
cbar = fig.colorbar(contours,shrink=0.4)#,ticks=levels)
cbar.set_label('Total Precipitation (mm)',fontsize=12)
# 获取国家基础地理数据，加载全国shp文件
#后台回复chinashp，获取全国的shapefile文件。
china = shpreader.Reader('H:/geo_data/国家基础地理数据400万/省.shp').geometries()
pro=shpreader.Reader('H:/geo_data/国家基础地理数据400万/国界线.shp').geometries()
# 绘制中国国界省界九段线等等
ax.add_geometries(china, proj,facecolor='none',edgecolor='black', zorder=1)
# 添加南海
sub_ax = fig.add_axes([0.115, 0.315, 0.1, 0.1], projection=proj)
sub_ax.set_extent([105, 125, 0,25],crs=ccrs.PlateCarree())
sub_ax.add_feature(cfeature.COASTLINE.