整层水汽通量和整层水汽通量散度计算及python绘图
一、公式推导
1、整层水汽通量：
（1）单层水汽通量：
在P坐标下，
单层水汽通量 = q·v/g
q的单位为kg/kg，v的单位为m/s。对于重力加速度g的单位要进行换算：
也就是说，重力加速度g的单位是10**-2·hPa·m**2/kg。
最终，单层水汽通量的单位为kg/m•hPa•s。
（2）整层水汽通量：
对单层水汽通量进行积分，采用np.trapz。
最终，整层水汽通量的单位为kg/m·s。
2、整层水汽通量散度
（1）单层水汽通量散度：
采用的是mpcalc.divergence。
即：metpy.calc.divergence(u, v, *, dx=None, dy=None, x_dim=- 1, y_dim=- 2)计算矢量的水平散度。
单层水汽通量散度单位为kg/m**2•hPa•s
（2）整层水汽通量散度：
对单层水汽通量散度进行积分，依然使用np.trapz。
为了显示好看，可将最终值提取10**-5或者10**-6次方。
因此整层水汽通量散度的最终单位为：10-5kg/(m**2·s)
二、程序
import matplotlib.pyplot as pltimport matplotlibimport xarray as xrimport cartopy.crs as ccrsimport cartopy.feature as cfimport cartopy.io.shapereader as shpreaderimport cartopy.mpl.ticker as ctickerfrom cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatterfrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERimport matplotlib.ticker as mtickerfrom datetime import datetime#科学计算的包from metpy.units import units#里面是单位import metpy.constants as constants#里面是常数import metpy.calc as mpcalc#里面有各种计算函数plt.rcParams['font.sans-serif'] = ['SimHei']# 用黑体字体显示中文plt.rcParams['axes.unicode_minus']=False# 正常显示负号位置matplotlib.get_cachedir()# Read Datafilename = r'D:\data\physic\201808_physic.nc'f=xr.open_dataset(filename)time = f.time[18:21] # 根据不同的个例选取时间lev = f.level[23:] # 读取气压层，单位为mb，即hPa，一维的14.lat = f.latitude # 读取纬度，一维的21lon = f.longitude# 读取经度，一维的41for i in range(18,21):u = f.u[i,23:,:,:] # U风分量，单位为m/sv = f.v[i,23:,:,:] # V风分量，单位为m/sq = f.q[i,23:,:,:] # 读取比湿，单位为kg/kg# # # 计算单层水汽通量和水汽通量散度qv_u = u*q/(constants.g*10**-2)# g的单位为m/s2，换算为N/kg，再换算为10-2hPa·m2/kg,最终单层水汽通量的单位是kg/m•hPa•sqv_v = v*q/(constants.g*10**-2)# 计算q*v/g,单位是kg/m•hPa•sdx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)# 将经纬度转换为格点距离div_qv = np.zeros((lev.shape[0],lat.shape[0],lon.shape[0]))for j in range(lev.shape[0]):div_qv[j] = mpcalc.divergence(u = qv_u[j],v = qv_v[j],dx = dx ,dy = dy) # 单位是kg/m2•hPa•s# # # 计算整层水汽通量散度total_div_qv = np.trapz(div_qv, lev, axis=0)*10**5#单位为10-5kg/(m**2*s)# # # 计算整层水汽通量total_q_u = np.trapz(qv_u,lev,axis=0) #将单位kg/(m*s)total_q_v = np.trapz(qv_v,lev,axis=0)a = np.sqrt(total_q_u * total_q_u + total_q_v * total_q_v)###绘图levs = np.arange(-1, 1+0.1, 0.1)fig = plt.figure(figsize=(12,9))ax = fig.add_axes([1,0,1,1],projection=ccrs.PlateCarree())ax.set_xticks(np.arange(114, 124, 1), crs=ccrs.PlateCarree())#x刻度值ax.set_yticks(np.arange(34, 39, 0.5), crs=ccrs.PlateCarree())#y刻度值ax.tick_params(axis='both', which='major', labelsize=15) #刻度修饰命令ax.set_extent([114,124,34,39],crs = ccrs.PlateCarree())#绘图范围限制，投影方式为ccrs.PlateCarree()#ax.coastlines('50m', linewidth=0.8)# 绘制水汽通量散度的阴影图,cmap颜色映射表。mfc_contourf = ax.contourf(lon, lat,total_div_qv, cmap='seismic',levels=levs, extend='both', transform=ccrs.PlateCarree())# 绘制水汽通量的箭头图h1 = ax.quiver(lon[::2],lat[::2],total_q_u[::2,::2],total_q_v[::2,::2], #X,Y,U,V 确定位置和对应的风速width = 0.002, #箭杆箭身宽度scale = 700, # 箭杆长度,参数scale越小箭头越长pivot = 'tail'#箭头的其实位置，这里表示从点起，还有点在中心的‘mid’)# 说明箭轴长度与风速的对应关系ax.quiverkey(h1,#传入quiver句柄X= 0.1, Y = -0.07, #确定 label 所在位置，都限制在[0,1]之间U = 20, #参考箭头长度 表示20。angle = 0,#参考箭头摆放角度。默认为0，即水平摆放 label='20m/s', #箭头的补充：label的内容+ labelpos='E',#label在参考箭头的哪个方向; S表示南边 color = 'k',labelcolor = 'k', #箭头颜色 + label的颜色 )# 绘制水汽通量的等值线ct=ax.contour(lon,lat,a,8,colors='k',linewidths=1,levels=range(0,28,2))# 标记ct每根线条的数值ax.clabel(ct,inline=True,fontsize=10,fmt='%.0f') # 绘制山东省界province = shpreader.Reader(r'D:\shp\Shandong-city-2020\Shandong-city-2020.shp')ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5,edgecolor='k',facecolor='none')ax.add_feature# 图例，颜色与数值的对应关系。orientation：colorbar摆放的横竖位置。cax：colorbar放在指定位置，最高优先级。position = fig.add_axes([ax.get_position().x0, ax.get_position().y0-0.08, ax.get_position().width, 0.02])cb = fig.colorbar(mfc_contourf, orientation='horizontal', cax=position)#cb.set_label('g/(m**2*s)', fontsize=12)#fig.savefig(r'D:\py_pic\整层水汽通量和整层水汽通量散度图.jpg', bbox_inches = 'tight')