如何用xarray合并ERA-5数据逐年阈值首次跨越日的计算结果？

我使用包含1991-2020年逐日数据的ERA-5 xarray数据集，其中变量"ends"的0值表示未跨越阈值，正值表示已跨越阈值。我需要找出每年首次跨越该阈值的年积日（DoY），并将得到的30个DataArray合并为一个带有year维度（存储年份数值）的单一DataArray或Dataset，效果类似于执行dsxr.groupby('time.year').mean()的结果。

若处理完整数据集，可通过在time维度上调用idxmin()方法获取首次跨越的年积日，但我尝试循环处理单年数据时，无法成功合并结果。以下是我目前编写的代码：

for year in range(1991, 2021):
    yr = str(year)
    hwt_da = dsxr.ends.sel(time=slice(yr,yr))  # 获取原数据集对应年份的"ends"变量子集
    dT = hwt_da.indexes['time'].strftime('%j').astype(int)  # 将cftime格式转为整数类型的年积日
    hwt_da['time'] = dT   # 把年积日赋值给time维度
    hh = hwt_da.where(hwt_da>0).idxmin(dim='time').assign_coords({'year': year})  # 获取首次跨越阈值的年积日
    if year == '1991':
        da = hh.copy()  # 复制第一年的结果
    else:
        xr.concat([da, hh], dim='year')  # 拼接后续年份的结果

问题分析与修正 #

原代码有两个关键错误导致合并失败：

1. 条件判断if year == '1991'类型不匹配：year是整数类型，应与整数1991比较，而非字符串
2. xr.concat未重新赋值：执行拼接后结果未保存回da，导致da始终只保留第一年的数据

另外，获取年积日可以用xarray内置的dt.dayofyear方法，比手动格式化更简洁可靠。

修正后的循环代码 #

import xarray as xr

# 初始化空列表存储每年计算结果
results = []

for year in range(1991, 2021):
    # 按年份筛选数据
    hwt_da = dsxr.ends.sel(time=slice(str(year), str(year)))
    # 用内置方法获取年积日，替换time维度坐标
    hwt_da = hwt_da.assign_coords(time=hwt_da.time.dt.dayofyear)
    # 筛选跨越阈值的记录，取首次出现的年积日并添加年份坐标
    first_doy = hwt_da.where(hwt_da > 0).idxmin(dim='time').assign_coords(year=year)
    # 将当前年份结果加入列表
    results.append(first_doy)

# 按year维度拼接所有年份结果
final_da = xr.concat(results, dim='year')

更高效的groupby方案（推荐） #

无需手动循环，直接利用xarray的groupby功能一步完成计算与合并，逻辑更简洁且性能更优：

# 为数据集添加年积日和年份坐标
ds_with_doy = dsxr.ends.assign_coords(
    doy=dsxr.ends.time.dt.dayofyear,
    year=dsxr.ends.time.dt.year
)

# 按年份分组，筛选阈值>0的记录，取每组最小的doy（即首次跨越的年积日）
final_da = ds_with_doy.where(ds_with_doy > 0).groupby('year').min(dim='time')['doy']

或者更直观的自定义分组处理：

def get_first_doy(group):
    # 找到每组中第一个满足阈值条件的时间点，提取对应的年积日
    first_time_idx = group.where(group > 0).idxmin(dim='time')
    return group.doy.sel(time=first_time_idx)

# 按年份分组并应用自定义函数
final_da = dsxr.ends.groupby('time.year').apply(get_first_doy)

以上两种groupby方法都会直接生成带有year维度的DataArray，效果与groupby('time.year').mean()完全一致。

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