python计算单层和整层水汽物理量

数学学得特别好

单层和整层的水汽通量及水汽通量散度在ncl下的计算在论坛里已经有很多讨论和经验了，这里总结分享一下python下的计算方法。计算使用的数据是era5，如果使用其他数据计算也都大同小异。主要用到的python库有metpy，geocat。

在ncl中垂直积分有dpres_plevel，vibeta两种计算方法，geocat作为ncar开发的库提供了dpres_plevel的python实现，因此这里使用的是dpres_plevel方法。计算过程和单位的处理参考这个链接：https://cloud.tencent.com/developer/article/1764064 。大家觉得不对的地方欢迎提出交流哈。

2023.09.21更新：geocat-comp库更新后，dpres_plev更新为delta_pressure，不再需要层顶；除此之外修改了图片的量级。

1. import xarray as xr
   
2. import numpy as np
   
3. import cfgrib
   
4. import matplotlib.pyplot as plt
   
5. import matplotlib.ticker as mticker
   
6. import cartopy.crs as ccrs
   
7. import cartopy.feature as cfeature
   
8. from cartopy.io.shapereader import Reader  # 读取 .shp
   
9. import metpy.calc as mpcalc
   
10. import geocat.comp as gc
    
11. import cmaps
    
12. from metpy.units import units
    
13. # 读取省界
    
14. shpName = 'Province_9.shp'
    
15. reader = Reader(shpName)
    
16. provinces = cfeature.ShapelyFeature(reader.geometries(), ccrs.PlateCarree(), edgecolor='k', facecolor='none')
    
17. # 需要的层次
    
18. levels = [1000.,  975.,  950.,  925.,  900.,  850.,  800.,  750.,  700.,  650.,
    
19.           600.,  550.,  500.,  450.,  400.,  350.,  300.,  250.,  200.]
    
20. # 常数
    
21. g = 9.8
    
22. # 读取文件
    
23. files = 'ERA5/' # 这里改成自己的数据路径
    
24. file = 'pl_2014070506.grib'
    
25. file_sl = "sl_2014070506.grib"
    
26. ds = cfgrib.open_datasets(files+file)
    
27. ds_sl = cfgrib.open_datasets(files+file_sl)
    
28. # 读取变量
    
29. u_levels = ds[0].u.sel(isobaricInhPa=levels)
    
30. v_levels = ds[0].v.sel(isobaricInhPa=levels)
    
31. q_levels = ds[0].q.sel(isobaricInhPa=levels) # kg/kg
    
32. psfc = ds_sl[4].sp # 为了在垂直积分时考虑地形，计算时要用到地表气压，单位Pa
    
33. # 接下来是数据的处理：
    
34. plev = u_levels.coords["isobaricInhPa"] * 100 # 单位变为Pa
    
35. plev.attrs['units'] = 'Pa'
    
36. q_levels = q_levels * 1000 # kg/kg -> g/kg
    
37. q_levels.attrs['units'] = 'g/kg'
    
38. # 计算各层次厚度，效果类似ncl中的dpres_plevel
    
39. dp = gc.meteorology.delta_pressure(pressure_lev=plev.data, surface_pressure=psfc.data)
    
40. # Layer Mass Weighting
    
41. dpg = dp / g
    
42. dpg = np.transpose(dpg, (2,0,1)) # 调整dpg维度的顺序，和下面的uq vq保持一致
    
43. # 水汽通量
    
44. uq = u_levels * q_levels
    
45. vq = v_levels * q_levels
    
46. uq.attrs["units"] = "("+u_levels.units+")("+q_levels.units+")"
    
47. vq.attrs["units"] = "("+v_levels.units+")("+q_levels.units+")"
    
48. uq_dpg = uq * dpg.data
    
49. vq_dpg = vq * dpg.data
    
50. # 整层水汽通量
    
51. iuq = uq_dpg.sum(axis=0)
    
52. ivq = vq_dpg.sum(axis=0)
    
53. iuq.attrs["units"] = "[m/s][g/kg]"
    
54. ivq.attrs["units"] = "[m/s][g/kg]"
    
55. # 水汽通量散度
    
56. duvq = mpcalc.divergence(uq, vq)
    
57. duvq.attrs["units"] = "g/(kg-s)"
    
58. duvq_dpg = duvq * dpg.data
    
59. # 整层水汽通量散度
    
60. iduvq = duvq_dpg.sum(axis=0)
    
61. iduvq = iduvq * units('kg/m^2')
    
62. iduvq.attrs["units"] = "g/(m2-s)"
    
63. # 画图
    
64. leftlon, rightlon, lowerlat, upperlat = 80, 100, 20, 40
    
65. img_extent = [leftlon, rightlon, lowerlat, upperlat]
    
66. levs = np.arange(-3, 3+0.5, 0.5)
    
67. 
    
68. fig = plt.figure(figsize=(12,9))
    
69. ax = plt.subplot(111, projection=ccrs.PlateCarree())
    
70. ax.set_extent(img_extent)
    
71. ax.coastlines('50m', linewidth=0.8)
    
72. 
    
73. # 整层水汽通量散度填色
    
74. mfc_contourf = ax.contourf(lons, lats,
    
75.                            iduvq,
    
76.                            cmap=cmaps.MPL_seismic,
    
77.                            levels=levs,
    
78.                            extend='both', transform=ccrs.PlateCarree()) # testcmap precip3_16lev levels=levs,
    
79. 
    
80. # 整层水汽通量流线
    
81. mf_stream = ax.streamplot(lons, lats,
    
82.                           iuq, ivq,
    
83.                           color='black', transform=ccrs.PlateCarree())
    
84. 
    
85. 
    
86. 
    
87. gl = ax.gridlines(draw_labels=True, dms=False,
    
88.                   linestyle='--', color='gray',
    
89.                   x_inline=False, y_inline=False)
    
90. gl.top_labels = False
    
91. gl.right_labels = False
    
92. ax.add_feature(provinces, linewidth=1)
    
93. 
    
94. # 让色标和图等宽
    
95. position = fig.add_axes([ax.get_position().x0,
    
96.                          ax.get_position().y0-0.08,
    
97.                          ax.get_position().width,
    
98.                          0.02])
    
99. cb = fig.colorbar(mfc_contourf, orientation='horizontal', cax=position)
    
100. cb.set_label('%s'%iduvq.units, fontsize=12)

复制代码

最后出的图：

vimfc

notebook文件：

游客，如果您要查看本帖隐藏内容请回复

AuuuTi

顶顶顶！顶顶顶！

zane

谢谢分享。请问ERA5的数据怎么下载的吗？ 谢谢！

随缘

谢谢分享。....

阳光不怕晒

hhhhhh

yuchen

学习了，感谢楼主！

Tao_Xu

道狼52203

学习了，感谢楼主！

cl22

太棒了哈哈哈

kingjames

可以的，很棒的分享！！！