乐于分享
好东西不私藏
夜雨聆风 > > 办公文件 > 代码实战 | 下载 ERA5 并计算 GPI(台风最大潜在强度)

代码实战 | 下载 ERA5 并计算 GPI(台风最大潜在强度)

当前时间： 2026-08-17 01:49:37 分类：办公文件 评论(0)
代码实战 | 下载 ERA5 并计算 GPI(台风最大潜在强度)

代码实战 | 下载 ERA5 并计算 GPI(台风最大潜在强度)

前言

如何计算 GPI 在之前的 skyborn 库有介绍，之前的示例对于 era5 的变量名有小小的瑕疵。

本期文章给出一份"开箱即用"的完整脚本，实现功能如下：

  • 自动从 CDS（Copernicus Data Store）下载 ERA5 单层/多层变量；
  • 把下载好的文件拼接成一场台风研究常用的 0.25°、逐小时或逐日数据；
  • 用 Skyborn（xarray 接口）一次性算出整个区域/时段的台风最大潜在强度（PI）。

脚本已经在 Linux / macOS / Win 上 Python≥3.9 实测通过，只需把 dates、area 改成自己需要的时间段/区域即可直接跑。

0. 环境准备

# 建议先建个干净环境conda create -n era5_PI python=3.11 -yconda activate era5_PI# 一次性装好pip install cdsapi skyborn xarray netcdf4 dask

1. 在 CDS 注册并获取 API key

  1. 注册 https://cds.climate.copernicus.eu/
  2. 登录后右上角 "My API key" -> 复制两行 url+key
  3. 写入本地文件 ~/.cdsapirc（Win 用户路径为 %USERPROFILE%\.cdsapirc）
url: https://cds.climate.copernicus.eu/api/v2key: 你的uid:你的key

2. 下载脚本：get_era5_for_PI.py

#!/usr/bin/env python"""一键下载 ERA5 计算台风 PI 所需变量单层：sst, msl多层：t, q  (37 层，1000 → 1 hPa，Skyborn 官方推荐)运行前请先在 ~/.cdsapirc 里放好 CDS API key"""import os, cdsapi, datetime as dtfrom pathlib import Path# ---------------- 用户区域 ----------------dates     = ["2023-08-01/to/2023-08-31"]   # 可写多段area      = [30, 110, 15, 140]            # N/W/S/Eoutdir    = Path("ERA5")# -----------------------------------------outdir.mkdir(exist_ok=True)c = cdsapi.Client()# 单层single = ['sea_surface_temperature', 'mean_sea_level_pressure']c.retrieve("reanalysis-era5-single-levels",    {"product_type": "reanalysis","variable"    : single,"date"        : dates,"area"        : area,"grid"        : "0.25/0.25","format"      : "netcdf","time"        : "00/to/23/by/1",   # 逐小时；如只需日平均改成 ["00:00", "12:00"]    },    outdir / "ERA5_single.nc")# 多层levels = [1000, 975, 950, 925, 900, 875, 850, 825, 800, 775, 750, 700,650, 600, 550, 500, 450, 400, 350, 300, 250, 200,175, 150, 125, 100, 70, 50, 30, 20, 10, 7, 5, 3, 2, 1]multi = ['temperature', 'specific_humidity']c.retrieve("reanalysis-era5-pressure-levels",    {"product_type"   : "reanalysis","variable"       : multi,"pressure_level" : levels,"date"           : dates,"area"           : area,"grid"           : "0.25/0.25","format"         : "netcdf","time"           : "00/to/23/by/1",    },    outdir / "ERA5_multilevel.nc")print("✅ 下载完成！文件位于", outdir.absolute())

运行

python get_era5_for_PI.py

即可得到

ERA5/ ├─ ERA5_single.nc └─ ERA5_multilevel.nc

3. 计算脚本：calc_PI_ERA5.py

#!/usr/bin/env python"""用 Skyborn 计算 ERA5 的台风最大潜在强度 PI支持 1D/3D/4D 场，自动识别单位 & 缺测"""import xarray as xrfrom skyborn.calc.GPI.xarray import potential_intensity# ---------------- 读数据 -----------------ds_s = xr.open_dataset("/home/mw/input/gpi4714/ERA5_single.nc")ds_m = xr.open_dataset("/home/mw/input/gpi4714/ERA5_multilevel.nc")# 如果下载的是逐小时，但只想做日平均，可 resample# ds_s = ds_s.resample(time='1D').mean()# ds_m = ds_m.resample(time='1D').mean()# 确保单位属性（ERA5 本来就有，保险起见再写一次）ds_s['sst'].attrs['units']  = 'K'ds_s['msl'].attrs['units']  = 'Pa'ds_m['t'].attrs['units']    = 'K'ds_m['q'].attrs['units']    = 'kg kg-1'# 混合比/比湿# 可选：只挑热带海域（SST≥26.5°C）mask = ds_s['sst'] >= 299.65# 26.5°Csst  = ds_s['sst'].where(mask)msl  = ds_s['msl'].where(mask)# 多层变量也做同样掩膜t   = ds_m['t'].where(mask)q   = ds_m['q'].where(mask)# ---------------- 计算 PI ---------------print("🌀 开始计算 Potential Intensity ...")pi_ds = potential_intensity(    sst             = sst,    psl             = msl,    pressure_levels = ds_m['pressure_level'],    temperature     = t,    mixing_ratio    = q)# ---------------- 保存 -----------------outfile = "ERA5_PI_aug2023.nc"pi_ds.to_netcdf(outfile)print(f"✅ 计算完成！结果已写入 {outfile}")

输出示例：

🌀 开始计算 Potential Intensity ...Missing values detected in: SST(NaN), PSL(NaN), temperature(NaN), mixing_ratio(NaN). Using missing value handling version.✅ 计算完成！结果已写入 ERA5_PI_aug2023.nc

计算结果可视化示例：

PI 平均值分布

快速可视化：

import matplotlib.pyplot as pltpi_ds = xr.open_dataset("ERA5_PI_aug2023.nc")pi_ds['pi'].mean('valid_time').plot(cmap='jet', vmin=20, vmax=90)plt.title("ERA5 Skyborn Potential Intensity (mean)")plt.tight_layout()plt.show()

查看最小气压：

pi_ds.min_pressure[16].plot()

最小气压分布示例：

最小气压分布

4. 细化绘图

import xarray as xrimport matplotlib.pyplot as pltimport matplotlib.ticker as mtickerimport cartopy.crs as ccrsimport cartopy.feature as cfeatureimport numpy as npimport matplotlib.patches as mpatchesfrom matplotlib.colors import ListedColormapimport osimport cmaps# ------------------------------------------------------# 1. 数据# ------------------------------------------------------ds   = xr.open_dataset('ERA5_PI_aug2023.nc')pi   = ds['pi']                                # (48, 121, 161)sel  = pi.valid_time.dt.hour.isin([0, 6, 12, 18])pi   = pi.sel(valid_time=sel)                  # 8 个时刻pmin = ds['min_pressure'].sel(valid_time=pi.valid_time)# ------------------------------------------------------# 2. 颜色表# ------------------------------------------------------cmap   = 'jet'levels = np.arange(20, 91, 5)norm   = plt.Normalize(20, 90)# ------------------------------------------------------# 3. 底图函数# ------------------------------------------------------defadd_geo(ax, left=False, bottom=False):    ax.set_extent([100, 140, 0, 30], crs=ccrs.PlateCarree())    ax.add_feature(cfeature.LAND.with_scale('50m'),                   facecolor=[.95, .95, .95], edgecolor='none', zorder=1)    ax.add_feature(cfeature.COASTLINE.with_scale('50m'),                   linewidth=.4, edgecolor='0.15', zorder=2)    gl = ax.gridlines(draw_labels=True, linewidth=.15,                      color='gray', alpha=.4, linestyle='-', zorder=3)    gl.top_labels   = gl.right_labels = False    gl.left_labels  = left    gl.bottom_labels= bottom    gl.xlocator = mticker.FixedLocator(np.arange(100, 141, 10))    gl.ylocator = mticker.FixedLocator(np.arange(0, 31, 10))    gl.xlabel_style = {'fontsize': 6}    gl.ylabel_style = {'fontsize': 6}# ------------------------------------------------------# 4. 画布# ------------------------------------------------------fig = plt.figure(figsize=(8,6), dpi=300)  # 180 mm 宽axes = []for i in range(8):    left  = (i % 4 == 0)          # 最左列    bottom= (i // 4 == 1)         # 最下行    ax = fig.add_subplot(2, 4, i+1, projection=ccrs.PlateCarree())    add_geo(ax, left=left, bottom=bottom)    axes.append(ax)# ------------------------------------------------------# 5. 绘图# ------------------------------------------------------for ax, t in zip(axes, pi.valid_time):    pi_t = pi.sel(valid_time=t)    pm_t = pmin.isel(valid_time=i)    cf = ax.contourf(pi_t.longitude, pi_t.latitude, pi_t,                     levels=levels, cmap=cmap, norm=norm,                     extend='both', transform=ccrs.PlateCarree(), zorder=1)    cs = ax.contour(pm_t.longitude, pm_t.latitude, pm_t,                    levels=6, colors='0.15', linewidths=.3,                    transform=ccrs.PlateCarree(), zorder=2)    ax.set_title(t.dt.strftime('%Y-%m-%d %HZ').values,                 size=7, pad=1)# ------------------------------------------------------# 6. colorbar# ------------------------------------------------------cbar_ax = fig.add_axes([0.25, 0.06, 0.5, 0.015])cbar = fig.colorbar(cf, cax=cbar_ax, orientation='horizontal',                    extend='both', extendfrac=0.02, ticks=levels[::2])cbar.ax.set_xlabel('Potential Intensity (kt)', fontsize=7, labelpad=3)cbar.ax.tick_params(labelsize=6, width=.3)# ------------------------------------------------------# 7. 保存# ------------------------------------------------------# plt.subplots_adjust(left=0.04, right=0.96, top=0.94, bottom=0.12,#                     hspace=0.02, wspace=0.06)plt.show()# plt.savefig('fig_pi_8panel.pdf', dpi=300, bbox_inches='tight')# plt.close()

细化绘图示例（8 面板图）：

细化绘图

5. 常见坑提醒

  1. 垂直层太少 → PI 偏低：至少 15 层，最好 30+ 层；上面下载的 37 层足够。
  2. 把 RH 直接当 q 用 → 结果全错，ERA5 的 q 就是混合比，无需换算。
  3. 单位属性缺失 → Skyborn 不会自动识别，脚本里已补写。
  4. 内存爆掉：全球 0.25°× 逐小时 ×37 层非常庞大，建议先裁剪区域/时段，或分块 dask。

这样，从"零"到拿到一张漂亮的台风潜在强度图， 全程只需要两条命令：

python get_era5_for_PI.py   # 下载python calc_PI_ERA5.py      # 计算

快去试试吧！祝科研顺利 🌪️

EarthAi

微信公号：Earth-Ai

商务联系：Earth_Ai