GPM 卫星三级产品的简单可视化
本帖最后由 气ython风雨 于 2024-4-30 12:36 编辑之前GPM数据写了下载教程，现在简单试试可视化，毕竟是nc格式数据（下载可选），用起来相对简单
1.1 导入库与查看数据
# 导入库
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
# 读取
ds = xr.open_dataset('/home/mw/input/GPM5633/gpmPython/3B-DAY.MS.MRG.3IMERG.20210519-S000000-E235959.V06.nc4')
print(ds)
<xarray.Dataset>
Dimensions: (time: 1, lon: 3600, lat: 1800, nv: 2)
Coordinates:
* lon (lon) float32 -179.9 -179.8 ... 179.9 179.9
* lat (lat) float32 -89.95 -89.85 ... 89.85 89.95
* time (time) object 2021-05-19 00:00:00
Dimensions without coordinates: nv
Data variables:
precipitationCal (time, lon, lat) float32 ...
precipitationCal_cnt (time, lon, lat) int8 ...
precipitationCal_cnt_cond(time, lon, lat) int8 ...
HQprecipitation (time, lon, lat) float32 ...
HQprecipitation_cnt (time, lon, lat) int8 ...
HQprecipitation_cnt_cond (time, lon, lat) int8 ...
randomError (time, lon, lat) float32 ...
randomError_cnt (time, lon, lat) int8 ...
time_bnds (time, nv) object ...
Attributes:
BeginDate: 2021-05-19
BeginTime: 00:00:00.000Z
EndDate: 2021-05-19
EndTime: 23:59:59.999Z
FileHeader: StartGranuleDateTime=2021-05-19T00:00:00.000Z;\nStopGran...
InputPointer: 3B-HHR.MS.MRG.3IMERG.20210519-S000000-E002959.0000.V06B....
title: GPM IMERG Final Precipitation L3 1 day 0.1 degree x 0.1 ...
DOI: 10.5067/GPM/IMERGDF/DAY/06
ProductionTime:2021-08-31T19:44:25.000Z
可见该数据为2021年5月19日的日降水数据，经纬度分辨率为0.1°
ProductionTime指的是产品生成时间，因为这是三级数据，在三个月后生成是正常的。具体可见官网解释产品分级。
1.2 IMERG降水可视化
加入修正cartopy无法正常显示经纬度的函数
'''
参考Andrew Dawson 提供的解决方法: <a href="https://gist.github.com/ajdawson/dd536f786741e987ae4e" target="_blank">https://gist.github.com/ajdawson/dd536f786741e987ae4e</a>
'''
# 给出一个假定为矩形的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)
# 在兰伯特投影的底部X轴上绘制刻度线
def lambert_xticks(ax, ticks):
"""Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
te = lambda xy: xy
lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b, b, n))).T
xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te)
ax.xaxis.tick_bottom()
ax.set_xticks(xticks)
ax.set_xticklabels()
# 在兰伯特投影的左侧y轴上绘制刻度线
def lambert_yticks(ax, ticks):
"""Draw ricks on the left y-axis of a Lamber Conformal projection."""
te = lambda xy: xy
lc = lambda t, n, b: np.vstack((np.linspace(b, b, 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()
# 获取兰伯特投影中底部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 =
else:
tick = tick_extractor(locs.xy)
_ticks.append(tick)
# 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
ps: 需要使用np.transpose()进行倒转维度，不然画不出来
#precipitationcal是指质控的降水数据
prep=ds.precipitationCal
lon=prep.lon
lat=prep.lat
prep=prep.transpose()
lon, lat = np.meshgrid(lon, lat)
# 设置地图投影和范围
proj = ccrs.LambertConformal(central_longitude=105, central_latitude=90, standard_parallels=(25, 47))
extent =
# 创建画布和子图
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()
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()
1.3 GMI(hdf5)
# 读取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)
prnsdata = np.ndarray(shape=prns.shape,dtype=float)
prns.read_direct(prnsdata)
prnsdata = np.nan
start = 0
end = 43
mysub = prnsdata
mylon = lon
mylat = lat
prnsdata.min()
# 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()
cbar = plt.colorbar(pp,cax=cbar_ax,orientation='horizontal')
1.4 GMI可视化nc版
In :
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 =
surfacePrecipitation = ds1.S1_surfacePrecipitation
lon=surfacePrecipitation.S1_Longitude
lat=surfacePrecipitation.S1_Latitude
#prep=prep.transpose()
levels_Rr =
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_
lon = lon
lat = lat
surfacePrecipitation = surfacePrecipitation
# 画出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...
完整代码与文件在此
文章来源于微信公众号：气ython风雨
666
页:
[1]