【气象编程】利用ERA5数据计算涡度平流并绘图

利用ERA5数据计算涡度平流并绘图

1. 官网示例（基于NCSS）

metpy给出的涡度平流计算绘图代码链接: 500hpa_vorticity_advection
整体流程是读取数据，计算涡度，计算涡度的平流，然后绘图。

存在问题

示例中使用的数据结构和要使用的ERA5略有不同，此外，由于版本问题，示例中使用的以下计算语句由于position argument设置的原因，导致无法使用。

avor = mpcalc.vorticity(uwnd_500, vwnd_500, dx, dy, dim_order='yx') + f

avor = ndimage.gaussian_filter(avor, sigma=3, order=0) * units('1/s')

vort_adv = mpcalc.advection(avor, [uwnd_500, vwnd_500], (dx, dy), dim_order='yx') * 1e9

2. 解决方法

• 导入库

import datetime as dt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
import numpy as np
import scipy.ndimage as ndimage
import xarray as xr
from metpy.units import units

• 读取指定日期的文件

target_time = dt.datetime(2024, 7, 31)  # 举个例子读24年7月31的数据 
tiem_str = target_time.strftime('%Y%m')

# 以下路径设置成自己数据的路径
ERA5_Dir = '../Data\\ERA5_Reanalysis\\PressureLevels'
f_z = xr.open_dataset(f'{ERA5_DIR}\\Hgt\\Hgt-{tiem_str}.nc')  # 高度
f_u = xr.open_dataset(f'{ERA5_DIR}\\Uwnd\\Uwnd-{tiem_str}.nc')  # 风速u分量
f_v = xr.open_dataset(f'{ERA5_DIR}\\Vwnd\\Vwnd-{tiem_str}.nc')  # 风速v分量

• 提取500hPa的变量，并加入单位

z = f_z['z'].sel(time=target_time, level=500)
u = f_u['u'].sel(time=target_time, level=500)
v = f_v['v'].sel(time=target_time, level=500) 
lon = f_z['longitude'].values
lat = f_z['latitude'].values

hght_500 = z.values / 98
hght_500 = ndimage.gaussian_filter(hght_500, sigma=3, order=0) * units.meter

uwnd_500 = units('m/s') * u.values
vwnd_500 = units('m/s') * v.values

• 计算涡度和涡度平流

dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)

f = mpcalc.coriolis_parameter(np.deg2rad(lat)).to(units('1/sec'))

avor = mpcalc.vorticity(u,v) + f

avor = ndimage.gaussian_filter(avor, sigma=3, order=0) * units('1/s')

vort_adv = mpcalc.advection(avor, uwnd_500, vwnd_500, dx=dx, dy=dy) * 1e9

使用xarray.dataarray变量类型作为mpcalc.vorticity的输入时，不需要另外输入dx, dy，通过这种方式，便捷地计算涡度。随后，另一版本的mpcalc.advection不再支持原来的dim_order的参数输入，其他参数的输入方式也略有不同，按照官方文档，进行相应的修改（如上述代码），即可完成平流计算。

• 画图
  根据样例的画图程序进行画图即可

# Set up Coordinate System for Plot and Transforms
plotcrs = ccrs.PlateCarree()
fig = plt.figure(1, figsize=(14., 12))
gs = gridspec.GridSpec(2, 1, height_ratios=[1, .02], bottom=.07, top=.99,
                       hspace=0.01, wspace=0.01)
ax = plt.subplot(gs[0], projection=plotcrs)

# Plot Titles
plt.title(r'500-hPa Heights (m), AVOR$*10^5$ ($s^{-1}$), AVOR Adv$*10^8$ ($s^{-2}$)',
          loc='left')
plt.title(f'VALID: {tiem_str}', loc='right')

# Plot Background
ax.set_extent([leftlon, rightlon, lowerlat, upperlat], ccrs.PlateCarree())
ax.coastlines('50m', edgecolor='black', linewidth=0.75)
ax.add_feature(cfeature.STATES, linewidth=.5)

# Plot Height Contours
clev500 = np.arange(5100, 6061, 60)
cs = ax.contour(lon, lat, hght_500.m, clev500, colors='black', linewidths=1.0,
                linestyles='solid', transform=ccrs.PlateCarree())
plt.clabel(cs, fontsize=10, inline=1, inline_spacing=10, fmt='%i',
           rightside_up=True, use_clabeltext=True)

# Plot Absolute Vorticity Contours
clevvort500 = np.arange(-9, 50, 5)
cs2 = ax.contour(lon, lat, avor*10**5, clevvort500, colors='grey',
                 linewidths=1.25, linestyles='dashed', transform=ccrs.PlateCarree())
plt.clabel(cs2, fontsize=10, inline=1, inline_spacing=10, fmt='%i',
           rightside_up=True, use_clabeltext=True)

# Plot Colorfill of Vorticity Advection
clev_avoradv = np.arange(-30, 31, 5)
cf = ax.contourf(lon, lat, vort_adv.m, clev_avoradv[clev_avoradv != 0], extend='both',
                 cmap='bwr', transform=ccrs.PlateCarree())
cax = plt.subplot(gs[1])
cb = plt.colorbar(cf, cax=cax, orientation='horizontal', extendrect='True', ticks=clev_avoradv)
cb.set_label(r'$1/s^2$', size='large')

# Plot Wind Barbs
# Transform Vectors and plot wind barbs.
ax.barbs(lon, lat, uwnd_500.m, vwnd_500.m, length=6, regrid_shape=20,
         pivot='middle', transform=ccrs.PlateCarree())

作图效果大致如图所示：