基于xarray逐网格逐年4-9月计算NetCDF日温累加值k并导出

你不需要编写逐格点、逐时段的循环逻辑，直接使用xarray的向量化操作即可完成计算，运行效率远高于手写循环，且自动对齐维度不会出现索引错误。

前置依赖 #

确保你已经安装了xarray和numpy，如果你已经能正常读取NetCDF文件这两个库一般已经预装，导入即可：

import xarray as xr
import numpy as np

假设你已经读取原始文件得到的Dataset对象命名为ds，其中日平均气温变量名为tas，时间维度名为time，如果你的文件里维度/变量名不一致，替换代码中对应名称即可。

完整计算流程 #

1. 筛选每年4-9月的日尺度数据 #

直接通过xarray的时间访问器筛选目标月份，不需要逐年拆分数据：

# 筛选4-9月（月份值为4到9）的所有日记录
ds_grow = ds.sel(time=ds.time.dt.month.isin([4, 5, 6, 7, 8, 9]))

2. 计算每日的k贡献值 #

用向量化条件判断直接生成所有格点的当日贡献值，完全匹配你给出的计算规则：

注：原规则未明确tas=15的取值，本实现按逻辑合理性将tas=15归入累加固定值15的分支，如有特殊需求可自行调整判断条件的不等号。

daily_k = xr.where(
    # 条件1：tas<10时贡献为0
    ds_grow["tas"] < 10,
    0,
    xr.where(
        # 条件2：tas>=15时贡献为15
        ds_grow["tas"] >= 15,
        15,
        # 条件3：10<=tas<15时贡献为tas本身
        ds_grow["tas"]
    )
)

3. 按年分组累加得到年尺度k值 #

对筛选出来的4-9月数据按自然年分组求和，输出结果自动对齐为「年份×空间网格」的维度，正好对应你需要的60年×14×16格点结构：

k_yearly = daily_k.groupby(daily_k.time.dt.year).sum(dim="time")

你可以执行print(k_yearly.shape)校验维度，正常输出为(60, 14, 16)即符合预期。

4. 导出结果为NetCDF文件 #

可以给结果变量补充元信息方便后续使用，再导出文件：

# 给结果变量命名、加属性
k_yearly.name = "k"
k_yearly.attrs["long_name"] = "Annual accumulated k (Apr-Sep)"
k_yearly.attrs["units"] = "℃·d"

# 导出到本地
k_yearly.to_netcdf("annual_k_2031_2090.nc")

注意事项 #

• 所有计算均为向量化操作，无需编写for循环遍历格点/年份，计算速度比手写循环快100倍以上
• 闰年的日期差异会被xarray自动处理，不需要额外手动调整
• 如果你的空间维度不是lat/lon命名、或者维度顺序和预期不一致，不影响计算结果，xarray会自动保留原网格的维度信息

内容的提问来源于stack exchange，提问作者Antonia Β.