基于指定公式优化南半球MCS指数计算代码的技术问询

目标公式 #

SA−MCS指数 = [(0–6 km风切变 –20.01)÷7.87] + [(775 hPa水平梯度 –4.84×10⁻⁵)÷(5.65×10⁻⁵)] + {–[(800 hPa垂直速度ω + 0.27)÷0.29]} + {–[(四层抬升指数LI + 2.17)÷2.23]}

修改后的完整代码 #

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import metpy.calc as mpcalc
from metpy.units import units

# 打开GRIB格式的数据文件
# 等压面层数据
ds_isobaric = xr.open_dataset(
    'cdas1_2022011700_.t00z.pgrbh00.grib2', 
    engine='cfgrib', 
    backend_kwargs={'filter_by_keys': {'typeOfLevel': 'isobaricInhPa'}}
).metpy.parse_cf()

# 地面层派生数据（抬升指数、可降水量）
ds_surface = xr.open_dataset(
    'cdas1_2022011700_.t00z.pgrbh00.grib2', 
    engine='cfgrib', 
    backend_kwargs={'filter_by_keys': {'typeOfLevel': 'pressureFromGroundLayer'}}
).metpy.parse_cf()

# ---------------------- 提取所需气象要素 ----------------------
# 0-6km风切变对应层：1000hPa（近地面）和400hPa（约6km高度）
u_1000 = ds_isobaric.u.sel(isobaricInhPa=1000)
v_1000 = ds_isobaric.v.sel(isobaricInhPa=1000)
u_400 = ds_isobaric.u.sel(isobaricInhPa=400)
v_400 = ds_isobaric.v.sel(isobaricInhPa=400)

# 775hPa要素（水平温度梯度相关）
u_775 = ds_isobaric.u.sel(isobaricInhPa=775)
v_775 = ds_isobaric.v.sel(isobaricInhPa=775)
t_775 = ds_isobaric.t.sel(isobaricInhPa=775)

# 800hPa垂直速度ω（注意：气象中ω向下为正，公式中已适配符号）
w_800 = ds_isobaric.w.sel(isobaricInhPa=800)

# 四层抬升指数LI和可降水量
pli = ds_surface.pli
pwat = ds_surface.pwat

# ---------------------- 计算SA-MCS指数核心变量 ----------------------
# 计算0-6km风切变：风矢量差的模（正确的风切变计算方式，而非风速差）
shear_vector_u = u_400 - u_1000
shear_vector_v = v_400 - v_1000
shear_0_6km = mpcalc.wind_speed(shear_vector_u, shear_vector_v)

# 计算775hPa水平温度平流（对应公式中的"水平梯度"项）
# 用numpy.gradient计算温度的空间梯度
dT_dx, dT_dy = np.gradient(t_775, axis=(0, 1))
# 温度平流公式：-V·∇T
gradiente_horizontal = - (u_775 * dT_dx + v_775 * dT_dy)

# 按照目标公式计算SA-MCS指数
sa_mcs = (
    (shear_0_6km - 20.01) / 7.87 
    + (gradiente_horizontal - 4.84e-5) / 5.65e-5 
    + (- (w_800 + 0.27) / 0.29) 
    + (- (pli + 2.17) / 2.23)
)

# ---------------------- 绘图部分 ----------------------
fig = plt.figure(figsize=(16, 14))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# 添加地理要素
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS)
ax.add_feature(cfeature.STATES)

# 添加经纬网格
ax.gridlines(
    draw_labels=True, 
    linewidth=0.5, 
    color='gray', 
    alpha=0.5, 
    linestyle='--'
)

# 设置南半球区域范围
extent = [-75, -45, -40, -16]
ax.set_extent(extent, crs=ccrs.PlateCarree())
ax.set_title("南半球SA-MCS指数", fontsize=16, fontweight='bold')

# 应用掩码：过滤可降水量<27、LI>0及指数在(-0.5,0.5)的区域
masked_sa_mcs = np.where(
    (pwat < 27) & (pli > 0) & (sa_mcs >= -0.5) & (sa_mcs <= 0.5), 
    np.nan, 
    sa_mcs
)

# 绘制填色图
clevs = np.arange(-1.5, 5)
contourf = ax.contourf(
    ds_isobaric['longitude'], 
    ds_isobaric['latitude'], 
    masked_sa_mcs, 
    levels=clevs,
    cmap='rainbow', 
    transform=ccrs.PlateCarree()
)

# 添加色标
cbar = plt.colorbar(
    contourf, 
    ax=ax, 
    orientation='horizontal', 
    aspect=40, 
    pad=0.02
)
cbar.set_label('SA-MCS指数', size=14)

plt.show()

关键修改点说明 #

• 风切变计算修正：原代码用上下层风速差，改为风矢量差的模，更符合0-6km风切变的物理定义
• 公式系数完全匹配：替换原代码中所有旧系数，严格对齐目标SA-MCS指数公式
• 变量名优化：将模糊的ds/ds1改为ds_isobaric/ds_surface，提升代码可读性
• 注释增强：补充关键计算步骤的气象意义说明，便于后续维护
• 输出内容明确：将绘图标题更新为"南半球SA-MCS指数"，清晰标识输出结果

内容的提问来源于stack exchange，提问作者leojim19