pythonERA5画水汽通量散度图地图：风速风向矢量图、叠加等高线、色彩分级、添加shp文件、添加位置点及备注

刘志忠

动机有个同事吧，写论文，让我帮忙出个图，就写了个代码，然后我的博客好久没更新了，就顺便贴上来了！很多人感兴趣风速的箭头怎样画，可能这种图使用NCL非常容易，很多没用过代码的小朋友，就有点犯怵，怕python画起来很困难。但是不然，看完我的代码，就会发现很简单，并且也可以批量，同时还能自定义国界等shp文件，这对于发sci等国际论文很重要，因为有时候内置的国界是有问题的。数据本次博客使用的数据为ERA5hourlydataonpressurelevelsfrom1940topresent数据，数据的下载方式及注册账号，我在前面的博客中都写过，详细可参考以下两篇博客：http://t.csdnimg.cn/657dghttp://t.csdnimg.cn/YDELh以下为我们数据介绍界面和需要下载的变量：数据介绍地址：https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels?tab=overview数据选择界面代码废话不多说，直接上代码。导入包importxarrayasxrimportnumpyasnpimportmatplotlib.pyplotaspltimportmatplotlibimportgeopandasasgpd#设置全局字体为新罗马plt.rcParams['font.family']='serif'plt.rcParams['font.serif']=['TimesNewRoman']#plt.rcParams['font.serif']=['SimSun']#设置全局字体权重为normalplt.rcParams['font.weight']='normal'#设置全局字体大小matplotlib.rcParams['font.size']=19#设置全局字体大小为121234567891011121314画水汽通量散度图#加载shapefilegdf=gpd.read_file(r'./shp/Pronvience.shp')#使用geopandas读取地理数据，这里我们手动创建一个GeoDataFramegdf_point=gpd.GeoDataFrame({'City':['MingfengStation','KalasaiStation'],'Latitude':[37.5,37],'Longitude':[80,81]},geometry=gpd.points_from_xy([80,81],[37.5,37]))#载入数据data_path=r'./20170731_case.nc'#替换为您的文件路径ds=xr.open_dataset(data_path)time='2017-07-30T22:00:00'#level_hPa=700#forlevel_hPain[200,500,700,850]:forlevel_hPain[600]:#选择特定时间和气压层ds_selected=ds.sel(time=time,level=level_hPa)#示例：2022年1月1日0时，850hPa#获取数据变量u=ds_selected['u']#东西向风速v=ds_selected['v']#南北向风速q=ds_selected['q']#比湿#获取经度和纬度，假设这些是坐标维度longitude=u.longitudelatitude=u.latitude#计算水汽通量qu=q*u#东西向水汽通量qv=q*v#南北向水汽通量#计算水汽通量散度单位为div_q=(qu.differentiate('longitude')+qv.differentiate('latitude'))*10#打印结果#print(div_q)#创建图形和轴对象fig,ax=plt.subplots(figsize=(6,6),dpi=500)#图形尺寸为10x6英寸#可视化散度结果contour=div_q.plot(add_colorbar=False,cmap="RdBu_r",vmin=-1,vmax=1)#使用黑色线条绘制20个等级的等高线##在ax上绘制等高线图div_q.plot.contour(levels=25,colors='black',linewidths=0.6)#添加颜色条fig.colorbar(contour,ax=ax,label='WaterVaporFluxDivergence(g/cm²/s)')#使用quiver函数需要确保数据的间隔，这里我们每隔5个点取样Q=ax.quiver(longitude[::5],latitude[::5],u[::5,::5],v[::5,::5],scale=300,color="red")#绘制shapefilegdf.plot(ax=ax,color='none',edgecolor='green',linewidths=0.7)#无填充，黑色边界#gdf_point.plot(ax=ax,color='red')#标记纽约的位置#绘制点ax.scatter(gdf_point['Longitude'],gdf_point['Latitude'],color='red',s=100)#标注城市名称forx,y,cityinzip(gdf_point['Longitude'],gdf_point['Latitude'],gdf_point['City']):ax.text(x,y,''+city,verticalalignment='center',fontsize=15)#设置经纬度范围ax.set_xlim(75,90)ax.set_ylim(30,45)ax.set_xlabel('Longitude')ax.set_ylabel('Latitude')ax.set_title('')#清除#添加在图片正下方#fig.suptitle('{}hPa{}'.format(level_hPa,time.replace("T","")),y=-0.01,va='bottom')#调整布局以避免重叠和裁剪fig.tight_layout()plt.savefig("./{}hPa{}.jpg".format(level_hPa,time.replace(":","")),dpi=500)plt.show()1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889水汽通量图#加载shapefilegdf=gpd.read_file(r'./shp/Pronvience.shp')#载入数据data_path=r'./20170731_case.nc'#替换为您的文件路径ds=xr.open_dataset(data_path)time='2017-07-30T22:00:00'forlevel_hPain[200,500,600,700,850]:#选择特定时间和气压层ds_selected=ds.sel(time=time,level=level_hPa)#示例：2022年1月1日0时，850hPa#获取数据变量u=ds_selected['u']#东西向风速v=ds_selected['v']#南北向风速q=ds_selected['q']#比湿#获取经度和纬度，假设这些是坐标维度longitude=u.longitudelatitude=u.latitude#计算水汽通量qu=q*u*100#东西向水汽通量qv=q*v*100#南北向水汽通量wvf=np.sqrt(qu**2+qv**2)#计算水汽通量散度单位为#div_q=(qu.differentiate('longitude')+qv.differentiate('latitude'))*10#打印结果#print(div_q)#创建图形和轴对象fig,ax=plt.subplots(figsize=(6,6),dpi=400)#图形尺寸为10x6英寸#可视化散度结果contour=wvf.plot(add_colorbar=False,cmap="RdBu_r",vmin=0,vmax=10)#使用黑色线条绘制20个等级的等高线##在ax上绘制等高线图wvf.plot.contour(levels=25,colors='black',linewidths=0.6)#添加颜色条fig.colorbar(contour,ax=ax,label='WaterVaporFlux(g/cm/s)')#使用quiver函数需要确保数据的间隔，这里我们每隔5个点取样Q=ax.quiver(longitude[::5],latitude[::5],u[::5,::5],v[::5,::5],scale=300,color="red")#绘制shapefilegdf.plot(ax=ax,color='none',edgecolor='green',linewidths=0.7)#无填充，黑色边界#设置经纬度范围ax.set_xlim(75,90)ax.set_ylim(30,45)ax.set_xlabel('Longitude')ax.set_ylabel('Latitude')ax.set_title('')#清除#添加在图片正下方#fig.suptitle('{}hPa{}'.format(level_hPa,time.replace("T","")),y=-0.01,va='bottom')#调整布局以避免重叠和裁剪fig.tight_layout()plt.savefig("./WVF_{}hPa{}.jpg".format(level_hPa,time.replace(":","")),dpi=500)plt.show()12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667结果图