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

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 plt
import matplotlib
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cf
import cartopy.io.shapereader as shpreader
import cartopy.mpl.ticker as cticker  
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
from 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 Data
filename = 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           # 读取纬度，一维的21
lon = f.longitude          # 读取经度，一维的41

for 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')