使用xarray处理带时间偏移的数据集求和时遇数值偏差问题

问题背景 #

两个NetCDF数据集，包含latitude、longitude、time三维维度，空间维度完全一致，时间维度存在1天偏移。目标是将两个数据集的tagged_precip变量求和生成第三个合并数据集，但操作中遇到以下问题：

操作1：单数据集总和计算 #

执行代码：

nd4 = xarray.open_mfdataset("../chec_prep/TAG_PREP/04/*.nc")
nd5 = xarray.open_mfdataset("../chec_prep/TAG_PREP/05/*.nc")
print(nd4['tagged_precip'].sum().values)
print(nd5['tagged_precip'].sum().values)

得到结果：2682132500.0 和 3773539800.0，理论总和为6455672300.0。

操作2：直接相加的异常 #

执行nd6 = nd4 + nd5后验证总和：

print(nd6['tagged_precip'].sum().values)
print(nd5['tagged_precip'].sum().values + nd4['tagged_precip'].sum().values)

结果为2710198300.0 vs 6455672300.0，偏差极大。

操作3：分组求和的局部偏差 #

批量打开所有文件并按时间分组求和：

np = xarray.open_mfdataset("../chec_prep/TAG_PREP/*/*.nc", combine='nested', concat_dim='time').groupby('time').sum(dim='time')
np['tagged_precip'].sum().values == nd4['tagged_precip'].sum().values + nd5['tagged_precip'].sum().values

总体验证返回True，但单个时间步存在-256.0的偏差：

print('nd4: ', nd4['tagged_precip'].sum(dim=['latitude','longitude']).values[10])
print('nd5: ', nd5['tagged_precip'].sum(dim=['latitude','longitude']).values[9])
print('np: ', np['tagged_precip'].sum(dim=['latitude','longitude']).values[10])
print('nd4 + nd5: ', (nd4['tagged_precip'].sum(dim=['latitude','longitude']).values[10] + nd5['tagged_precip'].sum(dim=['latitude','longitude']).values[9]))
print('nd4 + nd5 - np: ', nd4['tagged_precip'].sum(dim=['latitude','longitude']).values[10] + nd5['tagged_precip'].sum(dim=['latitude','longitude']).values[9] - np['tagged_precip'].sum(dim=['latitude','longitude']).values[10])

问题原因分析 #

1. 直接相加的偏差根源：
   xarray的算术运算默认按坐标对齐，仅保留两个数据集共有的时间步，非重叠的时间步会被直接丢弃，因此总和仅包含交集部分的数值，与全量求和的预期不符。

2. 单个时间步的-256偏差原因：

• 数据类型限制：tagged_precip大概率使用16位整数（int16）存储，其取值范围为-32768~32767，求和时若数值超过范围会发生溢出，-256是2的整数次幂，符合整数溢出的特征；若使用float32存储，有限的精度也会导致累加时出现舍入误差。
• 时间戳精度隐形差异：若两个数据集的时间戳精度不一致（如一个精确到天，另一个精确到小时/分钟），分组求和时可能出现隐形的分组错误，但总体验证通过，因此数据类型问题的概率更高。

解决方法 #

1. 正确合并所有时间步 #

使用xarray.concat拼接两个数据集，保留所有时间步后再按需分组求和：

import xarray as xr

# 拼接两个数据集，保留全部时间步
combined = xr.concat([nd4, nd5], dim='time')
# 按时间分组求和（合并同一时间戳的数值）
combined_sum = combined.groupby('time').sum(dim='time')
# 验证总和是否符合预期
print(combined_sum['tagged_precip'].sum().values == nd4['tagged_precip'].sum().values + nd5['tagged_precip'].sum().values)

2. 消除数据类型偏差 #

打开数据集时将tagged_precip转换为更高精度的类型，避免溢出或精度损失：

nd4 = xr.open_mfdataset("../chec_prep/TAG_PREP/04/*.nc", dtype={'tagged_precip': 'float64'})
nd5 = xr.open_mfdataset("../chec_prep/TAG_PREP/05/*.nc", dtype={'tagged_precip': 'float64'})

3. 对齐时间戳精度 #

确保两个数据集的时间戳精度一致，比如统一截断到天级别：

nd4['time'] = nd4['time'].dt.floor('D')
nd5['time'] = nd5['time'].dt.floor('D')

4. 缺失值填充后相加 #

若需按时间对齐并保留所有时间步，可合并时填充缺失值为0后相加：

merged = xr.merge([nd4['tagged_precip'], nd5['tagged_precip']], fill_value=0)
merged_precip = merged['tagged_precip_x'] + merged['tagged_precip_y']
# 验证总和
print(merged_precip.sum().values == nd4['tagged_precip'].sum().values + nd5['tagged_precip'].sum().values)

内容的提问来源于stack exchange，提问作者Luís Brito