ERA5 | 台风 | 基于ERA5数据的台风摩羯的气象动图制作

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)