python绘制GPM CSH潜热垂直剖面

灭火器

使用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. # 读取经纬度和潜热.
   
7. filename = "2B.GPM.DPRGMI.2HCSHv3-1.20170503-S070958-E084232.018057.V05A.HDF5"
   
8. group = "Swath"
   
9. with h5py.File(filename, "r") as infile:
   
10.     lat2d = np.array(infile[group+"/Latitude"])
    
11.     lon2d = np.array(infile[group+"/Longitude"])
    
12.     LH3d = np.array(infile[group+"/latentHeating"])
    
13.    
    
14.     nscan = LH3d.shape[0]
    
15.     nray = LH3d.shape[1]
    
16.     nlayer = LH3d.shape[2]
    
17. 
    
18. # 选取所需的经纬度范围,还有想要画的ray.
    
19. ray = int(nray/2) + 20
    
20. lonmin = 105.0
    
21. lonmax = 120.0
    
22. latmin = 25.0
    
23. latmax = 40.0
    
24. 
    
25. # 根据官方文档设置对应的垂直高度.
    
26. layerSize = 0.25    # Units is km.
    
27. height = (np.arange(1, nlayer+1)-0.5)*layerSize
    
28. 
    
29. # 利用条件判断索引所需的潜热截面数据.
    
30. logi = np.logical_and.reduce((lat2d[:,ray] >= latmin, lat2d[:,ray] <= latmax, \
    
31.                               lon2d[:,ray] >= lonmin, lon2d[:,ray] <= lonmax))
    
32. lat = lat2d[logi,ray]
    
33. LH = LH3d[logi,ray,:]
    
34. print(LH.min(), LH.max())
    
35. 
    
36. # 定义一个能保证midpoint在colorbar正中间的normalize类.
    
37. class MidpointNormalize(mpl.colors.Normalize):
    
38.     def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False):
    
39.         self.midpoint = midpoint
    
40.         mpl.colors.Normalize.__init__(self, vmin, vmax, clip)
    
41. 
    
42.     def __call__(self, value, clip=None):
    
43.         x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1]
    
44.         return np.ma.masked_array(np.interp(value, x, y))
    
45. 
    
46. # 做出潜热截面的填色图.
    
47. fig = plt.figure(dpi=300, figsize=(10,4))
    
48. ax = plt.subplot(111)
    
49. bounds = [-5.0, -2.0, -1.0, -0.1, 0.1, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0]
    
50. norm = mpl.colors.BoundaryNorm(boundaries=bounds, ncolors=256)
    
51. im = ax.contourf(lat, height, LH.T, levels=bounds, cmap="seismic", extend="both", \
    
52.                  norm=MidpointNormalize(midpoint=0.0, vmin=min(bounds), vmax=max(bounds)))
    
53. ax.axhline(y=4.0, color="k", linestyle="--")
    
54. 
    
55. ax.set_xlabel("Latitude", fontsize=12)
    
56. ax.set_ylabel("Altitude (km)", fontsize=12)
    
57. ax.set_yticks(np.arange(25, step=5))
    
58. ax.set_title("CSH Latent Heat (K/hr)", fontsize=15)
    
59. 
    
60. cbar = fig.colorbar(im, ax=ax, ticks=bounds, orientation="horizontal", \
    
61.                     shrink=0.7, pad=0.2)
    
62. 
    
63. # plt.show()
    
64. fig.savefig("CSH_plot_py")

复制代码

我的明天

你好，请问这个剖面是沿哪个经度剖的吗，没看懂你29-34行代码，求楼主解答

xiaofeixiazyh

不对称无所谓啊感觉，不过归一化是一种方法，但是觉得楼主你这个colorbar还是很奇怪，颜色用的

灭火器

我的明天 发表于 2019-10-20 12:17
你好，请问这个剖面是沿哪个经度剖的吗，没看懂你29-34行代码，求楼主解答

我是选取了一个包含卫星扫描的经纬度方框，画nray正中间的那个剖面。

我的明天

灭火器 发表于 2019-10-20 14:44
我是选取了一个包含卫星扫描的经纬度方框，画nray正中间的那个剖面。

谢谢楼主的解答，我还有些问题
1.nray正中间的剖面，不应该是ray = int(nray/2)吗，为什么是ray = int(nray/2) + 20，不知楼主的nray是不是等于49
2.我想把楼主的程序改了下套在我画GPM反射率的剖面上，但是画出来没数据，不知是怎么回事，顺便问下，楼主现在实现任意剖面了没

灭火器

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

呃，那我估计画的是别的剖面，反正ray可以任意取1~49之间的数。反射率我画过，没问题啊。任意剖面还没实现过（不过你可以考虑用最邻近插值来取点）

我的明天

灭火器 发表于 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")

灭火器

我的明天 发表于 2019-10-26 15:37
我只改了你的路径，经纬度选的这个
ray = int(nray/2)+20
lonmin = 112

不知道。但至少有一点，2A DPR产品中nbin从1到176，高度是从高到低。而SLH产品中nlayer从1到80，高度是从低到高。你这个会画反。

我的明天

灭火器 发表于 2019-10-27 09:48
不知道。但至少有一点，2A DPR产品中nbin从1到176，高度是从高到低。而SLH产品中nlayer从1到80，高度是从 ...

嗯嗯，谢谢楼主的耐心解答~

bibibol

想问问楼主现在知道怎么在垂直剖面图上添加地形了吗~之前画山脉附近的水汽垂直剖面图不知道咋加上地形