数值模式 | WRF 台风方位角平均半径-气压剖面绘制教程

用户11172986

发布于 2026-05-07 19:15:41

1150

文章被收录于专栏：气python风雨气python风雨

WRF 台风方位角平均半径-气压剖面绘制教程

前言

上期做了格点数据计算台风方位角平均的教程，后面留了一个坑：wrf的怎么搞

wrf当然没有现成的库直接用，直接手搓看看吧

本教程做了一些优化，核心思路是：

跳过垂直插值，直接在WRF原始模式面上计算，保留全部垂直分辨率；
Box裁剪：只取台风中心附近的子网格做水平插值，避免在全网格上浪费计算；
规则网格插值：利用wrf.ll_to_xy将目标经纬度映射到WRF的网格坐标索引，再用scipy.interpolate.RegularGridInterpolator在规则网格上插值，速度比griddata快一个数量级以上；
气象方位角公式：metpy的azimuth_range_to_lat_lon使用0°=北、顺时针的气象方位角，切向风/径向风分解公式需与之匹配。

1. 环境准备

需要安装以下Python包：

wrf-python
metpy
netCDF4
numpy、scipy、matplotlib

如果你使用conda环境，可以通过以下命令安装：

conda install -c conda-forge wrf-python metpy netCDF4

import sys, warnings
from pathlib import Path
from glob import glob
import numpy as np
from netCDF4 import Dataset
from scipy.interpolate import RegularGridInterpolator
import matplotlib.pyplot as plt
import wrf
import metpy.calc as mpcalc
from metpy.units import units

plt.rcParams['font.size'] = '14'

2. 读取WRF数据并定位台风中心

台风中心通过最低海平面气压（SLP）定位。wrf-python提供了wrf.getvar(nc, "slp")直接计算SLP，然后取二维平面上的最小值位置即可得到中心经纬度。

同时读取纬向风ua、经向风va以及各层气压pressure，这里我们直接在模式面（bottom_top层）上操作，不做垂直插值。

# 设置参数
MAX_R = 300.0       # 最大半径(km)
NR = 101            # 径向格点数
NAZ = 73            # 方位角格点数
INPATH = "/home/mw/input/typhoon9537"

# 定位文件
p = Path(INPATH)
if p.is_file():
    f = str(p)
elif p.is_dir():
    f = (sorted(glob(str(p / "wrfout_d01_*"))) or sorted(glob(str(p / "wrfout_*"))))[0]
else:
    f = sorted(glob(INPATH))[12]

print(f"读取文件: {f}")
nc = Dataset(f)

# 利用最低SLP定位台风中心
slp = wrf.getvar(nc, "slp").values
min_idx = np.unravel_index(np.argmin(slp), slp.shape)
y0, x0 = min_idx
lat_c = float(wrf.getvar(nc, "lat").values[min_idx])
lon_c = float(wrf.getvar(nc, "lon").values[min_idx])
print(f"台风中心: {lat_c:.2f}°N, {lon_c:.2f}°E  最低SLP: {slp[min_idx]:.1f} hPa")

读取文件: /home/mw/input/typhoon9537/wrfout_d01_2019-08-08_18_00_00
台风中心: 25.52°N, 124.48°E  最低SLP: 923.1 hPa

3. 构建极坐标目标网格

利用metpy的azimuth_range_to_lat_lon生成围绕台风中心的极坐标目标网格。该函数接受气象方位角（0°=北，顺时针）和距离，返回对应的经纬度坐标。

随后通过wrf.ll_to_xy将这些目标经纬度映射回WRF的网格坐标索引（x, y）。这一步非常关键：它把曲率的经纬度问题转化为了规则网格坐标问题，后续插值就可以在规则的（y, x）索引上进行。

# 极坐标网格: 0~360°, 0~300km
az = np.linspace(0, 360, NAZ)
r = np.linspace(0, MAX_R, NR) * 1000# 转换为米

# 生成目标点的经纬度
lon_a, lat_a = mpcalc.azimuth_range_to_lat_lon(
    az * units.degree, r * units.meter, lon_c, lat_c
)

# 映射到WRF网格坐标 (x, y)
x_t, y_t = wrf.ll_to_xy(
    nc, np.asarray(lat_a).ravel(), np.asarray(lon_a).ravel()
)
x_t = x_t.values.astype(float)
y_t = y_t.values.astype(float)

print(f"目标网格范围: x={x_t.min():.1f}~{x_t.max():.1f}, y={y_t.min():.1f}~{y_t.max():.1f}")

目标网格范围: x=173.0~371.0, y=65.0~263.0

4. Box裁剪：只取中心附近的子网格

我们不需要在整个WRF大网格（本例为530x755）上做插值。只需取目标点覆盖的最小外接矩形，再向外扩展几个格点（pad=3）作为缓冲，即可将计算量降低60%以上。

读取子网格的ua、va、pressure后，Y轴坐标用台风中心所在列的各层气压值表示，这样剖面图的纵坐标就是实际的气压高度。

ny, nx = len(nc.dimensions["south_north"]), len(nc.dimensions["west_east"])

# Box裁剪
pad = 3
xmin = max(0, int(x_t.min()) - pad)
xmax = min(nx, int(x_t.max()) + pad + 1)
ymin = max(0, int(y_t.min()) - pad)
ymax = min(ny, int(y_t.max()) + pad + 1)

# 将坐标转换为相对子网格的索引
x_t -= xmin
y_t -= ymin
ys, xs = ymax - ymin, xmax - xmin

print(f"Box裁剪: [{ymin}:{ymax}, {xmin}:{xmax}] ({ys}x{xs}) vs 全网格 {ny}x{nx}")

# 读取子网格数据（保留全部模式面层）
ua = wrf.getvar(nc, "ua")[:, ymin:ymax, xmin:xmax].values
va = wrf.getvar(nc, "va")[:, ymin:ymax, xmin:xmax].values
pres = wrf.getvar(nc, "pressure")[:, ymin:ymax, xmin:xmax].values
nz = ua.shape[0]

# Y轴用中心点的各层气压值
plevs = pres[:, y0 - ymin, x0 - xmin]
print(f"垂直层数: {nz}")

Box裁剪: [62:267, 170:375] (205x205) vs 全网格 437x447
垂直层数: 34

5. 水平插值与切向风/径向风分解

对每一层模式面，使用RegularGridInterpolator在裁剪后的规则网格（y, x）上插值到目标点。相比scipy.interpolate.griddata，RegularGridInterpolator不需要构建三角剖分，在规则网格上每次插值的时间复杂度为O(1)，速度快一个数量级以上。

插值后，用气象方位角公式将笛卡尔坐标下的(u, v)分解为切向风vt和径向风vr：

vt = v * sin(theta) - u * cos(theta)  
vr = u * sin(theta) + v * cos(theta)

其中theta为气象方位角（0°=北，顺时针）。注意该公式与数学极坐标（0°=东，逆时针）的sin/cos组合不同，两者不可混用。

最后对每个半径圈上的所有方位角做平均，得到方位角平均剖面。

pts = np.column_stack([y_t, x_t])
vt_am = np.zeros((nz, NR), dtype=np.float64)
vr_am = np.zeros((nz, NR), dtype=np.float64)

for k in range(nz):
    iu = RegularGridInterpolator(
        (np.arange(ys), np.arange(xs)), ua[k],
        bounds_error=False, fill_value=np.nan
    )
    iv = RegularGridInterpolator(
        (np.arange(ys), np.arange(xs)), va[k],
        bounds_error=False, fill_value=np.nan
    )
    uo = iu(pts).reshape(NAZ, NR)
    vo = iv(pts).reshape(NAZ, NR)

    vt = np.zeros((NAZ, NR))
    vr = np.zeros((NAZ, NR))
    for i in range(NAZ):
        th = az[i] * np.pi / 180.0
        vt[i, :] = vo[i, :] * np.sin(th) - uo[i, :] * np.cos(th)
        vr[i, :] = uo[i, :] * np.sin(th) + vo[i, :] * np.cos(th)

    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)
        vt_am[k, :] = np.nanmean(vt, axis=0)
        vr_am[k, :] = np.nanmean(vr, axis=0)

print(f"切向风范围: {np.nanmin(vt_am):.2f} ~ {np.nanmax(vt_am):.2f} m/s")
print(f"径向风范围: {np.nanmin(vr_am):.2f} ~ {np.nanmax(vr_am):.2f} m/s")

切向风范围: -2.32 ~ 56.58 m/s
径向风范围: -12.48 ~ 10.50 m/s

6. 绘图

使用matplotlib绘制双面板图。左图为切向风vt，右图为径向风vr。纵轴为各模式层在台风中心的气压值（hPa），横轴为半径（m）。

色标范围根据数据自动调整：取vt和vr绝对值的最大值，向上取整到最近的5的倍数，间隔为2 m/s。

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

vmax = max(np.nanmax(np.abs(vt_am)), np.nanmax(np.abs(vr_am)))
cmax = max(10.0, np.ceil(vmax / 5.0) * 5.0)

for ax, fld, ttl in [(ax1, vt_am, "Tangential Wind"), (ax2, vr_am, "Radial Wind")]:
    cf = ax.contourf(r, plevs, fld, levels=10, cmap="bwr", extend="both")
    plt.colorbar(cf, ax=ax, shrink=0.75)
    ax.set_title(ttl)
    ax.set_xlabel("Radius (m)")
    ax.set_ylabel("Pressure (hPa)")
    ax.invert_yaxis()

fig.suptitle(
    f"Center: ({lat_c:.2f}°N, {lon_c:.2f}°E)  SLP={slp[min_idx]:.1f} hPa"
)
plt.tight_layout()
plt.savefig("tc_azimuthal_profile.png", dpi=200, bbox_inches="tight")
plt.show()

output

7. 完整脚本

以下是将上述步骤整合后的完整脚本，约100行，可直接保存为.py文件在命令行运行。

#!/usr/bin/env python3
import sys, warnings
from pathlib import Path
from glob import glob
import numpy as np
from netCDF4 import Dataset
from scipy.interpolate import RegularGridInterpolator
import matplotlib.pyplot as plt
import wrf
import metpy.calc as mpcalc
from metpy.units import units

INPATH = sys.argv[1] if len(sys.argv) > 1else"."
OUTFILE = sys.argv[2] if len(sys.argv) > 2else"tc_azimuthal_profile.png"
MAX_R, NR, NAZ = 300.0, 101, 73

p = Path(INPATH)
if p.is_file(): f = str(p)
elif p.is_dir(): f = (sorted(glob(str(p / "wrfout_d01_*"))) or sorted(glob(str(p / "wrfout_*"))))[0]
else: f = sorted(glob(INPATH))[0]
print(f"读取: {f}")

nc = Dataset(f)
slp = wrf.getvar(nc, "slp").values
min_idx = np.unravel_index(np.argmin(slp), slp.shape)
y0, x0 = min_idx
lat_c = float(wrf.getvar(nc, "lat").values[min_idx])
lon_c = float(wrf.getvar(nc, "lon").values[min_idx])
print(f"中心: {lat_c:.2f}°N, {lon_c:.2f}°E  SLP={slp[min_idx]:.1f} hPa")

az = np.linspace(0, 360, NAZ)
r = np.linspace(0, MAX_R, NR) * 1000
lon_a, lat_a = mpcalc.azimuth_range_to_lat_lon(az * units.degree, r * units.meter, lon_c, lat_c)
x_t, y_t = wrf.ll_to_xy(nc, np.asarray(lat_a).ravel(), np.asarray(lon_a).ravel())
x_t = x_t.values.astype(float)
y_t = y_t.values.astype(float)

ny, nx = len(nc.dimensions["south_north"]), len(nc.dimensions["west_east"])
pad = 3
xmin = max(0, int(x_t.min()) - pad)
xmax = min(nx, int(x_t.max()) + pad + 1)
ymin = max(0, int(y_t.min()) - pad)
ymax = min(ny, int(y_t.max()) + pad + 1)
x_t -= xmin
y_t -= ymin
ys, xs = ymax - ymin, xmax - xmin
print(f"Box裁剪: [{ymin}:{ymax}, {xmin}:{xmax}] ({ys}x{xs}) vs 全网格 {ny}x{nx}")

ua = wrf.getvar(nc, "ua")[:, ymin:ymax, xmin:xmax].values
va = wrf.getvar(nc, "va")[:, ymin:ymax, xmin:xmax].values
pres = wrf.getvar(nc, "pressure")[:, ymin:ymax, xmin:xmax].values
nz = ua.shape[0]
plevs = pres[:, y0 - ymin, x0 - xmin]

pts = np.column_stack([y_t, x_t])
vt_am = np.zeros((nz, NR), dtype=np.float64)
vr_am = np.zeros((nz, NR), dtype=np.float64)
for k in range(nz):
    iu = RegularGridInterpolator((np.arange(ys), np.arange(xs)), ua[k], bounds_error=False, fill_value=np.nan)
    iv = RegularGridInterpolator((np.arange(ys), np.arange(xs)), va[k], bounds_error=False, fill_value=np.nan)
    uo = iu(pts).reshape(NAZ, NR)
    vo = iv(pts).reshape(NAZ, NR)
    vt = np.zeros((NAZ, NR))
    vr = np.zeros((NAZ, NR))
    for i in range(NAZ):
        th = az[i] * np.pi / 180.0
        vt[i, :] = vo[i, :] * np.sin(th) - uo[i, :] * np.cos(th)
        vr[i, :] = uo[i, :] * np.sin(th) + vo[i, :] * np.cos(th)
    with warnings.catch_warnings():
        warnings.simplefilter("ignore", category=RuntimeWarning)
        vt_am[k, :] = np.nanmean(vt, axis=0)
        vr_am[k, :] = np.nanmean(vr, axis=0)

print(f"vt: {np.nanmin(vt_am):.2f} ~ {np.nanmax(vt_am):.2f}  vr: {np.nanmin(vr_am):.2f} ~ {np.nanmax(vr_am):.2f}")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))
vmax = max(np.nanmax(np.abs(vt_am)), np.nanmax(np.abs(vr_am)))
cmax = max(10.0, np.ceil(vmax / 5.0) * 5.0)

for ax, fld, ttl in [(ax1, vt_am, "Tangential Wind"), (ax2, vr_am, "Radial Wind")]:
    cf = ax.contourf(r, plevs, fld, levels=10, cmap="bwr", extend="both")
    plt.colorbar(cf, ax=ax, shrink=0.75)
    ax.set_title(ttl)
    ax.set_xlabel("Radius (m)")
    ax.set_ylabel("Pressure (hPa)")
    ax.invert_yaxis()
fig.suptitle(f"Center: ({lat_c:.2f}°N, {lon_c:.2f}°E)  SLP={slp[min_idx]:.1f} hPa")
plt.tight_layout()
plt.savefig(OUTFILE, dpi=200, bbox_inches="tight")
print(f"输出: {OUTFILE}")
nc.close()

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-05-05，如有侵权请联系 cloudcommunity@tencent.com 删除