使用MetPy计算热力指数时出现IndexError的原因排查

问题场景 #

尝试使用MetPy计算LCL、CCL等热力指数，输出的探空观测数据如下：

气压(p): [34.5 30.0 29.1 25.6 20.0 19.9 16.4 15.1 12.0 10.3 10.1] hectopascal
气温(t): [-56.3 -56.1 -56.9 -46.9 -50.7 -50.9 -45.3 -46.7 -39.3 -44.3 -44.3] degree_Celsius
露点(td): [-88.3 -88.1 -88.9 -84.9 -85.7 -85.9 -84.3 -84.7 -81.3 -84.3 -84.3] degree_Celsius

运行代码时在调用mpcalc.ccl(p, t, td)处触发IndexError: index -1 is out of bounds for axis 0 with size 0，完整代码及报错追踪如下：

运行代码 #

import numpy as np
import metpy.calc as mpcalc
from metpy.units import units

print(p, t, td)

# LCL, CCL, LFC, EL
lcl_p, lcl_t = mpcalc.lcl(p[0], t[0], td[0])
ccl_p, ccl_t, t_c = mpcalc.ccl(p, t, td)
lfc_p, lfc_t = mpcalc.lfc(p, t, td, which='bottom')
el_p, el_t = mpcalc.el(p, t, td)


# CAPE, CIN
prof = mpcalc.parcel_profile(p, t[0], td[0]).to('degC')
cape, cin = mpcalc.cape_cin(p, t, td, prof)

print(f'LCL = {lcl_p.magnitude:.2f} {lcl_p.units}')
print(f'CCL = {ccl_p.magnitude:.2f} {ccl_p.units}')
print(f'LFC = {lfc_p.magnitude:.2f} {lfc_p.units}')
print(f'EL = {el_p.magnitude:.2f} {el_p.units}')
print(f'CAPE = {cape.magnitude:.2f} {cape.units}')
print(f'CIN = {cin.magnitude:.2f} {cin.units}')

报错追踪 #

[10.1 10.3 12.0 15.1 16.4 19.9 20.0 25.6 29.1 30.0 34.5] hectopascal [-44.3 -44.3 -39.3 -46.7 -45.3 -50.9 -50.7 -46.9 -56.9 -56.1 -56.3] degree_Celsius [-84.3 -84.3 -81.3 -84.7 -84.3 -85.9 -85.7 -84.9 -88.9 -88.1 -88.3] degree_Celsius
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[70], line 14
     12 # LCL, CCL, LFC, EL
     13 lcl_p, lcl_t = mpcalc.lcl(p[0], t[0], td[0])
---> 14 ccl_p, ccl_t, t_c = mpcalc.ccl(p, t, td)
     15 lfc_p, lfc_t = mpcalc.lfc(p, t, td, which='bottom')
     16 el_p, el_t = mpcalc.el(p, t, td)

File ~\AppData\Roaming\Python\Python310\site-packages\metpy\xarray.py:1328, in preprocess_and_wrap.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
   1325     _mutate_arguments(bound_args, units.Quantity, lambda arg, _: arg.m)
   1327 # Evaluate inner calculation
-> 1328 result = func(*bound_args.args, **bound_args.kwargs)
   1330 # Wrap output based on match and match_unit
   1331 if match is None:

File ~\AppData\Roaming\Python\Python310\site-packages\metpy\units.py:325, in check_units.<locals>.dec.<locals>.wrapper(*args, **kwargs)
    322 @functools.wraps(func)
    323 def wrapper(*args, **kwargs):
    324     _check_units_inner_helper(func, sig, defaults, dims, *args, **kwargs)
--> 325     return func(*args, **kwargs)

File ~\AppData\Roaming\Python\Python310\site-packages\metpy\calc\thermo.py:574, in ccl(pressure, temperature, dewpoint, height, mixed_layer_depth, which)
    572 # In the case of multiple CCLs, select which to return
    573 if which == 'top':
--> 574     x, y = x[-1], y[-1]
    575 elif which == 'bottom':
    576     x, y = x[0], y[0]

File ~\AppData\Roaming\Python\Python310\site-packages\pint\facets\numpy\quantity.py:238, in NumpyQuantity.__getitem__(self, key)
    236 def __getitem__(self, key):
    237     try:
--> 238         return type(self)(self._magnitude[key], self._units)
    239     except PintTypeError:
    240         raise

IndexError: index -1 is out of bounds for axis 0 with size 0

问题根源 #

1. 探空数据顺序不符合要求 #

MetPy的热力计算函数（包括ccl）强制要求气压数据按从高到低（地面到高空，即气压值递减）的顺序排列，但你的气压数组是[10.1 10.3 12.0 ... 34.5]，是气压递增的顺序（从高空到地面），这会导致函数无法正确识别大气垂直结构，直接干扰CCL的计算逻辑。

2. CCL不存在触发空数组索引 #

CCL的定义是地面气块干绝热上升后，与环境露点曲线的交点高度。当数据顺序错误，或者当前探空条件下环境露点始终低于干绝热抬升的气块温度时，ccl函数会生成空的结果数组。后续代码尝试访问空数组的[-1]索引（默认取最顶部的CCL），就会触发IndexError。

解决建议 #

• 修正数据顺序：将气压、气温、露点数组按气压递减的顺序反转：

  p = p[::-1]
  t = t[::-1]
  td = td[::-1]
  

• 添加异常处理：针对CCL不存在的情况捕获异常，避免程序崩溃：

  try:
      ccl_p, ccl_t, t_c = mpcalc.ccl(p, t, td)
  except IndexError:
      print("无法计算CCL：当前探空条件下不存在对流凝结高度")
      ccl_p, ccl_t, t_c = None, None, None
  

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