首页
学习
活动
专区
圈层
工具
发布
社区首页 >专栏 >代码实战 | 因为metpy函数太慢,一怒之下写了加速版dcape算法

代码实战 | 因为metpy函数太慢,一怒之下写了加速版dcape算法

作者头像
用户11172986
发布2026-05-26 20:43:23
发布2026-05-26 20:43:23
30
举报
文章被收录于专栏:气python风雨气python风雨

代码实战 | 因为metpy函数太慢,一怒之下写了加速版dcape算法

早有前人抱怨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 的实战用法。

一、安装与导入

代码语言:javascript
复制
%pip install skyborn -U 
代码语言:javascript
复制
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.
代码语言:javascript
复制
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
代码语言:javascript
复制
Skyborn version: 0.4.1

二、DCAPE 原理速览

2.1 什么是 DCAPE?

DCAPE 描述的是中层干冷空气下沉到地面过程中可获得的负浮力能量,本质上是 CAPE 的 「镜像」:

  • CAPE → 气块从地面上升到平衡高度释放的正浮力能量
  • DCAPE → 气块从中层(通常 700–500 hPa)下沉到地面释放的负浮力能量

2.2 Skyborn 计算流程

Skyborn 的 calculate_dcape 采用以下步骤(基于 Emanuel 1994 与 compiled backend):

  1. 寻找最干层:在 700–500 hPa 之间搜索 theta-e 最小值所在层
  2. 湿球温度追踪:沿该气块的 假绝热线(pseudoadiabat) 湿球温度下沉到地面
  3. 负浮力积分:计算气块温度与环境温度的差值,从起始高度积分到地面

其中 是沿假绝热线下沉的气块温度, 是环境温度。

2.3 经验阈值参考

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 + 露点温度计算

ERA5 直接提供的变量是 T (K) 和 **q (kg/kg)**,没有露点温度 Td。我们需要先通过 Magnus 公式从比湿反算露点。

输入要求

  • pressure:hPa
  • temperature:°C
  • dewpoint:°C

露点计算流程(Buck 1981,水面):

其中 为水汽压(hPa)。

代码语言:javascript
复制
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]}")
代码语言:javascript
复制
坐标 : ['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
代码语言:javascript
复制
# === 露点温度计算函数 ===
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) 是正常的——那里非常干燥。")
代码语言:javascript
复制
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) 是正常的——那里非常干燥。

四、单点剖面 DCAPE(1D)

先以单点 (40°N, 120°E) 为例,展示完整剖面诊断过程。

calculate_dcape 的 1D 输入返回一个 标量(单位 J/kg)。

代码语言:javascript
复制
# 提取单点剖面
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}")
代码语言:javascript
复制
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

代码语言:javascript
复制
Theta-e minimum in 700–500 hPa: 337.22 K at 550 hPa

五、全网格 DCAPE 空间分布(3D)

把 3D 场 (level, lat, lon) 直接传入 calculate_dcape,即可一次性计算整个区域的 DCAPE。

Skyborn 编译后端处理 241 × 361 × 37 层 数据仅需约 3–5 秒,效率极高。

代码语言:javascript
复制
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}")
代码语言:javascript
复制
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 空间分布 + 关键区域标注

代码语言:javascript
复制
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

计算结果的空间特征

  • 2023-08-02 00 UTC,研究区域内 DCAPE 最高值接近 2000 J/kg,出现在约 25–35°N、120–130°E 附近。
  • 30–40°N, 120–130°E 出现大面积 >800 J/kg 区域。
  • 高纬度(>50°N)DCAPE 普遍较低(<400 J/kg)。

📌 上述描述仅基于本次 ERA5 计算结果的空间分布,不做天气学归因解释。实际业务分析需结合雷达、卫星和探空进一步验证。

六、多时刻演变:DCAPE 的时间变化

利用 ERA5 的 4 小时间隔 数据,我们可以追踪 DCAPE 的日变化和对流系统过境信号。

下面选取 4 个关键时次(00, 08, 12, 16 UTC),快速对比空间格局演变。

代码语言:javascript
复制
# 选取 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.')
代码语言:javascript
复制
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
代码语言:javascript
复制
/tmp/ipykernel_96/751310536.py:12: RuntimeWarning: invalid value encountered in log
  alpha = np.log(e / 6.112)
代码语言:javascript
复制
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.
代码语言:javascript
复制
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 与 500 hPa 风场叠加

DCAPE 的高值区往往与 中层干侵入(dry slot)急流出口区 有关。将 DCAPE 与 500 hPa 风场叠加,可以直观看到这种关联。

代码语言:javascript
复制
# 取 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

叠加观察

  • 蓝色等值线:500 hPa 风速 ≥30 m/s 的区域。
  • 图中 DCAPE 高值带(>800 J/kg)大致位于急流轴南侧(约 30–35°N),与中层干区常出现的位置一致。
  • 风矢量显示该区域为西风急流控制。

📌 上述仅为空间叠加观察,急流与 DCAPE 的因果关系需通过动力诊断(如 Q 矢量、ω 方程)进一步验证。

八、参数速查表

8.1 API

代码语言:javascript
复制
skyborn.calc.calculate_dcape(
代码语言:javascript
复制
pressure,      # hPa, 1D (level,) or 3D (level, lat, lon)  
temperature,   # °C,  同上  
dewpoint,      # °C,  同上  

)

Returns:# 1D input → float (J/kg)# 3D input → np.ndarray (lat, lon) 或 xr.DataArray

8.2 单位换算 Checklist

ERA5 原始单位

calculate_dcape 要求

转换

气压

hPa (level coord)

hPa

无需转换

温度

K

°C

T - 273.15

露点

无直接变量

°C

从 q 或 r 计算

8.3 露点计算(从比湿 q)

代码语言:javascript
复制
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

8.4 常见陷阱

现象

原因

解决

返回 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 为强对流诊断提供了 高效、易用 的工具:

核心能力

  • 支持 1D 单点剖面(返回标量)和 3D 全网格(返回 2D 场)
  • Compiled backend:241×361×37 网格 ~4 秒 完成计算
  • 自动寻找 700–500 hPa theta-e 最小值层,沿假绝热线下沉积分

最佳实践

  1. 单位三件套:p (hPa)、T (°C)、Td (°C) —— 这是最常见的出错点
  2. 露点预计算:ERA5 无直接 Td,需从 q 经 Magnus 公式反算
  3. 空间分析:DCAPE >800 J/kg 区域在业务中常被关注,但具体阈值需本地化标定
  4. 时间演变:利用再分析数据追踪 DCAPE 随时间变化,识别潜势增强时段

关于 DCAPE 阈值的说明

教程中使用的分级(<400 / 400–800 / 800–1200 / >1200 J/kg)为经验参考区间,主要参考来源:

  • Emanuel, K. A. (1994). Atmospheric Convection. Oxford University Press.(DCAPE 原始定义)
  • 业务预报中常将 DCAPE >1000 J/kg 作为强下沉气流的参考指标(NOAA/SPC 培训材料)
  • 具体业务阈值应结合本地探空资料、雷达回波和灾情记录进行标定

教程数据:ERA5 2023-08,© ECMWF,CDS 开放数据 Skyborn 0.40 · 2026

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2026-05-26,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 气python风雨 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 代码实战 | 因为metpy函数太慢,一怒之下写了加速版dcape算法
    • 一、安装与导入
    • 二、DCAPE 原理速览
      • 2.1 什么是 DCAPE?
      • 2.2 Skyborn 计算流程
      • 2.3 经验阈值参考
    • 三、数据准备:ERA5 + 露点温度计算
    • 四、单点剖面 DCAPE(1D)
    • 五、全网格 DCAPE 空间分布(3D)
    • 六、多时刻演变:DCAPE 的时间变化
    • 七、DCAPE 与 500 hPa 风场叠加
    • 八、参数速查表
      • 8.1 API
  • Returns:# 1D input → float (J/kg)# 3D input → np.ndarray (lat, lon) 或 xr.DataArray
    • 8.2 单位换算 Checklist
    • 8.3 露点计算(从比湿 q)
    • 8.4 常见陷阱
    • 九、小结
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档