请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 115506|回复: 407

[经验总结] python计算单层和整层水汽物理量

  [复制链接]

新浪微博达人勋

发表于 2021-10-18 20:53:25 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 数学学得特别好 于 2023-10-6 19:38 编辑

单层和整层的水汽通量及水汽通量散度在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. fig = plt.figure(figsize=(12,9))
  68. ax = plt.subplot(111, projection=ccrs.PlateCarree())
  69. ax.set_extent(img_extent)
  70. ax.coastlines('50m', linewidth=0.8)

  71. # 整层水汽通量散度填色
  72. mfc_contourf = ax.contourf(lons, lats,
  73.                            iduvq,
  74.                            cmap=cmaps.MPL_seismic,
  75.                            levels=levs,
  76.                            extend='both', transform=ccrs.PlateCarree()) # testcmap precip3_16lev levels=levs,

  77. # 整层水汽通量流线
  78. mf_stream = ax.streamplot(lons, lats,
  79.                           iuq, ivq,
  80.                           color='black', transform=ccrs.PlateCarree())



  81. gl = ax.gridlines(draw_labels=True, dms=False,
  82.                   linestyle='--', color='gray',
  83.                   x_inline=False, y_inline=False)
  84. gl.top_labels = False
  85. gl.right_labels = False
  86. ax.add_feature(provinces, linewidth=1)

  87. # 让色标和图等宽
  88. position = fig.add_axes([ax.get_position().x0,
  89.                          ax.get_position().y0-0.08,
  90.                          ax.get_position().width,
  91.                          0.02])
  92. cb = fig.colorbar(mfc_contourf, orientation='horizontal', cax=position)
  93. cb.set_label('%s'%iduvq.units, fontsize=12)
复制代码

最后出的图:

vimfc

vimfc

notebook文件:
游客,如果您要查看本帖隐藏内容请回复







评分

参与人数 1金钱 +10 收起 理由
木木羊 + 10

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 成长值: 0
发表于 2021-10-18 21:05:24 | 显示全部楼层
顶顶顶!顶顶顶!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-10-19 07:51:50 | 显示全部楼层
谢谢分享。请问ERA5的数据怎么下载的吗? 谢谢!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-10-19 08:18:02 | 显示全部楼层
谢谢分享。....
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-10-19 08:36:27 | 显示全部楼层
hhhhhh
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-10-19 08:56:25 | 显示全部楼层
学习了,感谢楼主!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-10-19 09:09:10 | 显示全部楼层
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2021-10-19 09:12:28 | 显示全部楼层
学习了,感谢楼主!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-10-19 09:15:55 | 显示全部楼层
太棒了哈哈哈
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2021-10-19 09:43:46 | 显示全部楼层
可以的,很棒的分享!!!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表