Python可视化DPR（hdf5格式）降水

 前言：最近大创开始学习卫星数据了，建议先学习hdf5格式再来学习数据画图~

参考网页：Learning HDF5

（可以下载user guide）

本文使用的数据是CMB数据，时间是2014年7月3日

废话不多说，上代码：

import h5py
from collections import namedtuple
import numpy as np
import matplotlib.cm as cm
import matplotlib.colors as mcolors
import matplotlib.ticker as mticker
import matplotlib.pyplot as plt
import cmaps
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import pandas as pd
import warnings
warnings.filterwarnings('ignore')


def region_mask(vars, extents):
    '''返回切割好的经纬度和降水信息
    vars:带切割变量
    extents:画图经纬范围
    此方法是根据扫描中心线判断范围
    若中心线不在范围内，失败
    '''
    lonmin, lonmax, latmin, latmax = extents
    nscan, nray = vars[0].shape
    midray = nray // 2
    mask = (vars[0][:, midray] >= lonmin) & (vars[0][:, midray] <= lonmax) & (vars[1][:, midray] >= latmin) & (vars[1][:, midray] <= latmax)
    x = pd.DataFrame(mask)
    x.columns = ['Bool']
    var_new = []
    for i in range(len(vars)):
        var = pd.DataFrame(vars[i])
        var_s = pd.merge(x,var,left_index=True,right_index=True)
        var_new.append(var_s.loc[var_s['Bool']==True,~var_s.columns.isin(['Bool'])])
    return var_new


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


#28.5 108.5中心，5°*5°大小


path = r"\HDF5学习\1955.HDF5"
with h5py.File(str(path), 'r') as f:
    lon = f['NS/Longitude'][:]
    lat = f['NS/Latitude'][:]
    surfacePrecipitation = f['NS/surfPrecipTotRate'][:]
 
time = '2014-07-03 07:31:30-09:04:03'

extents = [106, 111, 26,31]
var=[lon,lat,surfacePrecipitation]
var_new= region_mask(var, extents)
lon_c = var_new[0]
lat_c = var_new[1]
surfacePrecipitation = var_new[2]

#设置缺测
surfacePrecipitation.replace(-9999.9,np.nan)

levels_Rr = [0.1, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0, 15.0, 20.0, 25]
cmap_Rr, norm_Rr = make_Rr_cmap(levels_Rr)

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.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(xlim=(extents[0], extents[1]),ylim=(extents[-2], extents[-1]))
'''
url = 'http://map1c.vis.earthdata.nasa.gov/wmts-geo/wmts.cgi'
layer = 'BlueMarble_ShadedRelief'
ax.add_wmts(url, layer)  
#使用上面三行的底图需先下载owslib
#conda install -c conda-forge owslib
'''

ax.grid()
ax.tick_params(labelsize='large')
im = ax.contourf(
        lon_c, lat_c, surfacePrecipitation, levels_Rr,
        cmap=cmap_Rr, norm=norm_Rr, extend='both', transform=proj
    )
cb = fig.colorbar(im, ax=ax, ticks=levels_Rr)
cb.set_label('Rain Rate (mm/hr)', fontsize='large')
cb.ax.tick_params(labelsize='large')
ax.set_xticks(np.arange(extents[0], extents[1]+1, 1), crs=proj)
ax.set_yticks(np.arange(extents[-2], extents[-1]+1, 1), crs=proj)
ax.set_title(f'Rain Rate on {time}', fontsize='x-large')
plt.show()

 出图：

本代码更新版见Python | GPM CMB 单个文件SurfPrecipTotRate 绘制_HHHH1014的博客-优快云博客