增强云图
FY-4A已经在上半年结束在我国上空的工作,FY-4B经过漂移变轨进入原FY-4A的泊位,现在是气象日常天气诊断分析主力静止卫星,但是Micaps系统及各种气象网站公开的卫星云图配色不是非常容易识别,这里学习了一下中气爱使用的一种配色方案,使用这种增强云图具有更好的识别性。
最近台风摩羯生成并即将登陆我国,台风系统特别适合使用用来表现增强云图,所以就近选取一个台风已经开眼时刻来绘图。
增强云图的配色方案参考了中气爱分享的图片配色,但色号为我手动生成。FY-4B卫星数据的读取及处理参考了墨大宝开发的处理程序段(https://github.com/Mo-Dabao)。
FY-4B数据下载
卫星数据开放性较好,有国家卫星气象中心的官网可供下载,这里下载FY-4B的4KM全圆盘数据。按照要求注册账号,提交订单后,数据中心会回调,回调结束后就可以按照方法下载数据。
增强云图色系的生成
中气爱的原图如下:
分析这个图像,可以看出这是一个多段拼接的色条,其中零下20℃至零上部分为灰色系;零下20℃至零下80℃为一个多色混合色条;更低温度部分为粉色系。
首先生成灰色色条,matplotlib有预置的灰色色条Greys,假定为1℃为一节,则灰色部分我们需要生成30个阶梯色:
然后借助LinearSegmentedColormap生成混合色条
最后通过粉色色条生成最后一段颜色:
Pink_part=plt.get_cmap('RdPu_r',256)(np.linspace(0,1,256))[10::10,:][0:10]
将生成好的多个颜色条堆叠在一起,生成整体颜色条:
解译卫星数据
墨大宝老师开发的这个模块特别好用,直接返回的是xarray的DataArray格式。
file='FY4B-_AGRI--_N_DISK_1050E_L1-_FDI-_MULT_NOM_20240905013000_20240905014459_4000M_V0001.HDF'
f=h5py.File(file)
extent=[0,70,50,160,0.05]
f=AGRI_L1(file,extent)
infrared_10_8=f.extract('Channel14',calibration='brightness_temperature')-273.15
定义地图绘图类
这个绘图类是我做的一个小项目里的副产品,可以给子图添加地图,节省调试时间。
import geopandas as gpd
import cartopy.feature as cfeature
class Map():
def __init__(self,path):
self.hubei=gpd.read_file(path+r'\hubei_quansheng\hubei_quansheng.shp')
self.china=gpd.read_file(path+r'\china\china.shp')
self.nineline=gpd.read_file(path+r'\nine_line\nine_line.shp')
self.hubeionly=gpd.read_file(path+r'\hubei\hubei.shp')
self.centralchina=gpd.read_file(path+r'\province\province.shp',encoding='utf-8')
def Lambert_Map(self,**kwargs):
ax=plt.gca()
ax.add_geometries(self.hubeionly.geometry,
lw=0.3,crs=ccrs.PlateCarree(),edgecolor='tab:red',facecolor='none',zorder=5,**kwargs)
ax.add_geometries(self.china.geometry,
lw=0.35,crs=ccrs.PlateCarree(),edgecolor='k',facecolor='none',zorder=5,**kwargs)
ax.add_geometries(self.nineline.geometry,
lw=0.5,crs=ccrs.PlateCarree(),edgecolor='k',facecolor='none',zorder=5,**kwargs)
ax.add_feature(cfeature.LAND.with_scale('50m'),color='#FFFFF9')
ax.add_feature(cfeature.OCEAN.with_scale('50m'),color='#EDFBFB')
gl=ax.gridlines(draw_labels=True,linewidth=0.3,linestyle='--',color='grey',alpha=0.75)
gl.xlocator=mticker.FixedLocator(np.arange(60,141,10))
gl.ylocator=mticker.FixedLocator(np.arange(10,51,10))
gl.top_labels=False
gl.right_labels=False
gl.bottom_labels=True
gl.rotate_labels=None
gl.xformatter=cmg.LONGITUDE_FORMATTER
gl.yformatter=cmg.LATITUDE_FORMATTER
gl.x_inline=False
gl.y_inline=False
gl.xlabel_style={'size':6,'fontfamily':'Times New Roman'}
gl.ylabel_style={'size':6,'fontfamily':'Times New Roman'}
ax.set_extent([80,136,13.5,56],crs=ccrs.PlateCarree())
nanhai=ax.inset_axes([0.7952,0,0.25,0.25],projection=ccrs.PlateCarree())
nanhai.coastlines(lw=0.3)
nanhai.add_geometries(self.china.geometry,
lw=0.3,crs=ccrs.PlateCarree(),edgecolor='k',facecolor='none',zorder=5)
nanhai.add_geometries(self.nineline.geometry,
lw=0.55,crs=ccrs.PlateCarree(),edgecolor='grey',facecolor='none',zorder=5)
nanhai.set_extent([105,125,0,25])
nanhai.add_feature(cfeature.LAND.with_scale('110m'),color='#FFFFF9')
Map_class=Map(r'C:\Users\Administrator\Desktop\ERA5lightning')
绘制图像
#定义温度等级
levs=np.concatenate((np.arange(-90,-20,1),np.arange(-20,41,2)))
#颜色与数值生成映射关系
norm=mcolors.BoundaryNorm(levs,brightness_temperature_cmap.N)
from mpl_toolkits.axes_grid1 import make_axes_locatable
import cmaps
plt.rcParams['font.sans-serif']=['FangSong']
plt.rcParams['axes.unicode_minus']=False
fig=plt.figure(figsize=(5,5),dpi=800)
ax=fig.add_axes([0,0,1,1],projection=ccrs.LambertConformal(central_longitude=110,central_latitude=30))
Map_class.Lambert_Map()
#主图画图
ac=ax.pcolormesh(infrared_10_8.lon,
infrared_10_8.lat,
infrared_10_8.values,
cmap=brightness_temperature_cmap,
norm=norm,
transform=ccrs.PlateCarree())
#南海小地图画图
ax.child_axes[0].pcolormesh(infrared_10_8.lon,
infrared_10_8.lat,
infrared_10_8.values,
cmap=brightness_temperature_cmap,
norm=norm,
transform=ccrs.PlateCarree())
divider=make_axes_locatable(ax)
cax=divider.new_horizontal(size="3%", pad=0.1,axes_class=plt.Axes)
fig.add_axes(cax)
cb=fig.colorbar(ac,cax=cax)
cb.ax.tick_params(labelsize=7)
cb.ax.set_yticks(np.concatenate((np.arange(-90,-21,10),np.arange(-20,41,20))))
fig.suptitle(x=0.08,y=0.93,t='FY-4B\n 2024-09-05 01:30 BJ',fontsize=6)
ax.set_title('增强红外云图',fontsize=15)
这里偷懒使用了现成的程序,所以台风还在角落里,不过已经能发现台风摩羯中心开眼,环流形状完整。预计将会给我国南部沿海省份带来一番暴风骤雨。