基于ERA5数据的台风摩羯的气象动图制作
前言
虽然之前介绍了geogif库制作动图，但是无法解决设置地图和图片精度的问题
于是现在出一期基于geocat和matplotlib的era5动图制作
项目链接：
https://www.heywhale.com/mw/project/66fec2d707a4464c3fc75d4b
导入库与读取数据
import cartopy.crs as ccrs
import matplotlib.animation as animation
import numpy as np
import xarray as xr
from matplotlib import pyplot as plt
import geocat.viz as gv
import cartopy.feature as cfeature
from meteva.base.tool.plot_tools import add_china_map_2basemap
ds = xr.open_dataset('/home/mw/input/era59051/era5_single_20240904_0908.nc')
print(ds)<xarray.Dataset> Size: 2GB
Dimensions: (valid_time: 96, latitude: 721, longitude: 1440)
Coordinates:
number int64 8B ...
* valid_time (valid_time) datetime64[ns] 768B 2024-09-04 ... 2024-09-07T23...
* latitude (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
* longitude (longitude) float64 12kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
expver (valid_time) <U4 2kB ...
Data variables:
u10 (valid_time, latitude, longitude) float32 399MB ...
v10 (valid_time, latitude, longitude) float32 399MB ...
d2m (valid_time, latitude, longitude) float32 399MB ...
msl (valid_time, latitude, longitude) float32 399MB ...
mtpr (valid_time, latitude, longitude) float32 399MB ...
tcrw (valid_time, latitude, longitude) float32 399MB ...
Attributes:
GRIB_centre: ecmf
GRIB_centreDescription: European Centre for Medium-Range Weather Forecasts
GRIB_subCentre: 0
Conventions: CF-1.7
institution: European Centre for Medium-Range Weather Forecasts
history: 2024-09-24T07:56 GRIB to CDM+CF via cfgrib-0.9.1...数据处理
slp = ds.msl/98
## 修改时间维度名称
slp = slp.rename({'valid_time':'time'}).isel(time=slice(0, 10))print(slp)<xarray.DataArray 'msl' (time: 10, latitude: 721, longitude: 1440)> Size: 42MB
array([[[1022.78827, 1022.78827, 1022.78827, ..., 1022.78827,
1022.78827, 1022.78827],
[1023.70917, 1023.70917, 1023.71173, ..., 1023.7066 ,
1023.7066 , 1023.7066 ],
[1024.2347 , 1024.2373 , 1024.2397 , ..., 1024.2296 ,
1024.2322 , 1024.2347 ],
...,
[1021.051 , 1021.051 , 1021.04846, ..., 1021.04846,
1021.04846, 1021.04846],
[1020.5969 , 1020.59436, 1020.59436, ..., 1020.5969 ,
1020.5969 , 1020.5969 ],
[1020.0434 , 1020.0434 , 1020.0434 , ..., 1020.0434 ,
1020.0434 , 1020.0434 ]],
[[1022.625 , 1022.625 , 1022.625 , ..., 1022.625 ,
1022.625 , 1022.625 ],
[1023.6071 , 1023.6071 , 1023.6071 , ..., 1023.6046 ,
1023.6046 , 1023.6046 ],
[1024.1403 , 1024.1428 , 1024.1454 , ..., 1024.1377 ,
1024.1377 , 1024.1403 ],
...
[1020.95154, 1020.95917, 1020.96173, ..., 1020.93115,
1020.9362 , 1020.94385],
[1020.69135, 1020.69385, 1020.69385, ..., 1020.68365,
1020.6862 , 1020.6888 ],
[1020.2857 , 1020.2857 , 1020.2857 , ..., 1020.2857 ,
1020.2857 , 1020.2857 ]],
[[1020.18494, 1020.18494, 1020.18494, ..., 1020.18494,
1020.18494, 1020.18494],
[1021.3049 , 1021.3049 , 1021.3074 , ..., 1021.3023 ,
1021.3023 , 1021.3023 ],
[1021.9171 , 1021.9171 , 1021.9196 , ..., 1021.912 ,
1021.912 , 1021.91455],
...,
[1021.1594 , 1021.1671 , 1021.1722 , ..., 1021.13904,
1021.1441 , 1021.1518 ],
[1020.94006, 1020.9426 , 1020.9451 , ..., 1020.9324 ,
1020.93494, 1020.9375 ],
[1020.7768 , 1020.7768 , 1020.7768 , ..., 1020.7768 ,
1020.7768 , 1020.7768 ]]], dtype=float32)
Coordinates:
number int64 8B 0
* time (time) datetime64[ns] 80B 2024-09-04 ... 2024-09-04T09:00:00
* latitude (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
* longitude (longitude) float64 12kB 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
expver (time) <U4 160B ...制作底图
fig = plt.figure(figsize=(10, 8))
# Generate axes using Cartopy and draw coastlines
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([100, 140, 10, 40], ccrs.PlateCarree())
add_china_map_2basemap(ax, name="nation",edgecolor='k', lw=0.5, encoding='gbk', grid0=None) # "国界"
# Use geocat-viz convenience function to set axes limits & tick values
gv.set_axes_limits_and_ticks(ax,
xlim=(100, 140),
ylim=(10,40),
xticks=np.linspace(100, 140, 5),
yticks=np.linspace(10, 40, 4))
# Use geocat-viz convenience function to add minor and major tick lines
gv.add_major_minor_ticks(ax, labelsize=10)
# Use geocat-viz convenience function to make latitude, longitude tick labels
gv.add_lat_lon_ticklabels(ax)
# Create initial plot
cplot= slp[0, :, :].plot.contourf(ax=ax,
transform=ccrs.PlateCarree(),
vmin=980,
vmax=1050,
levels=15,
cmap="inferno",
add_colorbar=False)
# Create a colorbar
cbar = fig.colorbar(cplot,
extendrect=True,
orientation="vertical",
# ticks=np.arange(220, 320, 5),
label="",
shrink=0.90)
# Remove minor ticks from colorbar that don't work well with other formatting
cbar.ax.minorticks_off()
动图制作
# Animate function for matplotlib FuncAnimation
def animate(i):
slp[i, :, :].plot.contourf(ax=ax,
transform=ccrs.PlateCarree(),
vmin=980,
vmax=1050,
levels=15,
cmap="inferno",
add_colorbar=False)
gv.set_titles_and_labels(
ax,
maintitle="typhoon YAGI SLP- TIME= " +
str(slp.coords['time'].values[i])[:13],
xlabel="",
ylabel="")
anim = animation.FuncAnimation(fig, animate, frames=10, interval=200)
anim.save('animate_1.gif', writer='pillow', fps=5)
在project文件夹下即可打开gif查看效果。博主还没搞懂怎么做到粘贴动图
