GPM 卫星IMERG降水数据的简单可视化

tosaka

公众号：气python风雨
和鲸在线运行链接：GPM卫星简单可视化 - Heywhale.com
GPM IMERG的简单可视化

GPM数据之前写了下载教程，现在简单试试可视化，毕竟是nc格式数据（下载可选），用起来相对简单

1.1 导入库与查看数据

1. # 导入库
   
2. import numpy as np
   
3. import xarray as xr
   
4. #可视化
   
5. import matplotlib.pyplot as plt
   
6. import cartopy.crs as ccrs
   
7. import cartopy.feature as cfeat
   
8. import cartopy.io.shapereader as shpreader
   
9. from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
   
10. import cmaps
    
11. from shapely import geometry as sgeom
    
12. # 辅助模块
    
13. from glob import glob
    
14. from copy import copy

复制代码

1. # 读取
   
2. ds = xr.open_dataset('/home/mw/input/GPM5633/gpmPython/3B-DAY.MS.MRG.3IMERG.20210519-S000000-E235959.V06.nc4')
   
3. print(ds)

复制代码

                               
登录/注册后可看大图

可见该数据为2021年5月19日的日降水数据，经纬度分辨率为0.1°
ProductionTime指的是产品生成时间，因为这是三级数据，在三个月后生成是正常的。具体可见官网解释产品分级。


1.2 IMERG降水可视化

加入修正cartopy无法正常显示经纬度的函数

1. '''
   
2. 参考Andrew Dawson 提供的解决方法: https://gist.github.com/ajdawson/dd536f786741e987ae4e
   
3. '''
   
4. # 给出一个假定为矩形的LineString，返回对应于矩形给定边的直线。
   
5. def find_side(ls, side):
   
6.     """
   
7.     Given a shapely LineString which is assumed to be rectangular, return the
   
8.     line corresponding to a given side of the rectangle.
   
9.     """
   
10.     minx, miny, maxx, maxy = ls.bounds
    
11.     points = {'left': [(minx, miny), (minx, maxy)],
    
12.               'right': [(maxx, miny), (maxx, maxy)],
    
13.               'bottom': [(minx, miny), (maxx, miny)],
    
14.               'top': [(minx, maxy), (maxx, maxy)],}
    
15.     return sgeom.LineString(points[side])
    
16. 
    
17. # 在兰伯特投影的底部X轴上绘制刻度线
    
18. def lambert_xticks(ax, ticks):
    
19.     """Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
    
20.     te = lambda xy: xy[0]
    
21.     lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T
    
22.     xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te)
    
23.     ax.xaxis.tick_bottom()
    
24.     ax.set_xticks(xticks)
    
25.     ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels])
    
26. 
    
27. # 在兰伯特投影的左侧y轴上绘制刻度线
    
28. def lambert_yticks(ax, ticks):
    
29.     """Draw ricks on the left y-axis of a Lamber Conformal projection."""
    
30.     te = lambda xy: xy[1]
    
31.     lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
    
32.     yticks, yticklabels = _lambert_ticks(ax, ticks, 'left', lc, te)
    
33.     ax.yaxis.tick_left()
    
34.     ax.set_yticks(yticks)
    
35.     ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels])
    
36. 
    
37. # 获取兰伯特投影中底部X轴或左侧y轴的刻度线位置和标签
    
38. def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):
    
39.     """Get the tick locations and labels for an axis of a Lambert Conformal projection."""
    
40.     outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist())
    
41.     axis = find_side(outline_patch, tick_location)
    
42.     n_steps = 30
    
43.     extent = ax.get_extent(ccrs.PlateCarree())
    
44.     _ticks = []
    
45.     for t in ticks:
    
46.         xy = line_constructor(t, n_steps, extent)
    
47.         proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1])
    
48.         xyt = proj_xyz[..., :2]
    
49.         ls = sgeom.LineString(xyt.tolist())
    
50.         locs = axis.intersection(ls)
    
51.         if not locs:
    
52.             tick = [None]
    
53.         else:
    
54.             tick = tick_extractor(locs.xy)
    
55.         _ticks.append(tick[0])
    
56.     # Remove ticks that aren't visible:
    
57.     ticklabels = copy(ticks)
    
58.     while True:
    
59.         try:
    
60.             index = _ticks.index(None)
    
61.         except ValueError:
    
62.             break
    
63.         _ticks.pop(index)
    
64.         ticklabels.pop(index)
    
65.     return _ticks, ticklabels

复制代码

ps: 需要使用np.transpose()进行倒转维度，不然画不出来

1. #precipitationcal是指质控的降水数据
   
2. prep=ds.precipitationCal
   
3. lon=prep.lon
   
4. lat=prep.lat
   
5. prep=prep[0,:,:].transpose()
   
6. lon, lat = np.meshgrid(lon, lat)
   
7. # 设置地图投影和范围
   
8. proj = ccrs.LambertConformal(central_longitude=105, central_latitude=90, standard_parallels=(25, 47))
   
9. extent = [80, 130, 17, 55]
   
10. 
    
11. # 创建画布和子图
    
12. fig, ax = plt.subplots(figsize=(10,5), subplot_kw={'projection': proj})
    
13. 
    
14. # 设置网格点属性
    
15. #gl = ax.gridlines(draw_labels=True, alpha=0.5, linewidth=1, color='k', linestyle='--')
    
16. #gl.xlabels_top = False  # 关闭顶端标签
    
17. #gl.ylabels_right = False  # 关闭右侧标签
    
18. #gl.xformatter = LONGITUDE_FORMATTER  # x轴设为经度格式
    
19. #gl.yformatter = LATITUDE_FORMATTER  # y轴设为纬度格式
    
20. # 绘制降水数据，设置colorbar
    
21. levels = np.arange(0, 50, 2)
    
22. cf = ax.pcolormesh(lon, lat, prep, cmap=cmaps.WhiteBlue, vmin=0, vmax=50, transform=ccrs.PlateCarree())
    
23. plt.title('gpm precipitation at 20220903(UTC)', fontsize=12)
    
24. 
    
25. # 添加经纬度网格线
    
26. ax.gridlines(color='black', linestyle='dotted')
    
27. 
    
28. # 添加经纬度标签
    
29. xticks = list(np.arange(80,130,5))
    
30. yticks = list(np.arange(17,55,5))
    
31. fig.canvas.draw()
    
32. ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
    
33. ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
    
34. lambert_xticks(ax, xticks)
    
35. lambert_yticks(ax, yticks)
    
36. 
    
37. # 添加colorbar
    
38. cbar_ax = fig.add_axes([0.80, 0.12, 0.02, 0.8])
    
39. cb = plt.colorbar(cf, cax=cbar_ax, ticks=levels)
    
40. cb.set_label('Total Precipitaion (mm)', fontdict={'size':10})
    
41. cb.ax.tick_params(labelsize=10)
    
42. # 添加省界数据
    
43. province = shpreader.Reader('/home/mw/project/cnhimap.shp')
    
44. ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5, edgecolor='k', facecolor='none')
    
45. # 设置地图范围和显示
    
46. ax.set_extent(extent, crs=ccrs.PlateCarree())
    
47. plt.show()

复制代码

                               
登录/注册后可看大图

1.3 GMI(hdf5)¶

1. # 读取GPM GMI降水产品
   
2. import numpy as np
   
3. import h5py                  
   
4. import matplotlib.pyplot as plt
   
5. import cartopy.crs as ccrs
   
6. import warnings
   
7. warnings.filterwarnings("ignore")
   
8. file = '/home/mw/input/GPM5633/gpmPython/2A.GPM.GMI.GPROF.20210829-S151345-E151504.042624.V05B.subset.HDF5'
   
9. f = h5py.File(file,"r")
   
10. print(f['S1'].keys())
    
11. prns = f['/S1/surfacePrecipitation']
    
12. lon  = np.array(f["/S1/Longitude"])
    
13. lat  = np.array(f["/S1/Latitude"])
    
14. print('precipRate dimension sizes are:',prns.shape)
    
15. 
    
16. prnsdata = np.ndarray(shape=prns.shape,dtype=float)
    
17. prns.read_direct(prnsdata)
    
18. prnsdata[prnsdata <= -9999] = np.nan
    
19. 
    
20. start = 0
    
21. end = 43
    
22. mysub = prnsdata[start:end,:]
    
23. mylon = lon[start:end,:]
    
24. mylat = lat[start:end,:]
    
25. prnsdata.min()
    
26. 
    
27. # Draw the subset of near-surface precipitation rate
    
28. fig = plt.figure(figsize=(8, 8))
    
29. ax = plt.axes(projection=ccrs.PlateCarree())
    
30. plt.title('test')
    
31. # Add coastlines and gridlines
    
32. ax.coastlines(resolution="50m",linewidth=0.5)
    
33. gl = ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,
    
34.                   linewidth=0.8,color='#555555',alpha=0.5,linestyle='--')
    
35. # Axis labels
    
36. gl.xlabels_top = False
    
37. gl.ylabels_right = False
    
38. gl.xlines = True
    
39. 
    
40. # Plot the scatter diagram
    
41. pp = plt.scatter(mylon, mylat, c=mysub, cmap=cmap_Rr, transform=ccrs.PlateCarree())
    
42. 
    
43. # Add a colorbar to the bottom of the plot.
    
44. fig.subplots_adjust(bottom=0.15,left=0.06,right=0.94)
    
45. cbar_ax = fig.add_axes([0.12, 0.08, 0.76, 0.025])  
    
46. cbar = plt.colorbar(pp,cax=cbar_ax,orientation='horizontal')
    
47. 
    

复制代码

                               
登录/注册后可看大图

1.4 GMI可视化nc版

1. ds1 = xr.open_dataset('/home/mw/project/2A.GPM.GMI.GPROF2021v1.20220902-S223342-E000616.048370.V07A.HDF5.nc4')
   
2. print(ds1)
   
3. #看文件描述的近地面降水 GMI
   
4. from pathlib import Path
   
5. from collections import namedtuple
   
6. 
   
7. import h5py
   
8. import numpy as np
   
9. import matplotlib.pyplot as plt
   
10. import matplotlib.ticker as mticker
    
11. import matplotlib.cm as cm
    
12. import matplotlib.colors as mcolors
    
13. import cmaps
    
14. import cartopy.crs as ccrs
    
15. import cartopy.feature as cfeature
    
16. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
    
17. 
    
18. def make_Rr_cmap(levels):
    
19.     '''制作降水的colormap.'''
    
20.     nbin = len(levels) - 1
    
21.     cmap = cm.get_cmap('jet', nbin)
    
22.     norm = mcolors.BoundaryNorm(levels, nbin)
    
23. 
    
24.     return cmap, norm
    
25. 
    
26. def region_mask(lon, lat, extents):
    
27.     '''返回用于索引经纬度方框内数据的Boolean数组.'''
    
28.     lonmin, lonmax, latmin, latmax = extent
    
29.     return (
    
30.         (lon >= lonmin) & (lon <= lonmax) &
    
31.         (lat >= latmin) & (lat <= latmax)
    
32.     )
    
33. 
    
34. extent = [80, 130, 17, 55]
    
35. surfacePrecipitation = ds1.S1_surfacePrecipitation
    
36. lon=surfacePrecipitation.S1_Longitude
    
37. lat=surfacePrecipitation.S1_Latitude
    
38. #prep=prep[0,:,:].transpose()
    
39. levels_Rr = [0.1, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 15.0, 20.0, 25, 30]
    
40. 
    
41. cmap_Rr,norm_Rr = make_Rr_cmap(levels_Rr)
    
42. #cmap_LH, norm_LH = make_LH_cmap(levels_LH)
    
43. 
    
44. #lon, lat = np.meshgrid(lon, lat)
    
45. nscan, npixel = lon.shape
    
46. midpixel = npixel // 2
    
47. mask = region_mask(lon[:, midpixel], lat[:, midpixel], extent)
    
48. index = np.s_[mask, :]
    
49. lon = lon[index]
    
50. lat = lat[index]
    
51. surfacePrecipitation = surfacePrecipitation[index]
    
52.     # 画出GMI降水.
    
53. 
    
54. proj = ccrs.PlateCarree()
    
55. fig = plt.figure(figsize=(10, 8))
    
56. ax = fig.add_subplot(111, projection=proj)
    
57. ax.coastlines(resolution='50m', lw=0.5)
    
58. ax.add_feature(cfeature.OCEAN.with_scale('50m'))
    
59. ax.add_feature(cfeature.LAND.with_scale('50m'))
    
60. ax.set_xticks(np.arange(-180, 181, 5), crs=proj)
    
61. ax.set_yticks(np.arange(-90, 91, 5), crs=proj)
    
62. ax.xaxis.set_minor_locator(mticker.AutoMinorLocator(2))
    
63. ax.yaxis.set_minor_locator(mticker.AutoMinorLocator(2))
    
64. ax.xaxis.set_major_formatter(LongitudeFormatter())
    
65. ax.yaxis.set_major_formatter(LatitudeFormatter())
    
66. ax.set_extent(extent, crs=proj)
    
67. ax.tick_params(labelsize='large')
    
68. cmap_Rr.set_under(color='lavender')
    
69. 
    
70. im = ax.contourf(
    
71.     lon, lat, surfacePrecipitation, levels_Rr,
    
72.     cmap=cmap_Rr, norm=norm_Rr, alpha=0.5,
    
73.     extend='both', transform=proj
    
74. )
    
75. cbar = fig.colorbar(im, ax=ax, ticks=levels_Rr)
    
76. cbar.set_label('Rain Rate (mm/hr)', fontsize='large')
    
77. cbar.ax.tick_params(labelsize='large')
    
78. plt.show()

复制代码

                               
登录/注册后可看大图

这画图范围是我一个个区域去试出来，感觉不太方便，不知有没更高效的方法

数据和程序都放在和鲸社区，想复现去帖子顶端的链接，Python环境问题的话自己查一下