使用Xarray(NetCDF)循环计算后，如何将结果追加至现有数据集

刚接触Python编程，正在使用Xarray处理NetCDF格式的气候模型数据，目标是对比未来投影数据（ssp子集）与历史数据（hr4子集），计算积雪（'snw'）的差值以展示变化趋势。目前已编写三重循环遍历时间、纬度、经度维度计算差值，代码如下：

def get_diff_snw(hist_avg, proj_value):
    return proj_value - hist_avg

for i, time in enumerate(ds_ssp_subset['snw']):
    # For each time, iterate through the latitudes    
    for j, lat in enumerate(time):
        # For each latitude, iterate through the longitudes
        for k, lon in enumerate(lat):
            print(
                get_diff_snw(
                    ds_hr4_running_subset['snw'][i, j, k].values,
                    lon.values,
                )
            )

注：ssp代表未来投影数据集，hr4代表历史数据集。

当前循环已输出一系列积雪差值，但不知如何将这些带对应坐标关联的结果追加至历史数据集中，恳请帮助。

1. 用Xarray矢量化操作替代三重循环 #

Xarray的核心优势之一是支持矢量化计算，无需手动遍历维度——它会自动匹配坐标并完成广播计算，效率远高于三重循环，还能直接保留所有坐标信息。

直接计算积雪差值：

# 直接计算未来投影与历史数据的积雪差值，自动保留时间、纬度、经度坐标
diff_snw = ds_ssp_subset['snw'] - ds_hr4_running_subset['snw']
# 给差值变量添加属性说明（可选，提升数据可读性）
diff_snw.attrs = {
    'long_name': 'Snow Depth Difference (Future - Historical)',
    'units': ds_hr4_running_subset['snw'].attrs.get('units', '')
}

2. 将差值追加到历史数据集 #

使用Xarray的assign()方法，把新的差值变量直接添加到历史数据集中：

# 把差值变量追加到历史数据集，生成新的数据集（原数据集不会被修改）
ds_hr4_updated = ds_hr4_running_subset.assign(snw_diff=diff_snw)

3. 验证结果 #

检查新数据集是否包含差值变量及对应坐标：

# 查看数据集变量列表
print(ds_hr4_updated.data_vars)
# 查看差值变量的坐标信息
print(ds_hr4_updated['snw_diff'].coords)

完整改进代码示例 #

# 计算积雪差值（矢量化操作）
diff_snw = ds_ssp_subset['snw'] - ds_hr4_running_subset['snw']
diff_snw.attrs = {
    'long_name': 'Snow Depth Difference (SSP Projection - Historical HR4)',
    'units': ds_hr4_running_subset['snw'].attrs.get('units', '')
}

# 追加到历史数据集
ds_hr4_updated = ds_hr4_running_subset.assign(snw_diff=diff_snw)

# 可选：保存更新后的数据集到NetCDF文件
ds_hr4_updated.to_netcdf('hr4_with_snw_diff.nc')

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