如何处理预报步骤的时间重叠？将观测数据维度转为(time,steps)

我已下载2021年ERA5再分析（reanalysis-era5-pressure-levels，变量observation）和ECMWF季节预报（seasonal-original-pressure-levels，变量forecast）的u、v风分量，计划用properscoring库的**Continuous Ranked Probability Score（CRPS）**评估预报性能。

评估思路是通过xr.apply_ufunc基于number（集合成员）维度，为每个时间、纬度、经度计算CRPS，代码示例如下：

import properscoring as ps

crps = xr.apply_ufunc(
    ps.crps_ensemble,
    observation,
    forecast,
    input_core_dims=[[], ["number"]],
    output_core_dims=[[]],
    vectorize=True,
)

但observation的维度为("time", "latitude", "longitude")，而forecast的维度为("time", "latitude", "longitude", "steps", "number")，二者无法兼容（xr.align(..., join = "exact")因time和steps维度不匹配失败）。需要将observation的time维度转换为("time", "steps")以适配预报数据。

已尝试方案 #

使用backend_kwargs #

在另一项目中，我曾用xr.open_dataset(..., backend_kwargs={"time_dims": ("valid_time",)})让valid_time匹配time维度，但当前项目中valid_time与time + steps存在重叠（比如time=2021-01-01 + steps=36:00:00 > time=2021-01-02），会丢失数据，示例代码：

forecast = xr.open_dataset(
    "data/raw/forecast.grib",
    engine="cfgrib",
    backend_kwargs=dict(time_dims=("valid_time",)),
)
# forecast丢失steps > 24:00:00对应的数据

拼接为新维度 #

尝试用滚动窗口拼接生成不同steps：为每个日期提取后续36小时作为steps并创建新维度，但既没生成目标维度（可能xr.concat不是正确函数），速度还极慢（两个月数据耗时超6分钟），示例代码：

xr.concat(
    [
        observation.sel(time = slice(date, date + pd.Timedelta("36H")))
        for date in observation.time.values
    ],
    dim = "steps"
)

正确解决方法 #

核心思路是为观测数据生成与预报数据匹配的(time, steps)维度，即对每个预报初始时间time，找到对应time + steps时刻的观测值，通过广播和索引实现高效匹配。

步骤1：计算预报的有效时间（valid_time） #

首先从预报数据中生成每个初始时间+步长对应的有效时间，方便后续匹配观测：

# 生成预报的valid_time维度：time + steps
forecast["valid_time"] = forecast.time + forecast.steps

步骤2：将观测数据重索引到预报的valid_time网格 #

利用xarray的索引和广播机制，将观测数据扩展为与预报匹配的维度：

# 将observation的time重命名为valid_time，统一匹配字段
obs_valid = observation.rename({"time": "valid_time"})

# 匹配每个预报(time, steps)对应的观测值，自动扩展为(time, latitude, longitude, steps)维度
obs_matched = obs_valid.sel(valid_time=forecast.valid_time)

步骤3：调整维度并计算CRPS #

确保观测和预报的维度逻辑对应，调用xr.apply_ufunc计算CRPS：

crps = xr.apply_ufunc(
    ps.crps_ensemble,
    obs_matched,  # 维度：(time, latitude, longitude, steps)
    forecast[["u", "v"]],  # 维度：(time, latitude, longitude, steps, number)
    input_core_dims=[[], ["number"]],
    output_core_dims=[[]],
    vectorize=True,
    dask="parallelized",  # 大数据场景下启用dask并行加速
    output_dtypes=[float],
)

# 最终crps的维度为(time, latitude, longitude, steps)

优化说明 #

• 无循环高效匹配：利用xarray的广播索引机制，避免显式循环，速度远快于滚动窗口拼接
• 完整保留数据：通过valid_time精准匹配，不会因跨天步长丢失数据
• 支持并行计算：结合dask加载数据，可进一步提升大数据量下的计算效率

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