使用scipy.stats计算NetCDF趋势时触发_UFuncBinaryResolutionError报错

scipy.stats.linregress 要求输入的自变量序列为数值型可计算格式，不支持直接接收datetime64[ns]类型的时间坐标。函数内部执行线性回归计算时需要对自变量做加减、求差等算术操作，遇到datetime64类型的输入就会触发dtype不匹配的ufunc add报错。

• 第一步：将原生datetime格式的时间坐标转换为浮点型数值刻度，禁止直接把xarray读取的原始time坐标传入回归函数。根据计算需求选择对应的转换方式：
  • 仅需要逐像素趋势的相对大小、显著性p值：直接生成长度和时间维度一致的整数索引作为自变量，time_numeric = np.arange(ntime)，计算得到的斜率单位为「每1个时间步长的GDI变化量」
  • 需要斜率具备实际物理单位（如每小时/每天的GDI变化幅度）：将时间戳转换为相对于序列起始时刻的间隔值，示例为转换为小时刻度：

    time_start = time.values[0]
    # 计算每个时间点距离起始时刻的小时数，输出为可计算的浮点类型
    time_numeric = (time.values - time_start) / np.timedelta64(1, 'h')
    

• 第二步：修正原代码的索引错位问题。原代码初始化的结果数组shape为(nlat, nlon)（来自GDI[0,:,:]的维度顺序），但赋值时使用了slope0[jlon, ilat]的经纬度倒序索引，会导致最终计算的结果和经纬度坐标对应错位，需要将索引调整为和数组维度顺序匹配的[ilat, jlon]。
• 第三步：替换回归函数的入参，将原始time参数替换为转换后的time_numeric序列，逐调用的返回值按对应像素位置存入结果数组。

import numpy as np
import scipy.stats

nlat = len(lat)
nlon = len(lon)
ntime = len(time)

# 初始化结果数组，维度顺序与GDI的经纬度维度保持一致
slope0      =  np.zeros((nlat, nlon))
intercept0  =  np.zeros((nlat, nlon))
r_value0    =  np.zeros((nlat, nlon))
p_value0    =  np.zeros((nlat, nlon))
std_err0    =  np.zeros((nlat, nlon))

# 预处理时间为可计算的数值格式
time_start = time.values[0]
time_numeric = (time.values - time_start) / np.timedelta64(1, 'h')

for ilat in range(nlat):
    for jlon in range(nlon): 
        reg_res = scipy.stats.linregress(time_numeric, GDI[:, ilat, jlon])
        slope0[ilat, jlon] = reg_res.slope
        intercept0[ilat, jlon] = reg_res.intercept
        r_value0[ilat, jlon] = reg_res.rvalue
        p_value0[ilat, jlon] = reg_res.pvalue
        std_err0[ilat, jlon] = reg_res.stderr

效率优化提示：逐像素双重循环计算速度较慢时，可直接使用xarray内置的polyfit方法批量完成全数组的线性趋势计算，无需手动编写循环：

# 沿time维度做1次多项式拟合，自动处理高维数组的逐像素计算
fit_coeff = GDI.polyfit(dim='time', deg=1)
# 提取得到的一阶项（即斜率）
slope_result = fit_coeff.sel(degree=1)['polyfit_coefficients']

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