数据说明:GPM的DPR降水产品与SLH潜热产品(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
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[key])
print_hdf5('/', f2)
输出目录太多,姑且隐藏了,可fork查看
In [6]:
import h5py
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
extents = [133, 140, 20, 30]
#读取变量
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_[mask]
#应用切片索引,对lon变量进行切片,提取筛选后的子集
lon = lon[index]
#同上
lat = lat[index]
#同上
precipRateNearSurface = precipRateNearSurface[index]
#对缺测值-9999处理
precipRateNearSurface[precipRateNearSurface <= -9999] = 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')
In [7]:
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 [127, 147, 17, 37]
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[precipRateNearSurface <= -9999] = 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[0]) & (lon[:, midray] < extents[1]) & (lat[:, midray] > extents[2]) & (lat[:, midray] < extents[3])
mask = region_mask(lon[:, midray], lat[:, midray], extents)
index = np.s_[mask]
lon = lon[index]
lat = lat[index]
precipRateNearSurface = precipRateNearSurface[index]
# 绘图
plot_latent_heating(lon, lat, precipRateNearSurface, extents)
(7935, 49) (7935, 49) (7935, 49)
7935 49
In [8]:
import h5py
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
extents = [133, 140, 20, 30]
#读取变量
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_[mask]
#应用切片索引,对lon变量进行切片,提取筛选后的子集
lon = lon[index]
#同上
lat = lat[index]
#同上
latentHeating_SLH = latentHeating_SLH[index]
#对缺测值-9999处理
latentHeating_SLH[latentHeating_SLH <= -9999] = 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')
In [9]:
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 [127, 147, 17, 37]
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[latentHeating_SLH <= -9999] = 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[0]) & (lon[:, midray] < extents[1]) & (lat[:, midray] > extents[2]) & (lat[:, midray] < extents[3])
#index = np.where(mask)[0]
index = np.s_[mask]
lon = lon[index]
lat = lat[index]
latentHeating_SLH = latentHeating_SLH[index]
# 绘图
plot_latent_heating(lon, lat, latentHeating_SLH, extents)
7935 49
(7933, 49) (7933, 49) (7933, 49)
(7933, 49) (7933, 49) (7933, 49)
(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)
(7933, 49) (7933, 49) (7933, 49, 176)
(7935, 49) (7935, 49) (7935, 49, 80)
(7935, 49) (7935, 49) (7935, 49, 80)
PS: 勇敢的少年,快去创造奇迹