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

# 建议先建个干净环境
conda create -n era5_PI python=3.11 -y
conda activate era5_PI

# 一次性装好
pip install cdsapi skyborn xarray netcdf4 dask

url: https://cds.climate.copernicus.eu/api/v2
key: 你的uid:你的key

#!/usr/bin/env python
"""
一键下载 ERA5 计算台风 PI 所需变量
单层：sst, msl
多层：t, q  (37 层，1000 → 1 hPa，Skyborn 官方推荐)
运行前请先在 ~/.cdsapirc 里放好 CDS API key
"""
import os, cdsapi, datetime as dt
from pathlib import Path

# ---------------- 用户区域 ----------------
dates     = ["2023-08-01/to/2023-08-31"]   # 可写多段
area      = [30, 110, 15, 140]            # N/W/S/E
outdir    = 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

#!/usr/bin/env python
"""
用 Skyborn 计算 ERA5 的台风最大潜在强度 PI
支持 1D/3D/4D 场，自动识别单位 & 缺测
"""
import xarray as xr
from 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°C
sst  = 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

import matplotlib.pyplot as plt
pi_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()

import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
import matplotlib.patches as mpatches
from matplotlib.colors import ListedColormap
import os
import 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. 底图函数
# ------------------------------------------------------
def add_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()

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