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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 31793|回复: 12

[源代码] python绘制GPM CSH潜热垂直剖面

[复制链接]

新浪微博达人勋

发表于 2019-7-13 10:48:45 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 灭火器 于 2019-7-13 10:48 编辑

使用h5py包来读取GPM CSH的潜热数据,然后用matplotlib来绘制填色图。
麻烦的地方在于,如何在潜热数据正负大小不对称的前提下,让潜热为零与红蓝colormap中的白色相对应。这里参考了matplotlib官网给出的解决方法,定义一个MidpointNormalize类,把原本在[vmin, vmax]内线性的normalize调整为在[vmin, midpoint]和[midpoint, vmax]两个区间内分别线性归一化。参考资料如下:
[1] https://www.eorc.jaxa.jp/GPM/doc/program/guide/GPM_read_program_guide_forPython_V5.1.pdf
[2] https://matplotlib.org/users/colormapnorms.html#custom-normalization-two-linear-ranges
仍然需要改进的地方:实现任意剖面、叠加地形……
  1. import numpy as np
  2. import matplotlib as mpl
  3. import matplotlib.pyplot as plt
  4. import h5py

  5. # 读取经纬度和潜热.
  6. filename = "2B.GPM.DPRGMI.2HCSHv3-1.20170503-S070958-E084232.018057.V05A.HDF5"
  7. group = "Swath"
  8. with h5py.File(filename, "r") as infile:
  9.     lat2d = np.array(infile[group+"/Latitude"])
  10.     lon2d = np.array(infile[group+"/Longitude"])
  11.     LH3d = np.array(infile[group+"/latentHeating"])
  12.    
  13.     nscan = LH3d.shape[0]
  14.     nray = LH3d.shape[1]
  15.     nlayer = LH3d.shape[2]

  16. # 选取所需的经纬度范围,还有想要画的ray.
  17. ray = int(nray/2) + 20
  18. lonmin = 105.0
  19. lonmax = 120.0
  20. latmin = 25.0
  21. latmax = 40.0

  22. # 根据官方文档设置对应的垂直高度.
  23. layerSize = 0.25    # Units is km.
  24. height = (np.arange(1, nlayer+1)-0.5)*layerSize

  25. # 利用条件判断索引所需的潜热截面数据.
  26. logi = np.logical_and.reduce((lat2d[:,ray] >= latmin, lat2d[:,ray] <= latmax, \
  27.                               lon2d[:,ray] >= lonmin, lon2d[:,ray] <= lonmax))
  28. lat = lat2d[logi,ray]
  29. LH = LH3d[logi,ray,:]
  30. print(LH.min(), LH.max())

  31. # 定义一个能保证midpoint在colorbar正中间的normalize类.
  32. class MidpointNormalize(mpl.colors.Normalize):
  33.     def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
  34.         self.midpoint = midpoint
  35.         mpl.colors.Normalize.__init__(self, vmin, vmax, clip)

  36.     def __call__(self, value, clip=None):
  37.         x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
  38.         return np.ma.masked_array(np.interp(value, x, y))

  39. # 做出潜热截面的填色图.
  40. fig = plt.figure(dpi=300, figsize=(10,4))
  41. ax = plt.subplot(111)
  42. bounds = [-5.0, -2.0, -1.0, -0.1, 0.1, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0]
  43. norm = mpl.colors.BoundaryNorm(boundaries=bounds, ncolors=256)
  44. im = ax.contourf(lat, height, LH.T, levels=bounds, cmap="seismic", extend="both", \
  45.                  norm=MidpointNormalize(midpoint=0.0, vmin=min(bounds), vmax=max(bounds)))
  46. ax.axhline(y=4.0, color="k", linestyle="--")

  47. ax.set_xlabel("Latitude", fontsize=12)
  48. ax.set_ylabel("Altitude (km)", fontsize=12)
  49. ax.set_yticks(np.arange(25, step=5))
  50. ax.set_title("CSH Latent Heat (K/hr)", fontsize=15)

  51. cbar = fig.colorbar(im, ax=ax, ticks=bounds, orientation="horizontal", \
  52.                     shrink=0.7, pad=0.2)

  53. # plt.show()
  54. fig.savefig("CSH_plot_py")
复制代码

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

新浪微博达人勋

发表于 2019-10-20 12:17:44 | 显示全部楼层
你好,请问这个剖面是沿哪个经度剖的吗,没看懂你29-34行代码,求楼主解答
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-10-20 13:13:43 | 显示全部楼层
不对称无所谓啊感觉,不过归一化是一种方法,但是觉得楼主你这个colorbar还是很奇怪,颜色用的
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-10-20 14:44:46 | 显示全部楼层
我的明天 发表于 2019-10-20 12:17
你好,请问这个剖面是沿哪个经度剖的吗,没看懂你29-34行代码,求楼主解答

我是选取了一个包含卫星扫描的经纬度方框,画nray正中间的那个剖面。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-10-24 15:11:42 | 显示全部楼层
灭火器 发表于 2019-10-20 14:44
我是选取了一个包含卫星扫描的经纬度方框,画nray正中间的那个剖面。

谢谢楼主的解答,我还有些问题
1.nray正中间的剖面,不应该是ray = int(nray/2)吗,为什么是ray = int(nray/2) + 20,不知楼主的nray是不是等于49
2.我想把楼主的程序改了下套在我画GPM反射率的剖面上,但是画出来没数据,不知是怎么回事,顺便问下,楼主现在实现任意剖面了没
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-10-24 16:39:28 | 显示全部楼层
我的明天 发表于 2019-10-24 15:11
谢谢楼主的解答,我还有些问题
1.nray正中间的剖面,不应该是ray = int(nray/2)吗,为什么是ray = int(n ...

呃,那我估计画的是别的剖面,反正ray可以任意取1~49之间的数。反射率我画过,没问题啊。任意剖面还没实现过(不过你可以考虑用最邻近插值来取点)
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-10-26 15:37:04 | 显示全部楼层
灭火器 发表于 2019-10-24 16:39
呃,那我估计画的是别的剖面,反正ray可以任意取1~49之间的数。反射率我画过,没问题啊。任意剖面还没实 ...

我只改了你的路径,经纬度选的这个
ray = int(nray/2)+20
lonmin = 112
lonmax = 113.2
latmin = 22.0
latmax = 23.0
把你的cmap去掉了normalize这类的改成了简单的cmap=plt.get_cmap('RdBu'),运行的时候,im = ax.contourf(lat, height, LH.T, cmap=plt.get_cmap('RdBu'), extend="both")这句error: Input z must be at least a 2x2 array.麻烦楼主可以帮我看下是怎么回事吗?
下边是我改的
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import h5py

# 读取经纬度和反射率.
filename = "/mnt/e/2A.GPM.DPR.V8-20180723.20160709-S214726-E231959.013431.V06A.HDF5"

with h5py.File(filename, "r") as f:
    name = '/NS/SLV/zFactorCorrected'
    LH3d = f[name][:,:,:]

    lat2d = f['/NS/Latitude'][:]
    lon2d = f['/NS/Longitude'][:]
   
    nscan = LH3d.shape[0]
    nray = LH3d.shape[1]
    nlayer = LH3d.shape[2]

ray = int(nray/2)
lonmin = 112
lonmax = 113.2
latmin = 22.0
latmax = 23.0

# 根据官方文档设置对应的垂直高度.
layerSize = 0.125    # Units is km.
height = (np.arange(0, nlayer))*layerSize

# 利用条件判断索引所需的潜热截面数据.
logi = np.logical_and.reduce((lat2d[:,ray] >= latmin, lat2d[:,ray] <= latmax, \
                              lon2d[:,ray] >= lonmin, lon2d[:,ray] <= lonmax))
lat = lat2d[logi,ray]
LH = LH3d[logi,ray,:]

fig = plt.figure()
ax = plt.subplot(111)

im = ax.contourf(lat, height, LH.T, cmap=plt.get_cmap('RdBu'), extend="both")
ax.axhline(y=4.0, color="k", linestyle="--")

ax.set_xlabel("Latitude", fontsize=12)
ax.set_ylabel("Altitude (km)", fontsize=12)
ax.set_yticks(np.arange(25, step=5))
ax.set_title("CSH Latent Heat (K/hr)", fontsize=15)

cbar = fig.colorbar(im, ax=ax,  orientation="horizontal", \
                    shrink=0.7, pad=0.2)

# plt.show()
fig.savefig("CSH_plot_py")
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2019-10-27 09:48:33 | 显示全部楼层
我的明天 发表于 2019-10-26 15:37
我只改了你的路径,经纬度选的这个
ray = int(nray/2)+20
lonmin = 112

不知道。但至少有一点,2A DPR产品中nbin从1到176,高度是从高到低。而SLH产品中nlayer从1到80,高度是从低到高。你这个会画反。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2019-10-28 10:00:29 | 显示全部楼层
灭火器 发表于 2019-10-27 09:48
不知道。但至少有一点,2A DPR产品中nbin从1到176,高度是从高到低。而SLH产品中nlayer从1到80,高度是从 ...

嗯嗯,谢谢楼主的耐心解答~
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-4-9 15:24:50 | 显示全部楼层
想问问楼主现在知道怎么在垂直剖面图上添加地形了吗~之前画山脉附近的水汽垂直剖面图不知道咋加上地形
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

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

本版积分规则

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

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

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