如何在一张图上同时绘制云图和降水
前言
需求：大家看到诸多文献使用卫星云图作为天气形势系统介绍时想必也想自己也为文章中加一张，那么卫星云图如何叠加降水图呢
面向群体：需要使用卫星云图进行天气学分析或天气系统阐释的小伙伴，当然你喜欢卫星云图做壁纸也可以画着玩
应用场景：汇报or写作
1. 导入库并可视化
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import metpy
import glob
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
#读取葵花8号卫星多颗netcdf数据
files = glob.glob('/home/mw/input/himawari7828/0800_20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_2000-HIMAWARI8-AHI/*.nc')
ds = xr.open_mfdataset(files)
#裁剪感兴趣的区域
ds1 = ds.isel(x=slice(1000, 4000), y=slice(1000, 3500), time=0)
#获取投影参考信息
central_longitude = ds1['geostationary'].longitude_of_projection_origin
satellite_height = ds1['geostationary'].satellite_height
#将km单位转换为m
mapx = ds1['x'].to_numpy() * 1000.
mapy = ds1['y'].to_numpy() * 1000.
##提取3个波段为RGB
rgb = np.stack((
ds1['channel_0001_scaled_radiance'].to_numpy(),
ds1['channel_0002_scaled_radiance'].to_numpy(),
ds1['channel_0003_scaled_radiance'].to_numpy(),
), axis=-1)
#扩充RGB值范围使对比度增强
rgb_stretched = np.clip((rgb/0.7)**0.7, 0, 1)
#定义投影坐标系
data_proj = ccrs.Geostationary(
central_longitude=central_longitude,
satellite_height=satellite_height,
)
map_proj = ccrs.PlateCarree()
#绘图
fig, ax = plt.subplots(
figsize=(6, 6), facecolor="w", dpi=200, layout='constrained',
subplot_kw=dict(projection=map_proj)
)
#将数据绘制到底图上
pcm = ax.pcolorfast(mapx, mapy, rgb_stretched, transform=data_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
gl.yformatter = LATITUDE_FORMATTER
#添加海岸线
ax.coastlines(color="w")
#设置地图范围
ax.set_extent()
#设置中文标题
ax.set_title('葵花8号卫星真彩色影像')
为什么修改RGB
原始的RGB值范围是0-1。直接使用会使图像对比度不高,颜色看起来比较浅淡。通过gamma校正等方法将这种线性关系转换为非线性,使较暗的区域变亮,较亮的区域保持不变。这样可以增加整个图像的对比度,使颜色更加饱和丰富
为什么修改单位km为m
图投影坐标系一般使用的是米为单位。直接拿千米单位的影像坐标去绘制地图,会造成非常严重的坐标错位。
因此需要提前将影像的坐标单位换算为与地图投影匹配的米单位,然后再传入投影变换,进行坐标转换到地图上。
为什么使用pcolorfast
对于绘制地图影像,pcolorfast能够提供更快速和直接的解决方案。它适合直接可视化大规模的不规则网格数据,比如常见的卫星影像等。是地图绘制过程中的一种非常有效和高效的方法
2.绘制era5小时降水
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr
pr = xr.open_dataset('/home/mw/project/2016081408_prep.grib',engine='pynio')
prhour=pr.TP_GDS0_SFC_acc1h * 1000
# 创建地图投影
projection = ccrs.PlateCarree()
# 创建地图和子图对象
fig = plt.figure(figsize=(10, 6),dpi=200)
ax = fig.add_subplot(1, 1, 1, projection=projection)
lat =pr.g0_lat_0
lon =pr.g0_lon_1
# 绘制地理数据
#cf = ax.pcolormesh(lon,lat,prhour,vmin=0,vmax=30,cmap=cmaps.WhiteBlue,transform=ccrs.PlateCarree())
cf = ax.contourf(lon,lat,prhour,levels=np.arange(0,30,2),cmap=plt.cm.Spectral_r,transform=ccrs.PlateCarree())
ax.set_extent()
cbar = plt.colorbar(cf)
ax.coastlines(resolution='50m', color='black', linewidth=1)
ax.add_feature(ccrs.cartopy.feature.STATES, linewidth=1)
ax.set_title('ERA5 hour pr')
# 显示图形
plt.show()
3. era5降水叠加葵花卫星云图
import xarray as xr
pr = xr.open_dataset('/home/mw/project/2016081408_prep.grib',engine='pynio')
prhour=pr.TP_GDS0_SFC_acc1h * 1000
fig, ax = plt.subplots(
figsize=(6, 6), facecolor="w", dpi=200, layout='constrained',
subplot_kw=dict(projection=map_proj)
)
lat =pr.g0_lat_0
lon =pr.g0_lon_1
# 绘制降水数据
#cf = ax.pcolormesh(lon,lat,prhour,vmin=0,vmax=30,cmap=cmaps.WhiteBlue,transform=ccrs.PlateCarree())
cf = ax.contourf(lon,lat,prhour,levels=np.arange(0,30,2),cmap=plt.cm.Spectral_r,transform=ccrs.PlateCarree(),alpha=0.7)
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轴设为纬度格式
ax.coastlines(color="w")
ax.set_extent()
#绘制卫星数据
pcm = ax.pcolorfast(mapx, mapy, rgb_stretched, transform=data_proj,alpha=0.75)
#设置中文标题
ax.set_title('葵花8号卫星真彩色影像叠加降水数据')
叠加图关键在于填色图的alpha参数，可以理解为透明度，觉得效果不佳可以自行调整。
文章来源于微信公众号：气ython风雨
页:
[1]