GPM卫星数据hdf5格式读取与绘图
数据说明：GPM的DPR降水产品与SLH潜热产品（hdf5格式）
1、导入python库和hdf5文件结构探索
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
f = h5py.File('/home/mw/input/GPM5859/2A.GPM.DPR.GPM-SLH.20210929-S142530-E155804.043106.V07A.HDF5', 'r')
你刚开始拿到数据多半不知怎么看结构，一定很疑惑f['Swath/latentHeating'][:]怎么来的
hdf5数据逻辑和nc不太一样，
且看我下面如何操作
#查看f下目录
for key in f.keys():
print(key)
AlgorithmRuntimeInfo
Swath
好，我们继续看那个Swath下面的目录
print(f['Swath'].keys())
<KeysViewHDF5 ['ScanTime', 'Latitude', 'Longitude', 'sunLocalTime', 'latentHeating', 'Q1minusQR', 'Q2', 'rainTypeSLH', 'stormTopHeight', 'nearMeltLevel', 'nearSurfLevel', 'topoLevel', 'meltLevel', 'levelConvUpper', 'nearSurfacePrecipRate', 'precipRateNearMelt', 'precipRateConvUpper', 'rainType2ADPR', 'surfaceType']>
或者直接用一下函数
/
//AlgorithmRuntimeInfo
//Swath
//Swath/ScanTime
//Swath/ScanTime/Year
//Swath/ScanTime/Month
//Swath/ScanTime/DayOfMonth
//Swath/ScanTime/Hour
//Swath/ScanTime/Minute
//Swath/ScanTime/Second
//Swath/ScanTime/MilliSecond
//Swath/ScanTime/DayOfYear
//Swath/ScanTime/SecondOfDay
//Swath/Latitude
//Swath/Longitude
//Swath/sunLocalTime
//Swath/latentHeating
//Swath/Q1minusQR
//Swath/Q2
//Swath/rainTypeSLH
//Swath/stormTopHeight
//Swath/nearMeltLevel
//Swath/nearSurfLevel
//Swath/topoLevel
//Swath/meltLevel
//Swath/levelConvUpper
//Swath/nearSurfacePrecipRate
//Swath/precipRateNearMelt
//Swath/precipRateConvUpper
//Swath/rainType2ADPR
//Swath/surfaceType
2、绘图降水图（二级降水数据）
f2 = h5py.File('/home/mw/input/GPM5859/2A.GPM.DPR.V9-20211125.20210929-S142530-E155804.043106.V07A.HDF5', 'r')
def print_hdf5(name, obj):
print(name)
if isinstance(obj, h5py.Group):
for key in obj.keys():
subname = '{}/{}'.format(name, key)
print_hdf5(subname, obj)
print_hdf5('/', f2)输出目录太多，姑且隐藏了，可fork查看
2.1 普通版
In :
import h5py
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
extents =
#读取变量
f2 = h5py.File('/home/mw/input/GPM5859/2A.GPM.DPR.V9-20211125.20210929-S142530-E155804.043106.V07A.HDF5', 'r')
lon = f2['FS/Longitude'][:]# FS: nscan: 7935, nray: 49 HS: nscan: 7935, nrayHS: 24
lat = f2['FS/Latitude'][:]
precipRateNearSurface = f2['FS/SLV/precipRateNearSurface'][:]# 近地面降水率
def region_mask(lon, lat, extents):
minlon, maxlon, minlat, maxlat = extents
return (lon > minlon) & (lon < maxlon) & (lat > minlat) & (lat < maxlat)
#获取了lon形状,nscan为轨道数,nray为射束数。
nscan, nray = lon.shape
#计算射束的中间索引midray。
midray = nray // 2
#调用region_mask函数,以中间射束的经纬度数据为基准,利用给定的经纬度范围extents生成筛选mask
mask = region_mask(lon[:, midray], lat[:, midray], extents)
#用np.s_[]将mask封装成切片索引,方便对多个变量切片
index = np.s_
#应用切片索引,对lon变量进行切片,提取筛选后的子集
lon = lon
#同上
lat = lat
#同上
precipRateNearSurface = precipRateNearSurface
#对缺测值-9999处理
precipRateNearSurface = np.nan
#画图
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(20, 12))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
im = ax.contourf(lon, lat, precipRateNearSurface,cmaps=plt.cm.Spectral_r)
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)
plt.colorbar(im, extend='both')
plt.title('PS DPR')
/opt/conda/lib/python3.7/site-packages/cartopy/mpl/geoaxes.py:1508: UserWarning: The following kwargs were not used by contour: 'cmaps'
result = matplotlib.axes.Axes.contourf(self, *args, **kwargs)
Text(0.5, 1.0, 'PS DPR')
2.2 封装优化版
In :
import h5py
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
def get_extents():
return
def region_mask(lon, lat, extents):# 筛选目标区域
'''返回用于索引经纬度方框内数据的Boolean数组.'''
lonmin, lonmax, latmin, latmax = extents
return (
(lon >= lonmin) & (lon <= lonmax) &
(lat >= latmin) & (lat <= latmax)
)
def read_data(filename):
with h5py.File(filename, 'r') as f:
precipRateNearSurface = f['FS/SLV/precipRateNearSurface'][:]
lon = f['FS/Longitude'][:]
lat = f['FS/Latitude'][:]
return precipRateNearSurface, lon, lat
def process_missing(precipRateNearSurface):
precipRateNearSurface = np.nan
return precipRateNearSurface
def plot_latent_heating(lon, lat, var, extents):
fig = plt.figure(figsize=(20, 12))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
contourf = ax.contourf(lon, lat, var)
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=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 5), crs=ccrs.PlateCarree())
ax.xaxis.set_minor_locator(plt.MultipleLocator(2))
ax.yaxis.set_minor_locator(plt.MultipleLocator(2))
ax.set_extent(extents)
plt.colorbar(contourf, extend='both')
plt.title('precipRateNearSurface DPR')
plt.show()
# 读取数据并处理缺失
filename = '/home/mw/input/GPM5859/2A.GPM.DPR.V9-20211125.20210929-S142530-E155804.043106.V07A.HDF5'
precipRateNearSurface, lon, lat = read_data(filename)
precipRateNearSurface = process_missing(precipRateNearSurface)
print(lon.shape,lat.shape,precipRateNearSurface.shape)
# 预处理区域mask
extents = get_extents()
nscan, nray = precipRateNearSurface.shape
print(nscan, nray)
midray = nray // 2
#mask = (lon[:, midray] > extents) & (lon[:, midray] < extents) & (lat[:, midray] > extents) & (lat[:, midray] < extents)
mask = region_mask(lon[:, midray], lat[:, midray], extents)
index = np.s_
lon = lon
lat = lat
precipRateNearSurface = precipRateNearSurface
# 绘图
plot_latent_heating(lon, lat, precipRateNearSurface, extents)
(7935, 49) (7935, 49) (7935, 49)
7935 49
3、二级潜热图
3.1 普通版
In :
import h5py
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
extents =
#读取变量
f = h5py.File('/home/mw/input/GPM5859/2A.GPM.DPR.GPM-SLH.20210929-S142530-E155804.043106.V07A.HDF5', 'r')
latentHeating_SLH = f['Swath/latentHeating'][:]
lon = f['Swath/Longitude'][:]
lat = f['Swath/Latitude'][:]
def region_mask(lon, lat, extents):
minlon, maxlon, minlat, maxlat = extents
return (lon > minlon) & (lon < maxlon) & (lat > minlat) & (lat < maxlat)
#获取了lon形状,nscan为轨道数,nray为射束数。
nscan, nray = lon.shape
#计算射束的中间索引midray。
midray = nray // 2
#调用region_mask函数,以中间射束的经纬度数据为基准,利用给定的经纬度范围extents生成筛选mask
mask = region_mask(lon[:, midray], lat[:, midray], extents)
#用np.s_[]将mask封装成切片索引,方便对多个变量切片
index = np.s_
#应用切片索引,对lon变量进行切片,提取筛选后的子集
lon = lon
#同上
lat = lat
#同上
latentHeating_SLH = latentHeating_SLH
#对缺测值-9999处理
latentHeating_SLH = np.nan
#画图
fig = plt.figure(figsize=(20, 12))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
im = ax.contourf(lon, lat, latentHeating_SLH[:,:,10],)
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)
plt.colorbar(im, extend='both')
plt.title('Latent Heating SLH')
Text(0.5, 1.0, 'Latent Heating SLH')
3.2 封装优化版
In :
import h5py
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
def get_extents():
return
def read_data(filename):
with h5py.File(filename, 'r') as f:
latentHeating_SLH = f['Swath/latentHeating'][:]
lon = f['Swath/Longitude'][:]
lat = f['Swath/Latitude'][:]
return latentHeating_SLH, lon, lat
def process_missing(latentHeating_SLH):
latentHeating_SLH = np.nan
return latentHeating_SLH
def plot_latent_heating(lon, lat, latentHeating_SLH, extents):
fig = plt.figure(figsize=(20, 12))
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree())
contourf = ax.contourf(lon, lat, latentHeating_SLH[:,:,10])
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=ccrs.PlateCarree())
ax.set_yticks(np.arange(-90, 91, 5), crs=ccrs.PlateCarree())
ax.xaxis.set_minor_locator(plt.MultipleLocator(2))
ax.yaxis.set_minor_locator(plt.MultipleLocator(2))
ax.set_extent(extents)
plt.colorbar(contourf, extend='both')
plt.title('Latent Heating SLH')
plt.show()
# 读取数据并处理缺失
filename = '/home/mw/input/GPM5859/2A.GPM.DPR.GPM-SLH.20210929-S142530-E155804.043106.V07A.HDF5'
latentHeating_SLH, lon, lat = read_data(filename)
latentHeating_SLH = process_missing(latentHeating_SLH)
# 预处理区域mask
extents = get_extents()
nscan, nray, _ = latentHeating_SLH.shape
print(nscan, nray)
midray = nray // 2
mask = (lon[:, midray] > extents) & (lon[:, midray] < extents) & (lat[:, midray] > extents) & (lat[:, midray] < extents)
#index = np.where(mask)
index = np.s_
lon = lon
lat = lat
latentHeating_SLH = latentHeating_SLH
# 绘图
plot_latent_heating(lon, lat, latentHeating_SLH, extents)
7935 49
4、全球DPR降水示例图及局部
4.1 全球示例图
(7933, 49) (7933, 49) (7933, 49)
4.2 中国附近区域可视化
(7933, 49) (7933, 49) (7933, 49)
4.3 降水剖面图
4.3.1 y方向
(7933, 49) (7933, 49) (7933, 49, 176)
/opt/conda/lib/python3.7/site-packages/cmaps/cmaps.py:2173: UserWarning: Trying to register the cmap 'NCV_jaisnd' which already exists.
matplotlib.cm.register_cmap(name=cname, cmap=cmap)
4.3.2 x方向
(7933, 49) (7933, 49) (7933, 49, 176)
4.4 潜热剖面图
4.4.1 y方向
(7935, 49) (7935, 49) (7935, 49, 80)
4.4.2 x方向
(7935, 49) (7935, 49) (7935, 49, 80)
小结
1、画个全球图有助于理解卫星数据分布
2、代码进行封装后复用更方便
3、本文的剖面是沿着x轴与y轴的，可以研究如何自定义剖线起点与终点进行剖面
4、看着剖面的横纵坐标很难受吧，研究GPM中对应高度的参数并将剖面的横纵坐标进行修正
PS: 勇敢的少年，快去创造奇迹
完整代码与文件在此
文章来源于微信公众号：气ython风雨
页:
[1]