python读取GPM二级降水等卫星数据及思路

import xarray as xr
from collections import namedtuple
import numpy as np
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.ticker as mticker
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as mcolors


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


Point = namedtuple('Point', ['x', 'y'])  # 命名元祖, 等号左边Point是对象， 里面属性是x, y 赋予每一个Point对象x和y的属性
Pair = namedtuple('Pair', ['start', 'end'])  # 命名元祖， 等号左边Pair是对象， 里面属性是x, y 赋予每一个Point对象x和y的属性
# 读取文件，习惯使用nc格式，HDF格式修改对应变量表示，一样可以使用，变量通过variables查看
# 时间字符串，表示所取文件时间
time = '2021-09-29'
filepath_GMI = r"E:\GPM\代码例子\2A-CLIM.GPM.GMI.GPROF2021v1.20210929-S142530-E155804.043106.V07A.HDF5.nc4"
filepath_DPR = r"E:\GPM\代码例子\2A.GPM.DPR.V9-20211125.20210929-S142530-E155804.043106.V07A.HDF5.nc4"


# 地图范围, 最后画图呈现的区域。 表示minlon, maxlon, minlat, maxlat
extents = [127, 147, 17, 37]

# 读取物理量
with xr.open_dataset(filepath_GMI) as f:  # 读取GMI数据.
    lon_GMI = f['S1_Longitude'][:]  # nscan: 2963, npixel: 221，默认值221
    lat_GMI = f['S1_Latitude'][:]  # nscan: 2963, npixel: 221, 经纬度一致
    surfacePrecipitation = f['S1_surfacePrecipitation'][:]


with xr.open_dataset(filepath_DPR) as f:  # 读取DPR数据.
    lon_DPR = f['FS_Longitude'][:]  # FS: nscan: 7935, nray: 49 HS: nscan: 7935, nrayHS: 24
    lat_DPR = f['FS_Latitude'][:]
    precipRateNearSurface = f['FS_SLV_precipRateNearSurface'][:]  # 近地面降水率

# 截取GMI数据.
# 原来经度: (2963, 221), 纬度： (2963, 221)，现在经度: (187, 221), 纬度： (187, 221)
nscan, npixel = lon_GMI.shape
midpixel = npixel // 2  # 中间是110数值
# 将lon_GMI， lat_GMI放到目标区域，定一维， 二维筛选不出，取中间轨道，筛选其符合的经纬度, 取其他与中心轨道区域相平行的区域，显示出来为卫星扫描轨道
mask = region_mask(lon_GMI[:, midpixel], lat_GMI[:, midpixel], extents)
index = np.s_[mask]  # 转化格式，方便后面切片 np.s_[mask, :],第二位置表示后面筛选几个轨道，小于221
lon_GMI = lon_GMI[index]  # (187, 221)
lat_GMI = lat_GMI[index]  # 筛选出目标区域
surfacePrecipitation = surfacePrecipitation[index]

# 截取DPR数据.
nscan, nray = lon_DPR.shape
midray = nray // 2
mask = region_mask(lon_DPR[:, midray], lat_DPR[:, midray], extents)
index = np.s_[mask]
lon_DPR = lon_DPR[index]
lat_DPR = lat_DPR[index]
precipRateNearSurface = precipRateNearSurface[index]

# 设置缺测. 注意到GMI的缺测是-9999,缺测值看文件说明
for data in [
    surfacePrecipitation, precipRateNearSurface,
]:
   data.values[data <= -9999] = np.nan

# 画出地图.
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(extents, crs=proj)
ax.tick_params(labelsize='large')

# 设置颜色标签
def make_Rr_cmap(levels):
    '''制作降水的colormap.'''
    nbin = len(levels) - 1
    cmap = cm.get_cmap('jet', nbin)
    norm = mcolors.BoundaryNorm(levels, nbin)
    return cmap, norm

# 设置level值
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)
# 画出GMI降水.
cmap_Rr.set_under(color='lavender')  # lavender
im = ax.contourf(
    lon_GMI, lat_GMI, surfacePrecipitation,  # 三个物理量的shape都是(187, 221)就是在187*221的格点上赋予这三个物理量
    levels_Rr,  cmap=cmap_Rr, norm=norm_Rr, alpha=0.5,  # alpha透明度
    extend='both', transform=proj
)

# 画出DPR降水.
im = ax.contourf(
    lon_DPR, lat_DPR, precipRateNearSurface, levels_Rr,  # 三个物理量为 (500, 49)
    cmap=cmap_Rr, norm=norm_Rr, 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')
# 设置标题.
ax.set_title(f'GMI and DPR Rain Rate on {time}', fontsize='x-large')
plt.show()

卫星扫描数据与网格数据最大不同，是在数据处理方面。

对于代码结果如图所示，之前网上找了很久没找到，在气象家园看到了大佬写的。气象家园链接如下：Python绘制GPM DPR二级降水和潜热-编程作图-气象家园_气象人自己的家园 (06climate.com)

有些人可能没注册，无法访问。本代码进行了详细注释。文章数据获取链接如下：

GES DISC Search: Showing 1 - 25 of 1747 datasets (nasa.gov)

如果不想下载可私信我，我将文件发送给您。

文章以转载形式表达对原文的敬意。