整层水汽通量和整层水汽通量散度计算及python绘图

整层水汽通量和整层水汽通量散度计算及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。

整层水汽通量散度
（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 cticker  from 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/s    v = f.v[i,23:,:,:]         # V风分量，单位为m/s    q = 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•s    qv_v = v*q/(constants.g*10**-2)# 计算q*v/g,单位是kg/m•hPa•s    dx, 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')

来源地址：https://blog.csdn.net/wdbhysszjswn/article/details/129044910