GPM卫星 IMERG 降水产品的简单可视化

# 导入库
import numpy as np
import xarray as xr
#可视化
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeat
import cartopy.io.shapereader as shpreader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cmaps
from shapely import geometry as sgeom
# 辅助模块
from glob import glob
from copy import copy

Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1

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

'''
参考Andrew Dawson 提供的解决方法: https://gist.github.com/ajdawson/dd536f786741e987ae4e
'''
# 给出一个假定为矩形的LineString，返回对应于矩形给定边的直线。
def find_side(ls, side):
    """
    Given a shapely LineString which is assumed to be rectangular, return the
    line corresponding to a given side of the rectangle.
    """
    minx, miny, maxx, maxy = ls.bounds
    points = {'left': [(minx, miny), (minx, maxy)],
              'right': [(maxx, miny), (maxx, maxy)],
              'bottom': [(minx, miny), (maxx, miny)],
              'top': [(minx, maxy), (maxx, maxy)],}
    return sgeom.LineString(points[side])

# 在兰伯特投影的底部X轴上绘制刻度线
def lambert_xticks(ax, ticks):
    """Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
    te = lambda xy: xy[0]
    lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T
    xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te)
    ax.xaxis.tick_bottom()
    ax.set_xticks(xticks)
    ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels])

# 在兰伯特投影的左侧y轴上绘制刻度线
def lambert_yticks(ax, ticks):
    """Draw ricks on the left y-axis of a Lamber Conformal projection."""
    te = lambda xy: xy[1]
    lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
    yticks, yticklabels = _lambert_ticks(ax, ticks, 'left', lc, te)
    ax.yaxis.tick_left()
    ax.set_yticks(yticks)
    ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels])

# 获取兰伯特投影中底部X轴或左侧y轴的刻度线位置和标签
def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):
    """Get the tick locations and labels for an axis of a Lambert Conformal projection."""
    outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist())
    axis = find_side(outline_patch, tick_location)
    n_steps = 30
    extent = ax.get_extent(ccrs.PlateCarree())
    _ticks = []
    for t in ticks:
        xy = line_constructor(t, n_steps, extent)
        proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1])
        xyt = proj_xyz[..., :2]
        ls = sgeom.LineString(xyt.tolist())
        locs = axis.intersection(ls)
        if not locs:
            tick = [None]
        else:
            tick = tick_extractor(locs.xy)
        _ticks.append(tick[0])
    # Remove ticks that aren't visible: 
    ticklabels = copy(ticks)
    while True:
        try:
            index = _ticks.index(None)
        except ValueError:
            break
        _ticks.pop(index)
        ticklabels.pop(index)
    return _ticks, ticklabels

#precipitationcal是指质控的降水数据 
prep=ds.precipitationCal 
lon=prep.lon
lat=prep.lat
prep=prep[0,:,:].transpose()
lon, lat = np.meshgrid(lon, lat)
# 设置地图投影和范围
proj = ccrs.LambertConformal(central_longitude=105, central_latitude=90, standard_parallels=(25, 47))
extent = [80, 130, 17, 55]

# 创建画布和子图
fig, ax = plt.subplots(figsize=(10,5), subplot_kw={'projection': proj})

# 设置网格点属性
#gl = ax.gridlines(draw_labels=True, alpha=0.5, linewidth=1, color='k', linestyle='--')
#gl.xlabels_top = False  # 关闭顶端标签
#gl.ylabels_right = False  # 关闭右侧标签
#gl.xformatter = LONGITUDE_FORMATTER  # x轴设为经度格式
#gl.yformatter = LATITUDE_FORMATTER  # y轴设为纬度格式
# 绘制降水数据，设置colorbar
levels = np.arange(0, 50, 2)
cf = ax.pcolormesh(lon, lat, prep, cmap=cmaps.WhiteBlue, vmin=0, vmax=50, transform=ccrs.PlateCarree())
plt.title('gpm precipitation at 20220903(UTC)', fontsize=12)

# 添加经纬度网格线
ax.gridlines(color='black', linestyle='dotted')

# 添加经纬度标签
xticks = list(np.arange(80,130,5))
yticks = list(np.arange(17,55,5))
fig.canvas.draw()
ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER) 
ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lambert_xticks(ax, xticks)
lambert_yticks(ax, yticks)

# 添加colorbar
cbar_ax = fig.add_axes([0.80, 0.12, 0.02, 0.8])
cb = plt.colorbar(cf, cax=cbar_ax, ticks=levels)
cb.set_label('Total Precipitaion (mm)', fontdict={'size':10})
cb.ax.tick_params(labelsize=10)
# 添加省界数据
province = shpreader.Reader('/home/mw/project/cnhimap.shp')
ax.add_geometries(province.geometries(), crs=ccrs.PlateCarree(), linewidths=0.5, edgecolor='k', facecolor='none')
# 设置地图范围和显示
ax.set_extent(extent, crs=ccrs.PlateCarree())
plt.show()

/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:40: DeprecationWarning: The outline_patch property is deprecated. Use GeoAxes.spines['geo'] or the default Axes properties instead.

# 读取GPM GMI降水产品
import numpy as np
import h5py                  
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import warnings
warnings.filterwarnings("ignore")
file = '/home/mw/input/GPM5633/gpmPython/2A.GPM.GMI.GPROF.20210829-S151345-E151504.042624.V05B.subset.HDF5'
f = h5py.File(file,"r")
print(f['S1'].keys())
prns = f['/S1/surfacePrecipitation']
lon  = np.array(f["/S1/Longitude"])
lat  = np.array(f["/S1/Latitude"])
print('precipRate dimension sizes are:',prns.shape)

<KeysViewHDF5 ['ScanTime', 'Latitude', 'Longitude', 'probabilityOfPrecip', 'surfacePrecipitation', 'frozenPrecipitation', 'convectivePrecipitation']>
precipRate dimension sizes are: (43, 221)

lon.shape

(43, 221)

prnsdata = np.ndarray(shape=prns.shape,dtype=float)
prns.read_direct(prnsdata)
prnsdata[prnsdata <= -9999] = np.nan

start = 0
end = 43
mysub = prnsdata[start:end,:]
mylon = lon[start:end,:]
mylat = lat[start:end,:]
prnsdata.min()

nan

# Draw the subset of near-surface precipitation rate 
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
plt.title('test')
# Add coastlines and gridlines
ax.coastlines(resolution="50m",linewidth=0.5)
gl = ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,
                  linewidth=0.8,color='#555555',alpha=0.5,linestyle='--')
# Axis labels
gl.xlabels_top = False
gl.ylabels_right = False
gl.xlines = True

# Plot the scatter diagram 
pp = plt.scatter(mylon, mylat, c=mysub, cmap=cmap_Rr, transform=ccrs.PlateCarree())

# Add a colorbar to the bottom of the plot.
fig.subplots_adjust(bottom=0.15,left=0.06,right=0.94)
cbar_ax = fig.add_axes([0.12, 0.08, 0.76, 0.025])  
cbar = plt.colorbar(pp,cax=cbar_ax,orientation='horizontal')

ds1 = xr.open_dataset('/home/mw/project/2A.GPM.GMI.GPROF2021v1.20220902-S223342-E000616.048370.V07A.HDF5.nc4')
print(ds1)
#看文件描述的近地面降水 GMI
from pathlib import Path
from collections import namedtuple

import h5py
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.cm as cm
import matplotlib.colors as mcolors
import cmaps
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

def make_Rr_cmap(levels):
    '''制作降水的colormap.'''
    nbin = len(levels) - 1
    cmap = cm.get_cmap('jet', nbin)
    norm = mcolors.BoundaryNorm(levels, nbin)

    return cmap, norm

def region_mask(lon, lat, extents):
    '''返回用于索引经纬度方框内数据的Boolean数组.'''
    lonmin, lonmax, latmin, latmax = extent
    return (
        (lon >= lonmin) & (lon <= lonmax) &
        (lat >= latmin) & (lat <= latmax)
    )

extent = [80, 130, 17, 55]
surfacePrecipitation = ds1.S1_surfacePrecipitation
lon=surfacePrecipitation.S1_Longitude
lat=surfacePrecipitation.S1_Latitude
#prep=prep[0,:,:].transpose()
levels_Rr = [0.1, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 15.0, 20.0, 25, 30]

cmap_Rr,norm_Rr = make_Rr_cmap(levels_Rr)
#cmap_LH, norm_LH = make_LH_cmap(levels_LH)

#lon, lat = np.meshgrid(lon, lat)
nscan, npixel = lon.shape
midpixel = npixel // 2
mask = region_mask(lon[:, midpixel], lat[:, midpixel], extent)
index = np.s_[mask, :]
lon = lon[index]
lat = lat[index]
surfacePrecipitation = surfacePrecipitation[index]
    # 画出GMI降水.

proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection=proj)
ax.coastlines(resolution='50m', lw=0.5)
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.set_xticks(np.arange(-180, 181, 5), crs=proj)
ax.set_yticks(np.arange(-90, 91, 5), crs=proj)
ax.xaxis.set_minor_locator(mticker.AutoMinorLocator(2))
ax.yaxis.set_minor_locator(mticker.AutoMinorLocator(2))
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.set_extent(extent, crs=proj)
ax.tick_params(labelsize='large')
cmap_Rr.set_under(color='lavender')

im = ax.contourf(
    lon, lat, surfacePrecipitation, levels_Rr,
    cmap=cmap_Rr, norm=norm_Rr, alpha=0.5,
    extend='both', transform=proj
)
cbar = fig.colorbar(im, ax=ax, ticks=levels_Rr)
cbar.set_label('Rain Rate (mm/hr)', fontsize='large')
cbar.ax.tick_params(labelsize='large')
plt.show()

<xarray.Dataset>
Dimensions:                  (nscan: 2963, npixel: 221)
Coordinates:
    S1_Longitude             (nscan, npixel) float32 ...
    S1_Latitude              (nscan, npixel) float32 ...
Dimensions without coordinates: nscan, npixel
Data variables:
    S1_surfacePrecipitation  (nscan, npixel) float32 ...
Attributes:
    FileHeader:                      DOI=10.5067/GPM/GMI/GPM/GPROF/2A/07;\nDO...
    InputRecord:                     InputFileNames=1C-R.GPM.GMI.XCAL2016-C.2...
    NavigationRecord:                LongitudeOnEquator=-41.792646;\nUTCDateT...
    FileInfo:                        DataFormatVersion=7a;\nTKCodeBuildVersio...
    GprofInfo:                       Satellite=GPM;\nSensor=GMI;\nPreProcesso...
    S1.SwathHeader:                  NumberScansInSet=1;\nMaximumNumberScansT...
    S1.fullnamepath:                 /S1
    DODS_EXTRA.Unlimited_Dimension:  nscan
    history:                         2023-05-25 11:10:38 GMT Hyrax-1.16.3 htt...