WRF | 如何基于wrfout绘制Stüve（施蒂韦）图（内含三种探空图代码）

用户11172986

发布于 2026-04-24 19:11:33

420

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

WRF | 如何基于wrfout绘制Stüve（施蒂韦）图（内含三种探空图代码）

前言

有读者提出对于艾玛图与skew图的区别傻傻分不清。本期就再增加一个Stüve图

并将三个图放一起对比，并修复上期作图的bug（cape填色显示错误）。

Stüve图与Skew-T逻辑相似，但纵坐标使用 p^(R/cp) 线性化，坐标轴垂直，读图更直观。

WRF模式输出的wrfout文件再一次成为我们的“金矿”。

本文给出一条从wrfout到出版级Stüve图的“一键流水线”，依旧依赖wrf-python + MetPy 1.7，代码不足80行，复制即可跑通。

!pip install --upgrade metpy -i https://pypi.mirrors.ustc.edu.cn/simple/

Stuve

import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from wrf import getvar, ll_to_xy
import metpy.calc as mpcalc
from metpy.plots import Stuve
from metpy.units import units

# 1. 读数据
ncfile = Dataset("/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_19_00_00")
lat0, lon0 = 25.0, 120.4
x, y = ll_to_xy(ncfile, lat0, lon0)

p  = getvar(ncfile, "pressure")[:, y, x].values * units.hPa
tc = getvar(ncfile, "tc")[:, y, x].values       * units.degC
td = getvar(ncfile, "td")[:, y, x].values       * units.degC
u  = getvar(ncfile, "ua", units="kt")[:, y, x].values * units.knots
v  = getvar(ncfile, "va", units="kt")[:, y, x].values * units.knots

# 2. 计算
prof = mpcalc.parcel_profile(p, tc[0], td[0])
lcl_p, lcl_t = mpcalc.lcl(p[0], tc[0], td[0])
cape, cin = mpcalc.surface_based_cape_cin(p, tc, td)

# 3. 绘图
fig = plt.figure(figsize=(8, 8))
stuve = Stuve(fig)

stuve.plot(p, tc, color='tab:red', linewidth=2, label='Temperature')
stuve.plot(p, td, color='tab:green', linewidth=2, label='Dewpoint')
stuve.plot(p, prof, color='k', linewidth=2, linestyle='--', label='Parcel')

stuve.shade_cape(p, tc, prof, alpha=0.3, color='orange')  # 不再报错
stuve.shade_cin(p, tc, prof, alpha=0.3, color='gray')

idx = np.arange(0, p.size, 3)
stuve.plot_barbs(p[idx], u[idx], v[idx])

stuve.plot_dry_adiabats(alpha=0.3, color='orangered')
stuve.plot_moist_adiabats(alpha=0.3, color='royalblue')
stuve.plot_mixing_lines(alpha=0.3, color='purple')

stuve.plot(lcl_p, lcl_t, 'ko', markersize=8)

stuve.ax.set_ylim(1050, 100)
stuve.ax.set_xlim(-110, 50)
plt.title(f'Stüve diagram @ ({lat0}, {lon0})', loc='left')
plt.title(f'CAPE={cape.m:.0f} J/kg  CIN={cin.m:.0f} J/kg', loc='right')
plt.legend()
plt.tight_layout()
plt.show()

image

skew

import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from wrf import getvar, ll_to_xy
import metpy.calc as mpcalc
from metpy.plots import SkewT
from metpy.units import units

# 数据路径
ncfile = Dataset("/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_19_00_00")
lat0, lon0 = 25.0, 120.4
x, y = ll_to_xy(ncfile, lat0, lon0)

# ===== 关键修复：加 .values 去掉 xarray 壳 =====
p  = getvar(ncfile, "pressure")[:, y, x].values * units.hPa
tc = getvar(ncfile, "tc")[:, y, x].values       * units.degC
td = getvar(ncfile, "td")[:, y, x].values       * units.degC
u  = getvar(ncfile, "ua", units="kt")[:, y, x].values * units.knots
v  = getvar(ncfile, "va", units="kt")[:, y, x].values * units.knots

# 热力学计算
prof = mpcalc.parcel_profile(p, tc[0], td[0])
lcl_p, lcl_t = mpcalc.lcl(p[0], tc[0], td[0])
cape, cin = mpcalc.surface_based_cape_cin(p, tc, td)

# 绘图
fig = plt.figure(figsize=(9, 9))
skew = SkewT(fig, rotation=45)

skew.plot(p, tc, 'r',  lw=2, label='T')
skew.plot(p, td, 'g',  lw=2, label='Td')
skew.plot(p, prof, 'k', lw=2, ls='--', label='Parcel')

skew.shade_cape(p, tc, prof, alpha=0.3, color='orange') 
skew.shade_cin(p, tc, prof,  alpha=0.3, color='gray')

idx = np.arange(0, p.size, 3)
skew.plot_barbs(p[idx], u[idx], v[idx])

skew.plot_dry_adiabats(alpha=0.25, color='orangered')
skew.plot_moist_adiabats(alpha=0.25, color='royalblue')
skew.plot_mixing_lines(alpha=0.25, color='purple')

skew.plot(lcl_p, lcl_t, 'ko', markersize=8)

skew.ax.set_ylim(1050, 100)
skew.ax.set_xlim(-110, 50)
plt.title(f'Skew-T log-P @ ({lat0}, {lon0})', loc='left')
plt.title(f'CAPE={cape.m:.0f} J/kg  CIN={cin.m:.0f} J/kg', loc='right')
plt.legend()
plt.tight_layout()
plt.show()

image

Emagram

import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from wrf import getvar, ll_to_xy
import metpy.calc as mpcalc
from metpy.plots import Emagram   # ① 关键：引入 Emagram
from metpy.units import units

# 1. 读数据
ncfile = Dataset("/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_19_00_00")
lat0, lon0 = 25.0, 120.4
x, y = ll_to_xy(ncfile, lat0, lon0)

# ===== 👇 关键修复：加 .values 去掉 xarray 壳 👇 =====
p  = getvar(ncfile, "pressure")[:, y, x].values * units.hPa
tc = getvar(ncfile, "tc")[:, y, x].values       * units.degC
td = getvar(ncfile, "td")[:, y, x].values       * units.degC
u  = getvar(ncfile, "ua", units="kt")[:, y, x].values * units.knots
v  = getvar(ncfile, "va", units="kt")[:, y, x].values * units.knots

# 2. 计算热力学量
prof = mpcalc.parcel_profile(p, tc[0], td[0])
lcl_p, lcl_t = mpcalc.lcl(p[0], tc[0], td[0])
cape, cin = mpcalc.surface_based_cape_cin(p, tc, td)

# 3. 绘图（Emagram）
fig = plt.figure(figsize=(8, 8))
ema = Emagram(fig)

# 基本曲线
ema.plot(p, tc, color='tab:red', linewidth=2, label='Temperature')
ema.plot(p, td, color='tab:green', linewidth=2, label='Dewpoint')
ema.plot(p, prof, color='k', linewidth=2, linestyle='--', label='Parcel')

# 能量区填色
ema.shade_cape(p, tc, prof, alpha=0.3, color='orange')
ema.shade_cin(p, tc, prof, alpha=0.3, color='gray')

# 风羽（隔3层）
idx = np.arange(0, p.size, 3)
ema.plot_barbs(p[idx], u[idx], v[idx])

# 参考线
ema.plot_dry_adiabats(alpha=0.25, color='orangered')
ema.plot_moist_adiabats(alpha=0.25, color='royalblue')
ema.plot_mixing_lines(alpha=0.25, color='purple')

# 标记 LCL
ema.plot(lcl_p, lcl_t, 'ko', markersize=8)

# 坐标范围 & 标题
ema.ax.set_ylim(1050, 100)
ema.ax.set_xlim(-110, 50)
plt.title(f'Emagram @ ({lat0}, {lon0})', loc='left')
plt.title(f'CAPE={cape.m:.0f} J/kg  CIN={cin.m:.0f} J/kg', loc='right')
plt.legend()
plt.tight_layout()
plt.show()