Python计算单层、整层水汽物理量(水汽通量和水汽通量散度)

                                                                                                   

文章转载自气象家园：

http://bbs.06climate.com/forum.php?mod=viewthread&tid=100595&extra=&page=1

---

作者联系邮箱：

ozziesun@outlook.com

---

单层和整层的水汽物理量在ncl下的计算在论坛里已经有很多讨论和经验了，这里总结分享一下python下的计算方法。计算使用的数据是era5，如果使用其他数据计算也都是大同小异的。主要用到的python库有metpy，geocat。
, Y& F5 l2 k  a/ Z5 k& Q在ncl中垂直积分有dpres_plevel，vibeta两种计算方法，geocat作为ncar开发的库提供了dpres_plevel的python实现，因此这里使用的是dpres_plevel方法。计算过程和单位的处理参考这个链接：https://cloud.tencent.com/developer/article/1764064 。大家觉得不对的地方欢迎提出交流哈。


* U! P4 H0 h# B

[Python] 纯文本查看 复制代码

import xarray as xr
import numpy as np
import cfgrib
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader  # 读取 .shp
import metpy.calc as mpcalc
import geocat.comp
import cmaps
from metpy.units import units
# 读取省界
shpName = 'Province_9.shp'
reader = Reader(shpName)
provinces = cfeature.ShapelyFeature(reader.geometries(), ccrs.PlateCarree(), edgecolor='k', facecolor='none')
# 需要的层次
levels = [1000.,  975.,  950.,  925.,  900.,  850.,  800.,  750.,  700.,  650.,
          600.,  550.,  500.,  450.,  400.,  350.,  300.,  250.,  200.]
# 常数
ptop = 200*100 # 垂直积分层顶，单位Pa
g = 9.8
# 读取文件
files = 'ERA5/' # 这里改成自己的数据路径
file = 'pl_2014070506.grib'
file_sl = "sl_2014070506.grib"
ds = cfgrib.open_datasets(files+file)
ds_sl = cfgrib.open_datasets(files+file_sl)
# 读取变量
u_levels = ds[0].u.sel(isobaricInhPa=levels)
v_levels = ds[0].v.sel(isobaricInhPa=levels)
q_levels = ds[0].q.sel(isobaricInhPa=levels) # kg/kg
psfc = ds_sl[4].sp # 地表气压，单位Pa
# 接下来是数据的处理：
plev = u_levels.coords["isobaricInhPa"] * 100 # 单位变为Pa
plev.attrs['units'] = 'Pa'
q_levels = q_levels * 1000 # kg/kg -> g/kg
q_levels.attrs['units'] = 'g/kg'
# 垂直积分的关键，计算各层次厚度，效果类似ncl中的dpres_plevel
dp = geocat.comp.dpres_plevel(plev, psfc, ptop)
# Layer Mass Weighting
dpg = dp / g
dpg.attrs["units"] = "kg/m2"
# 水汽通量计算：
# 水汽通量
uq = u_levels * q_levels
vq = v_levels * q_levels
uq.attrs["units"] = "("+u_levels.units+")("+q_levels.units+")" 
vq.attrs["units"] = "("+v_levels.units+")("+q_levels.units+")" 
uq_dpg = uq * dpg.data
vq_dpg = vq * dpg.data
# 整层水汽通量
iuq = uq_dpg.sum(axis=0)
ivq = vq_dpg.sum(axis=0)
iuq.attrs["units"] = "[m/s][g/kg]"
ivq.attrs["units"] = "[m/s][g/kg]"
# 水汽通量散度：
# 计算散度要用到dx，dy
lons = u_levels.coords['longitude'][:]
lats = v_levels.coords['latitude'][:]
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)
# 水汽通量散度
duvq = mpcalc.divergence(uq/g, vq/g, dx=dx[np.newaxis,:,:], dy=dy[np.newaxis,:,:]) 
duvq.attrs["units"] = "g/(kg-s)"
duvq_dpg = duvq * dpg.data
# 整层水汽通量散度
iduvq = duvq_dpg.sum(axis=0)
iduvq = iduvq * units('kg/m^2')
iduvq.attrs["units"] = "g/(m2-s)"
# 画图：
leftlon, rightlon, lowerlat, upperlat = 80, 100, 20, 40
img_extent = [leftlon, rightlon, lowerlat, upperlat]
levs = np.arange(-0.3, 0.3+0.05, 0.05)
fig = plt.figure(figsize=(12,9))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
ax.set_extent(img_extent)
ax.coastlines('50m', linewidth=0.8)
# 整层水汽通量散度填色
mfc_contourf = ax.contourf(lons, lats, 
                           iduvq,
                           cmap=cmaps.MPL_seismic, 
                           levels=levs,
                           extend='both', transform=ccrs.PlateCarree())
# 整层水汽通量流线
mf_stream = ax.streamplot(lons, lats,
                          iuq.data, ivq.data,
                          color='black', transform=ccrs.PlateCarree())
gl = ax.gridlines(draw_labels=True, dms=False, 
                  linestyle='--', color='gray',
                  x_inline=False, y_inline=False)
gl.top_labels = False
gl.right_labels = False
ax.add_feature(provinces, linewidth=1)
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('%s'%iduvq.units, fontsize=12)
# 保存图片
plt.savefig('./plo99.png',dpi=800,bbox_inches='tight',pad_inches=0)
plt.show()

3 G- D+ W. }6 r2 x8 q3 }+ j
- s" P* X+ ?3 g! n& g) C$ T

! m  W% w, `3 C, {9 t& }2 L9 L

---

+ @9 E7 j& L" p

如有问题,进入论坛与原作者交流。