代码实战 | 如何绘制台风方位角平均要素

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib.ticker import LogLocator, NullFormatter
from besttracks import parse_TCs
from xvortices import load_cylind, project_to_cylind

# ========================
# 1. 加载热带气旋最佳路径数据
# ========================
haima = parse_TCs('*.txt', agency='CMA',
                  tc_cond=(lambda tc: tc.ID == '200421'))[0]

# ========================
# 2. 读取格点气象数据
# ========================
dset = xr.open_dataset('/home/mw/project/Haima2004.nc')  # 请改为实际路径

# ========================
# 3. 插值到移动的柱坐标系
# ========================
olon = haima.get_as_xarray('LON')[1:]   # 时间对齐
olat = haima.get_as_xarray('LAT')[1:]

[u, v, w, h], lons, lats, etas = load_cylind(
    dset, olon, olat,
    azimNum=72, radiNum=31, radMax=6,
    lonname='lon', latname='lat'
)

# ========================
# 4. 计算风暴相对风场并投影到径向/切向
# ========================
haima.translate_velocity()
uo = haima.get_as_xarray('uo')
vo = haima.get_as_xarray('vo')

u_r = u - uo          # 相对纬向风
v_r = v - vo          # 相对经向风

uaz_r, vra_r = project_to_cylind(u_r, v_r, etas)  # 切向风、径向风

# ========================
# 5. 方位角平均
# ========================
uaz_rm = uaz_r.mean('azim')
vra_rm = vra_r.mean('azim')
wm = w.mean('azim')
hm = h.mean('azim')

# ========================
# 6. 按气压升序排序（关键修正！）
# ========================
uaz_rm = uaz_rm.sortby('lev')
vra_rm = vra_rm.sortby('lev')
wm = wm.sortby('lev')
hm = hm.sortby('lev')

# ========================
# 7. 绘制方位角平均结构（Matplotlib）
# ========================
tlev = 2   # 时间索引
fig, axes = plt.subplots(2, 2, figsize=(9.5, 7), sharex=True, sharey=True)

# ---- 图1：平均切向风 ----
ax = axes[0, 0]
cf = ax.contourf(uaz_rm.radi, uaz_rm.lev, uaz_rm.isel(time=tlev),
                 levels=np.linspace(-14, 14, 29), cmap='seismic', extend='both')
ax.set_title('mean azimuthal wind')
fig.colorbar(cf, ax=ax)

# ---- 图2：平均径向风 ----
ax = axes[0, 1]
cf = ax.contourf(vra_rm.radi, vra_rm.lev, vra_rm.isel(time=tlev),
                 levels=np.linspace(-4, 4, 17), cmap='seismic', extend='both')
ax.set_title('mean radial wind (relative)')
fig.colorbar(cf, ax=ax)

# ---- 图3：平均垂直速度 ----
ax = axes[1, 0]
cf = ax.contourf(wm.radi, wm.lev, wm.isel(time=tlev),
                 levels=np.linspace(-1, 0, 21), cmap='Reds_r', extend='both')
ax.set_title('mean vertical wind')
fig.colorbar(cf, ax=ax)

# ---- 图4：平均位势高度异常（减去外边界值） ----
ax = axes[1, 1]
anomaly = hm.isel(time=tlev) - hm.isel(time=tlev).isel(radi=-1)
cf = ax.contourf(anomaly.radi, anomaly.lev, anomaly,
                 levels=np.linspace(-50, 50, 21), cmap='seismic', extend='both')
ax.set_title('mean geopotential height anomaly')
fig.colorbar(cf, ax=ax)

# ========================
# 坐标轴设置（对数压力轴，上低下高）
# ========================
for ax in axes.flat:
    ax.set_xlabel('radius (degree)')
    ax.set_ylabel('pressure (hPa)')
    ax.set_yscale('log')
    ax.yaxis.set_major_locator(LogLocator(base=10))
    ax.yaxis.set_minor_locator(LogLocator(base=10, subs='all'))
    ax.yaxis.set_minor_formatter(NullFormatter())
    ax.invert_yaxis()                    # 让 1000 在下，100 在上
    ax.set_ylim(1000, 100) 
    ax.grid(True, linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()