整層水汽通量和整層水汽通量散度計(jì)算及python繪圖
一、公式推導(dǎo)
1、整層水汽通量：
（1）單層水汽通量：
在P坐標(biāo)下，
單層水汽通量 = q·v/g
q的單位為kg/kg，v的單位為m/s。對(duì)于重力加速度g的單位要進(jìn)行換算：

也就是說，重力加速度g的單位是10**-2·hPa·m**2/kg。
最終，單層水汽通量的單位為kg/m?hPa?s。

（2）整層水汽通量：
對(duì)單層水汽通量進(jìn)行積分，采用np.trapz。
最終，整層水汽通量的單位為kg/m·s。

2、整層水汽通量散度
（1）單層水汽通量散度：

采用的是mpcalc.divergence。
即：metpy.calc.divergence(u, v, *, dx=None, dy=None, x_dim=- 1, y_dim=- 2)計(jì)算矢量的水平散度。
單層水汽通量散度單位為kg/m**2?hPa?s

（2）整層水汽通量散度：
對(duì)單層水汽通量散度進(jìn)行積分，依然使用np.trapz。
為了顯示好看，可將最終值提取10**-5或者10**-6次方。
因此整層水汽通量散度的最終單位為：10-5kg/(m**2·s)

二、程序文章來源地址http://www.zghlxwxcb.cn/news/detail-432039.html

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
#科學(xué)計(jì)算的包
from metpy.units import units      #里面是單位            
import metpy.constants as constants  #里面是常數(shù)
import metpy.calc as mpcalc          #里面有各種計(jì)算函數(shù)
plt.rcParams['font.sans-serif'] = ['SimHei']                                    # 用黑體字體顯示中文
plt.rcParams['axes.unicode_minus']=False                                        # 正常顯示負(fù)號(hào)位置
matplotlib.get_cachedir()

# Read Data
filename = r'D:\data\physic\201808_physic.nc'
f=xr.open_dataset(filename)

time = f.time[18:21]       # 根據(jù)不同的個(gè)例選取時(shí)間
lev = f.level[23:]         # 讀取氣壓層，單位為mb，即hPa，一維的14.
lat = f.latitude           # 讀取緯度，一維的21
lon = f.longitude          # 讀取經(jīng)度，一維的41

for i in range(18,21):
    u = f.u[i,23:,:,:]         # U風(fēng)分量，單位為m/s
    v = f.v[i,23:,:,:]         # V風(fēng)分量，單位為m/s
    q = f.q[i,23:,:,:]         # 讀取比濕，單位為kg/kg
# # # 計(jì)算單層水汽通量和水汽通量散度
    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)                            # 計(jì)算q*v/g,單位是kg/m?hPa?s
    dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)              # 將經(jīng)緯度轉(zhuǎn)換為格點(diǎn)距離

    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
        
# # # 計(jì)算整層水汽通量散度

    total_div_qv = np.trapz(div_qv, lev, axis=0)*10**5    #單位為10-5kg/(m**2*s)

# # # 計(jì)算整層水汽通量
    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 確定位置和對(duì)應(yīng)的風(fēng)速
                width = 0.002, #箭桿箭身寬度  
                scale = 700, # 箭桿長度,參數(shù)scale越小箭頭越長
                pivot = 'tail'#箭頭的其實(shí)位置，這里表示從點(diǎn)起，還有點(diǎn)在中心的‘mid’
              )
# 說明箭軸長度與風(fēng)速的對(duì)應(yīng)關(guān)系
ax.quiverkey(h1,                      #傳入quiver句柄
              X= 0.1, Y = -0.07,       #確定 label 所在位置，都限制在[0,1]之間
              U = 20,                 #參考箭頭長度 表示20。
              angle = 0,              #參考箭頭擺放角度。默認(rèn)為0，即水平擺放
             label='20m/s',           #箭頭的補(bǔ)充：label的內(nèi)容  +
             labelpos='E',            #label在參考箭頭的哪個(gè)方向; S表示南邊
             color = 'k',labelcolor = 'k', #箭頭顏色 + label的顏色
             )


# 繪制水汽通量的等值線
ct=ax.contour(lon,lat,a,8,colors='k',linewidths=1,levels=range(0,28,2))
# 標(biāo)記ct每根線條的數(shù)值
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

# 圖例，顏色與數(shù)值的對(duì)應(yīng)關(guān)系。orientation：colorbar擺放的橫豎位置。cax：colorbar放在指定位置，最高優(yōu)先級(jí)。
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')

到了這里，關(guān)于整層水汽通量和整層水汽通量散度計(jì)算及python繪圖的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！