早有前人抱怨metpy的cape等函数仅支持一维数据。开发者称正在研究如何支持多维
但是小编不想等下去了,于是
Skyborn 0.4 新增 DCAPE (Downdraft Convective Available Potential Energy, 下沉对流有效位能) 计算模块,基于编译后端实现,支持 1D 单点剖面 与 3D 全网格 两种输入模式。
DCAPE 是衡量对流风暴中 冷池强度 / 下沉气流能量 的重要指标,对短时强降水、雷暴大风、下击暴流等强对流天气的预警具有重要参考价值。
本教程使用 ERA5 (2023-08-02 ~ 2023-08-05) 再分析数据,从 露点温度计算 → 单点 DCAPE → 全网格空间分布 → 多时刻演变分析,完整演示 skyborn.calc.calculate_dcape 的实战用法。
%pip install skyborn -U
Requirement already satisfied: skyborn in /opt/conda/lib/python3.11/site-packages (0.4.1)
Requirement already satisfied: numpy in /opt/conda/lib/python3.11/site-packages (from skyborn) (1.26.4)
Requirement already satisfied: pandas in /opt/conda/lib/python3.11/site-packages (from skyborn) (2.2.3)
Requirement already satisfied: xarray in /opt/conda/lib/python3.11/site-packages (from skyborn) (2024.3.0)
Requirement already satisfied: matplotlib in /opt/conda/lib/python3.11/site-packages (from skyborn) (3.9.1)
Requirement already satisfied: netCDF4 in /opt/conda/lib/python3.11/site-packages (from skyborn) (1.6.3)
Requirement already satisfied: metpy in /opt/conda/lib/python3.11/site-packages (from skyborn) (1.6.3)
Requirement already satisfied: tqdm in /opt/conda/lib/python3.11/site-packages (from skyborn) (4.67.0)
Requirement already satisfied: statsmodels in /opt/conda/lib/python3.11/site-packages (from skyborn) (0.14.4)
Requirement already satisfied: scipy in /opt/conda/lib/python3.11/site-packages (from skyborn) (1.14.1)
Requirement already satisfied: scikit-learn in /opt/conda/lib/python3.11/site-packages (from skyborn) (1.6.0)
Requirement already satisfied: dask in /opt/conda/lib/python3.11/site-packages (from skyborn) (2024.8.1)
Requirement already satisfied: click>=8.1 in /opt/conda/lib/python3.11/site-packages (from dask->skyborn) (8.1.7)
Requirement already satisfied: cloudpickle>=3.0.0 in /opt/conda/lib/python3.11/site-packages (from dask->skyborn) (3.1.0)
Requirement already satisfied: fsspec>=2021.09.0 in /opt/conda/lib/python3.11/site-packages (from dask->skyborn) (2025.2.0)
Requirement already satisfied: packaging>=20.0 in /opt/conda/lib/python3.11/site-packages (from dask->skyborn) (24.1)
Requirement already satisfied: partd>=1.4.0 in /opt/conda/lib/python3.11/site-packages (from dask->skyborn) (1.4.2)
Requirement already satisfied: pyyaml>=5.3.1 in /opt/conda/lib/python3.11/site-packages (from dask->skyborn) (6.0.2)
Requirement already satisfied: toolz>=0.10.0 in /opt/conda/lib/python3.11/site-packages (from dask->skyborn) (1.0.0)
Requirement already satisfied: importlib-metadata>=4.13.0 in /opt/conda/lib/python3.11/site-packages (from dask->skyborn) (8.5.0)
Requirement already satisfied: contourpy>=1.0.1 in /opt/conda/lib/python3.11/site-packages (from matplotlib->skyborn) (1.3.1)
Requirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.11/site-packages (from matplotlib->skyborn) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /opt/conda/lib/python3.11/site-packages (from matplotlib->skyborn) (4.55.3)
Requirement already satisfied: kiwisolver>=1.3.1 in /opt/conda/lib/python3.11/site-packages (from matplotlib->skyborn) (1.4.7)
Requirement already satisfied: pillow>=8 in /opt/conda/lib/python3.11/site-packages (from matplotlib->skyborn) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /opt/conda/lib/python3.11/site-packages (from matplotlib->skyborn) (3.2.0)
Requirement already satisfied: python-dateutil>=2.7 in /opt/conda/lib/python3.11/site-packages (from matplotlib->skyborn) (2.9.0.post0)
Requirement already satisfied: pint>=0.17 in /opt/conda/lib/python3.11/site-packages (from metpy->skyborn) (0.24.4)
Requirement already satisfied: pooch>=1.2.0 in /opt/conda/lib/python3.11/site-packages (from metpy->skyborn) (1.8.2)
Requirement already satisfied: pyproj>=3.0.0 in /opt/conda/lib/python3.11/site-packages (from metpy->skyborn) (3.5.0)
Requirement already satisfied: traitlets>=5.0.5 in /opt/conda/lib/python3.11/site-packages (from metpy->skyborn) (5.14.3)
Requirement already satisfied: pytz>=2020.1 in /opt/conda/lib/python3.11/site-packages (from pandas->skyborn) (2024.1)
Requirement already satisfied: tzdata>=2022.7 in /opt/conda/lib/python3.11/site-packages (from pandas->skyborn) (2024.2)
Requirement already satisfied: cftime in /opt/conda/lib/python3.11/site-packages (from netCDF4->skyborn) (1.6.4)
Requirement already satisfied: joblib>=1.2.0 in /opt/conda/lib/python3.11/site-packages (from scikit-learn->skyborn) (1.4.2)
Requirement already satisfied: threadpoolctl>=3.1.0 in /opt/conda/lib/python3.11/site-packages (from scikit-learn->skyborn) (3.5.0)
Requirement already satisfied: patsy>=0.5.6 in /opt/conda/lib/python3.11/site-packages (from statsmodels->skyborn) (1.0.1)
Requirement already satisfied: zipp>=3.20 in /opt/conda/lib/python3.11/site-packages (from importlib-metadata>=4.13.0->dask->skyborn) (3.21.0)
Requirement already satisfied: locket in /opt/conda/lib/python3.11/site-packages (from partd>=1.4.0->dask->skyborn) (1.0.0)
Requirement already satisfied: platformdirs>=2.1.0 in /opt/conda/lib/python3.11/site-packages (from pint>=0.17->metpy->skyborn) (4.3.6)
Requirement already satisfied: typing_extensions>=4.0.0 in /opt/conda/lib/python3.11/site-packages (from pint>=0.17->metpy->skyborn) (4.12.2)
Requirement already satisfied: flexcache>=0.3 in /opt/conda/lib/python3.11/site-packages (from pint>=0.17->metpy->skyborn) (0.3)
Requirement already satisfied: flexparser>=0.4 in /opt/conda/lib/python3.11/site-packages (from pint>=0.17->metpy->skyborn) (0.4)
Requirement already satisfied: requests>=2.19.0 in /opt/conda/lib/python3.11/site-packages (from pooch>=1.2.0->metpy->skyborn) (2.32.3)
Requirement already satisfied: certifi in /opt/conda/lib/python3.11/site-packages (from pyproj>=3.0.0->metpy->skyborn) (2024.12.14)
Requirement already satisfied: six>=1.5 in /opt/conda/lib/python3.11/site-packages (from python-dateutil>=2.7->matplotlib->skyborn) (1.17.0)
Requirement already satisfied: charset-normalizer<4,>=2 in /opt/conda/lib/python3.11/site-packages (from requests>=2.19.0->pooch>=1.2.0->metpy->skyborn) (3.4.0)
Requirement already satisfied: idna<4,>=2.5 in /opt/conda/lib/python3.11/site-packages (from requests>=2.19.0->pooch>=1.2.0->metpy->skyborn) (3.10)
Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/conda/lib/python3.11/site-packages (from requests>=2.19.0->pooch>=1.2.0->metpy->skyborn) (2.2.3)
Note: you may need to restart the kernel to use updated packages.
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import skyborn as skb
import skyborn.calc as sc
print(f"Skyborn version: {skb.__version__}")
plt.rcParams['font.family'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
Skyborn version: 0.4.1
DCAPE 描述的是中层干冷空气下沉到地面过程中可获得的负浮力能量,本质上是 CAPE 的 「镜像」:
Skyborn 的 calculate_dcape 采用以下步骤(基于 Emanuel 1994 与 compiled backend):
其中 是沿假绝热线下沉的气块温度, 是环境温度。
DCAPE 范围 | 一般经验描述 |
|---|---|
< 400 J/kg | 弱下沉能量 |
400–800 J/kg | 中等能量 |
800–1200 J/kg | 较强能量 |
> 1200 J/kg | 强下沉能量 |
⚠️ 上述分级仅为经验参考区间,非严格业务标准。主要参考:Emanuel (1994) 原始定义;NOAA/SPC 培训材料中常将 DCAPE >1000 J/kg 作为强下沉气流关注指标。实际预报需结合垂直风切变、0–3 km CAPE 等参数综合判断,本地化阈值应通过本地探空和灾情记录标定。
ERA5 直接提供的变量是 T (K) 和 **q (kg/kg)**,没有露点温度 Td。我们需要先通过 Magnus 公式从比湿反算露点。
输入要求:
pressure:hPatemperature:°Cdewpoint:°C露点计算流程(Buck 1981,水面):
其中 为水汽压(hPa)。
DATA_FILE = '/home/mw/input/era58091/ERA5-2023-08_pl.nc'
ds = xr.open_dataset(DATA_FILE)
print('坐标 :', list(ds.coords))
print('变量 :', list(ds.data_vars))
print(f"层次 : {ds.level.values}")
print(f"时间 : {ds.time.values[0]} → {ds.time.values[-1]}")
坐标 : ['longitude', 'latitude', 'level', 'time']
变量 : ['z', 'r', 't', 'q', 'u', 'v', 'w']
层次 : [ 1 2 3 5 7 10 20 30 50 70 100 125 150 175
200 225 250 300 350 400 450 500 550 600 650 700 750 775
800 825 850 875 900 925 950 975 1000]
时间 : 2023-08-02T00:00:00.000000000 → 2023-08-05T16:00:00.000000000
# === 露点温度计算函数 ===
def calc_dewpoint_from_q(pressure_hPa, q_kgkg):
"""
从比湿 q 和压力 p 计算露点温度 (°C).
pressure_hPa: hPa, 可广播
q_kgkg : kg/kg, 比湿
"""
epsilon = 0.622
# 水汽压 e (hPa)
e = pressure_hPa * q_kgkg / (epsilon + q_kgkg)
# Magnus 反算 (Buck 1981, over water)
alpha = np.log(e / 6.112)
td = (243.5 * alpha) / (17.67 - alpha)
return td
# 取分析时刻
ds_t0 = ds.isel(time=0)
# 构造 3D 压力场 (hPa), 与 T/q 同 shape
p_3d = ds_t0['level'].values.astype('float64')[:, None, None]
p_3d = np.broadcast_to(p_3d, ds_t0['t'].shape)
# 温度 (°C) 和 比湿
t_3d = ds_t0['t'].values - 273.15 # K → °C
q_3d = ds_t0['q'].values # kg/kg
# 计算露点
td_3d = calc_dewpoint_from_q(p_3d, q_3d)
print(f"p_3d shape: {p_3d.shape}, range: {p_3d.min():.0f}–{p_3d.max():.0f} hPa")
print(f"t_3d shape: {t_3d.shape}, range: {t_3d.min():.1f}–{t_3d.max():.1f} °C")
print(f"td_3d shape: {td_3d.shape}, range: {td_3d.min():.1f}–{td_3d.max():.1f} °C")
print(f"\n注意: 高空露点很低 (<< -50°C) 是正常的——那里非常干燥。")
p_3d shape: (37, 241, 361), range: 1–1000 hPa
t_3d shape: (37, 241, 361), range: -85.6–37.6 °C
td_3d shape: (37, 241, 361), range: -107.0–27.5 °C
注意: 高空露点很低 (<< -50°C) 是正常的——那里非常干燥。
先以单点 (40°N, 120°E) 为例,展示完整剖面诊断过程。
calculate_dcape 的 1D 输入返回一个 标量(单位 J/kg)。
# 提取单点剖面
loc = dict(latitude=40, longitude=120, method='nearest')
p_prof = ds_t0['level'].values.astype('float64')
t_prof = ds_t0['t'].sel(**loc).values - 273.15
q_prof = ds_t0['q'].sel(**loc).values
td_prof = calc_dewpoint_from_q(p_prof, q_prof)
print(f"Profile @ (40°N, 120°E), {ds_t0.time.values}")
print(f" Levels: {len(p_prof)} ({p_prof.min():.0f} – {p_prof.max():.0f} hPa)")
print(f" T : {t_prof.min():.1f} – {t_prof.max():.1f} °C")
print(f" Td: {td_prof.min():.1f} – {td_prof.max():.1f} °C")
# 计算 DCAPE
dcape_val = sc.calculate_dcape(p_prof, t_prof, td_prof)
print(f"\n>>> DCAPE = {dcape_val:.2f} J/kg")
if dcape_val < 400:
level = 'Weak'
elif dcape_val < 800:
level = 'Moderate'
elif dcape_val < 1200:
level = 'Strong — gust potential'
else:
level = 'Very strong — downburst risk'
print(f" Level: {level}")
Profile @ (40°N, 120°E), 2023-08-02T00:00:00.000000000
Levels: 37 (1 – 1000 hPa)
T : -69.5 – 27.1 °C
Td: -106.6 – 24.5 °C
>>> DCAPE = 1001.16 J/kg
Level: Strong — gust potential
绘图:T / Td / Theta-e 剖面,直观理解 DCAPE 来源
DCAPE 的核心是 700–500 hPa 之间 theta-e 最小值。我们叠加 theta-e 剖面来展示这一点。
output
Theta-e minimum in 700–500 hPa: 337.22 K at 550 hPa
把 3D 场 (level, lat, lon) 直接传入 calculate_dcape,即可一次性计算整个区域的 DCAPE。
Skyborn 编译后端处理 241 × 361 × 37 层 数据仅需约 3–5 秒,效率极高。
import time
print("Running calculate_dcape on full 3D grid...")
t0 = time.time()
dcape_2d = sc.calculate_dcape(p_3d, t_3d, td_3d)
t1 = time.time()
print(f"Done in {t1-t0:.2f} s")
print(f"Output shape: {dcape_2d.shape}")
print(f"DCAPE range : {np.nanmin(dcape_2d):.2f} – {np.nanmax(dcape_2d):.2f} J/kg")
print(f"Mean DCAPE : {np.nanmean(dcape_2d):.2f} J/kg")
print(f"NaN count : {np.sum(np.isnan(dcape_2d))} / {dcape_2d.size}")
Running calculate_dcape on full 3D grid...
Done in 4.63 s
Output shape: (241, 361)
DCAPE range : 0.00 – 1983.52 J/kg
Mean DCAPE : 587.29 J/kg
NaN count : 0 / 87001
绘图:DCAPE 空间分布 + 关键区域标注
fig, ax = plt.subplots(figsize=(12, 6))
# 填色
levels = [0, 200, 400, 600, 800, 1000, 1200, 1500, 2000]
cf = ax.contourf(ds.longitude, ds.latitude, dcape_2d,
levels=levels, cmap='YlOrRd', extend='max')
# 等高线
cs = ax.contour(ds.longitude, ds.latitude, dcape_2d,
levels=[400, 800, 1200], colors='k', linewidths=0.8)
ax.clabel(cs, inline=True, fontsize=9, fmt='%d')
# 标记 (40°N, 120°E)
ax.plot(120, 40, 'b*', ms=15, markeredgecolor='white', markeredgewidth=0.5,
label=f'Profile point: {dcape_val:.0f} J/kg')
ax.set_title(f'DCAPE Distribution · {str(ds_t0.time.values)[:16]} UTC')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_xlim(90, 180)
ax.set_ylim(10, 70)
plt.colorbar(cf, ax=ax, label='DCAPE (J kg⁻¹)')
ax.legend(loc='upper left')
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
output
计算结果的空间特征:
📌 上述描述仅基于本次 ERA5 计算结果的空间分布,不做天气学归因解释。实际业务分析需结合雷达、卫星和探空进一步验证。
利用 ERA5 的 4 小时间隔 数据,我们可以追踪 DCAPE 的日变化和对流系统过境信号。
下面选取 4 个关键时次(00, 08, 12, 16 UTC),快速对比空间格局演变。
# 选取 4 个时刻 (ERA5 为 4 小时间隔)
time_indices = [0, 2, 4, 6]
times_sel = [str(ds.time.values[i]) for i in time_indices]
print('Selected times:', times_sel)
dcape_series = {}
for idx, t_str in zip(time_indices, times_sel):
ds_t = ds.isel(time=idx)
p = ds_t['level'].values.astype('float64')[:, None, None]
p = np.broadcast_to(p, ds_t['t'].shape)
t = ds_t['t'].values - 273.15
q = ds_t['q'].values
td = calc_dewpoint_from_q(p, q)
dc = sc.calculate_dcape(p, t, td)
dcape_series[t_str] = dc
print(f'{t_str[11:16]}: DCAPE max = {np.nanmax(dc):.0f} J/kg, mean = {np.nanmean(dc):.0f} J/kg')
print('\nAll 4 times computed successfully.')
Selected times: ['2023-08-02T00:00:00.000000000', '2023-08-02T08:00:00.000000000', '2023-08-02T16:00:00.000000000', '2023-08-03T04:00:00.000000000']
00:00: DCAPE max = 1984 J/kg, mean = 587 J/kg
08:00: DCAPE max = 2558 J/kg, mean = 680 J/kg
/tmp/ipykernel_96/751310536.py:12: RuntimeWarning: invalid value encountered in log
alpha = np.log(e / 6.112)
16:00: DCAPE max = 2345 J/kg, mean = 614 J/kg
04:00: DCAPE max = 2225 J/kg, mean = 634 J/kg
All 4 times computed successfully.
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
levels = [0, 200, 400, 600, 800, 1000, 1200, 1500, 2000]
for ax, t_str in zip(axes.flat, times_sel):
dc = dcape_series[t_str]
cf = ax.contourf(ds.longitude, ds.latitude, dc,
levels=levels, cmap='YlOrRd', extend='max')
cs = ax.contour(ds.longitude, ds.latitude, dc,
levels=[400, 800, 1200], colors='k', linewidths=0.6)
ax.clabel(cs, inline=True, fontsize=8, fmt='%d')
ax.set_title(f'{t_str[11:16]} UTC · max={np.nanmax(dc):.0f} J/kg')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_xlim(90, 180)
ax.set_ylim(10, 70)
plt.colorbar(cf, ax=ax, label='DCAPE (J/kg)')
ax.grid(alpha=0.3)
plt.suptitle('Skyborn 0.40 · DCAPE Diurnal Evolution (2023-08-02)',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
output
各时次计算结果对比:
时次 | DCAPE 最大值 | DCAPE 平均值 | 高值区位置变化 |
|---|---|---|---|
00 UTC | ~1984 J/kg | ~587 J/kg | 高值区主要位于 25–35°N |
08 UTC | ~2558 J/kg | ~680 J/kg | 高值区向北扩展至 30–40°N |
12 UTC | ~2345 J/kg | ~614 J/kg | 高值区东移,覆盖黄海–日本海 |
16 UTC | ~2225 J/kg | ~634 J/kg | 高值区略有南压 |
📌 各时次差异反映了 DCAPE 随时间的演变,具体原因需结合当日天气形势(高空槽、副高位置、低空急流等)综合分析,此处不做归因。
DCAPE 的高值区往往与 中层干侵入(dry slot) 和 急流出口区 有关。将 DCAPE 与 500 hPa 风场叠加,可以直观看到这种关联。
# 取 500 hPa 风场
u500 = ds_t0['u'].sel(level=500)
v500 = ds_t0['v'].sel(level=500)
speed500 = np.sqrt(u500**2 + v500**2)
fig, ax = plt.subplots(figsize=(13, 6.5))
# DCAPE 填色
levels = [0, 200, 400, 600, 800, 1000, 1200, 1500, 2000]
cf = ax.contourf(ds.longitude, ds.latitude, dcape_2d,
levels=levels, cmap='YlOrRd', extend='max')
# 500 hPa 风速等值线
cs = ax.contour(ds.longitude, ds.latitude, speed500,
levels=[20, 30, 40, 50], colors='blue', linewidths=0.7, alpha=0.7)
ax.clabel(cs, inline=True, fontsize=9, fmt='%d', colors='blue')
# 风矢量(抽稀)
skip = 8
ax.quiver(ds.longitude.values[::skip], ds.latitude.values[::skip],
u500.values[::skip, ::skip], v500.values[::skip, ::skip],
scale=400, width=0.002, color='navy', alpha=0.7)
ax.set_title('DCAPE + 500 hPa Wind Speed & Vector')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
plt.colorbar(cf, ax=ax, label='DCAPE (J/kg)')
ax.set_xlim(90, 180)
ax.set_ylim(10, 70)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
output
叠加观察:
📌 上述仅为空间叠加观察,急流与 DCAPE 的因果关系需通过动力诊断(如 Q 矢量、ω 方程)进一步验证。
skyborn.calc.calculate_dcape(
pressure, # hPa, 1D (level,) or 3D (level, lat, lon)
temperature, # °C, 同上
dewpoint, # °C, 同上
)
量 | ERA5 原始单位 | calculate_dcape 要求 | 转换 |
|---|---|---|---|
气压 | hPa (level coord) | hPa | 无需转换 |
温度 | K | °C | T - 273.15 |
露点 | 无直接变量 | °C | 从 q 或 r 计算 |
epsilon = 0.622
e = pressure_hPa * q_kgkg / (epsilon + q_kgkg) # 水汽压 (hPa) alpha = np.log(e / 6.112) Td = (243.5 * alpha) / (17.67 - alpha) # °C
现象 | 原因 | 解决 |
|---|---|---|
返回 NaN | 有效层数 <4,或 500–700 hPa 无 theta-e 最小值 | 检查 T/Td 是否有缺失,确保包含 500–700 hPa |
全图 NaN | 压力单位错误(用了 Pa 而非 hPa) | 确认 pressure 为 hPa |
DCAPE 异常高 (>3000) | 温度单位错误(用了 K 而非 °C) | T = T_K - 273.15 |
3D 输入报 shape mismatch | p/T/Td 维度不一致 | 用 np.broadcast_to 统一 shape |
计算慢 | 用纯 Python 循环而非 compiled backend | 确保安装的是编译版 skyborn |
Skyborn 0.40 的 calculate_dcape 为强对流诊断提供了 高效、易用 的工具:
核心能力
最佳实践
q 经 Magnus 公式反算关于 DCAPE 阈值的说明
教程中使用的分级(<400 / 400–800 / 800–1200 / >1200 J/kg)为经验参考区间,主要参考来源:
教程数据:ERA5 2023-08,© ECMWF,CDS 开放数据 Skyborn 0.40 · 2026