如何识别Metpy风寒计算结果的--标记值并转为xarray的NaN

使用metpy.calc.windchill计算风寒值时，超出温度与风速阈值的区域会显示为--，但尝试用掩码、.where()或filled()方法将这些值转为np.nan时失败：调用ma.getmask(windchill_ma)得到的掩码全为False，查看wind_chill.values发现标记--的位置实际存在数值，无法正确识别无效区域。

所用数据为谷歌云ERA5数据集的t2m、uwind、vwind，均为包含lat、lon、time维度的3D数组，代码如下：

t2m = xr.open_dataarray('ERA5_t2m_hourly_2004_2024.nc')
vwind = xr.open_dataarray('ERA5_vwind_hourly_2004_2024.nc')
uwind = xr.open_dataarray('ERA5_uwind_hourly_2004_2024.nc')

def wind_tot(uwind, vwind):
    wind_mag = np.sqrt(vwind**2 + uwind**2)
    wind_dir = np.arctan2(vwind/wind_mag, uwind/wind_mag)
    wind_dir = wind_dir * 180/np.pi
    
    return wind_mag, wind_dir

wind_mag, wind_dir = wind_tot(uwind, vwind)

wind_chill = metpy.calc.windchill(
                                 t2m * units.K,
                                 wind_mag * units('m/s'),
                                    )
windchill_ma = wind_chill.to_masked_array()
mask = ma.getmask(windchill_ma)

解决方案 #

MetPy的windchill函数不会自动生成掩码，而是对超出阈值的情况返回原气温值，所以无法通过返回结果的掩码识别无效区域。需要基于风寒值的有效计算条件手动创建掩码：

1. 明确有效阈值
   MetPy风寒值的有效计算条件为：

• 气温 ≤ 10℃（即283.15K）
• 风速 ≥ 4.8 km/h（换算为国际单位是≈1.333 m/s）

2. 创建无效区域掩码
   用原始的t2m和wind_mag生成布尔掩码，标记不符合条件的区域：

# 直接用m/s阈值判断
valid_mask = (t2m <= 283.15) & (wind_mag >= 1.333)

3. 替换无效值为NaN并保留xarray结构
   用xarray的.where()方法，将无效区域的值替换为np.nan：

# 先去除单位方便处理，若需保留可后续重新添加
wind_chill_no_units = wind_chill.metpy.dequantify()
# 应用掩码替换无效值
wind_chill_clean = wind_chill_no_units.where(valid_mask, np.nan)

完整修改后的代码 #

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

t2m = xr.open_dataarray('ERA5_t2m_hourly_2004_2024.nc')
vwind = xr.open_dataarray('ERA5_vwind_hourly_2004_2024.nc')
uwind = xr.open_dataarray('ERA5_uwind_hourly_2004_2024.nc')

def wind_tot(uwind, vwind):
    wind_mag = np.sqrt(vwind**2 + uwind**2)
    wind_dir = np.arctan2(vwind/wind_mag, uwind/wind_mag)
    wind_dir = wind_dir * 180/np.pi
    
    return wind_mag, wind_dir

wind_mag, wind_dir = wind_tot(uwind, vwind)

# 计算风寒值
wind_chill = mpcalc.windchill(t2m * units.K, wind_mag * units('m/s'))

# 创建有效区域掩码
valid_mask = (t2m <= 283.15) & (wind_mag >= 1.333)  # 1.333m/s≈4.8km/h

# 去除单位并替换无效值为NaN，保留xarray结构
wind_chill_clean = wind_chill.metpy.dequantify().where(valid_mask, np.nan)

# 若需保留单位，可执行以下步骤
# wind_chill_clean = wind_chill_clean * units.K

验证方式 #

查看无效值转换情况：

# 统计NaN的数量
print(wind_chill_clean.isnull().sum())

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