如何计算网格化数据集与时间序列的Spearman秩相关系数？

针对你遇到的网格化xarray数据与单时间序列的逐格点Spearman相关计算问题，以下是两种可靠的实现方式：

方法一：用apply_ufunc结合scipy.stats.spearmanr（推荐） #

这种方式能精确控制逐格点的计算逻辑，同时可一并输出相关系数和p值，避免维度匹配问题：

1. 导入依赖库

import xarray as xr
from scipy.stats import spearmanr

2. 定义计算单个格点相关的包装函数

def calc_spearman(x, y):
    # x：单个格点的时间序列（一维）；y：全局时间序列（一维）
    corr_coef, p_value = spearmanr(x, y)
    # 返回包含相关系数和p值的数组，新增stat维度区分两者
    return xr.DataArray([corr_coef, p_value], dims=["stat"])

3. 应用到整个网格化数据集

# 对每个(lat, lon)格点，取对应time序列与nino_djf的time序列计算相关
result = xr.apply_ufunc(
    calc_spearman,
    spi_djf_cor,          # 网格化数据：(time, lat, lon)
    nino_djf,             # 单时间序列：(time,)
    input_core_dims=[["time"], ["time"]],  # 指定每个输入的核心计算维度为time
    output_core_dims=[["stat"]],           # 指定输出新增stat维度
    vectorize=True,                        # 开启向量化，逐格点计算
    dask="allowed"                         # 若使用dask分块数据可保留此参数
)

# 拆分结果为相关系数和p值
spearman_corr = result.sel(stat=0)
spearman_pval = result.sel(stat=1)

方法二：扩展单时间序列维度后使用内置函数 #

如果坚持使用xs.spearman_r，需要先将单时间序列扩展为与网格化数据一致的维度（新增lat和lon），确保维度匹配后再计算：

# 将nino_djf扩展lat和lon维度，与spi_djf_cor维度对齐
nino_expanded = nino_djf.expand_dims(lat=spi_djf_cor.lat, lon=spi_djf_cor.lon)

# 逐time维度计算相关
spearman_corr = xs.spearman_r(spi_djf_cor, nino_expanded, dim="time")

原代码可能出错的原因 #

你之前的代码未处理维度匹配问题：网格化数据是三维（time, lat, lon），而单时间序列是一维（time），xs.spearman_r可能无法自动识别需要对每个(lat, lon)格点单独计算，导致结果不符合预期。

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