Python重采样降水数据与固定指数的Bootstrap相关分析报错排查

问题 #

我想用Bootstrap重采样方法检验降水数据和某指数的相关性显著性，需要把降水数据（cru_djf_uv）重采样/打乱1000次，指数数据（nino_djf）保持固定，然后计算每个打乱后的降水数据集和固定指数的相关性，生成经纬度区域的1000个随机相关值。

我写了如下代码：

# Number of bootstrap samples
num_bootstrap_samples = 1000

# Function to calculate spearman correlation
def calculate_correlation(data1, data2):
    return xs.spearman_r(data1, data2, dim='time')

# Calculate the observed correlation
observed_correlation = calculate_correlation(cru_djf_uv, nino_djf)

# Initialize an array to store bootstrap samples
bootstrap_correlations = np.zeros((num_bootstrap_samples, cru_djf_uv.shape[1], cru_djf_uv.shape[2]))

# Perform bootstrap resampling
for i in range(num_bootstrap_samples):
    # Resample with replacement
    resampled_indices = np.random.choice(cru_djf_uv.shape[0], cru_djf_uv.shape[0], replace=True)
    resampled_cru = cru_djf_uv.isel(time=resampled_indices)
    
    # Calculate correlation for the resampled data
    bootstrap_correlations[i, :, :] = calculate_correlation(resampled_cru, nino_djf)

运行时出现错误：ValueError: indexes along dimension 'time' are not equal。我知道这是因为重采样后降水数据的时间索引和指数的不匹配，但这正是我想要的打乱效果，请问该怎么解决？

数据集详情 #

降水数据（cru_djf_uv）

<xarray.DataArray 'pre' (time: 30, lat: 18, lon: 6)>
array([[[243.      , 294.80002 , 309.40002 , 308.5     , 358.4     ,
         414.7     ],
        [240.1     , 254.      , 305.3     , 328.7     , 392.90002 ,
         401.9     ],
        [224.70001 , 271.7     , 282.1     , 393.4     , 487.9     ,
         478.5     ],
        ...,
        [181.09999 , 206.30002 , 234.20001 , 245.5     , 223.1     ,
         278.9     ],
        [176.8     , 207.      , 213.3     , 208.1     , 196.2     ,
         229.30002 ],
        [195.6     , 216.      , 211.7     , 186.6     , 172.4     ,
         165.3     ]],

       [[322.80002 , 370.9     , 427.7     , 393.7     , 389.09998 ,
         425.5     ],
        [274.9     , 292.7     , 425.8     , 398.80002 , 437.2     ,
         361.80002 ],
        [229.5     , 281.6     , 368.9     , 455.6     , 577.8     ,
         380.7     ],
...
        [223.      , 221.40001 , 237.4     , 238.20001 , 241.9     ,
         290.1     ],
        [236.2     , 237.4     , 232.5     , 215.      , 216.5     ,
         252.8     ],
        [253.6     , 261.19998 , 239.40001 , 203.1     , 191.7     ,
         201.5     ]],

       [[304.8     , 360.2     , 365.90002 , 338.80002 , 360.40002 ,
         412.90002 ],
        [326.6     , 343.1     , 366.      , 358.4     , 398.6     ,
         361.90002 ],
        [313.80002 , 359.9     , 349.3     , 410.5     , 494.4     ,
         404.      ],
        ...,
        [273.3     , 260.2     , 264.3     , 261.2     , 263.9     ,
         315.2     ],
        [258.1     , 251.6     , 240.20001 , 221.1     , 222.4     ,
         256.6     ],
        [259.30002 , 266.30002 , 241.5     , 201.40001 , 188.90001 ,
         197.5     ]]], dtype=float32)
Coordinates:
  * time     (time) datetime64[ns] 1984-12-01 1985-12-01 ... 2013-12-01
  * lon      (lon) float32 27.25 27.75 28.25 28.75 29.25 29.75
  * lat      (lat) float32 -30.75 -30.25 -29.75 -29.25 ... -23.25 -22.75 -22.25

指数数据（nino_djf）

<xarray.DataArray 'ANOM' (time: 30)>
array([-0.94333333, -0.45333333,  1.26333333,  0.84333333, -1.66666667,
        0.16333333,  0.52      ,  1.90333333,  0.26666667,  0.24      ,
        1.14      , -0.75333333, -0.38333333,  2.37      , -1.43      ,
       -1.54666667, -0.62666667, -0.06666667,  1.00666667,  0.44333333,
        0.72333333, -0.7       ,  0.79      , -1.51333333, -0.71      ,
        1.64333333, -1.28333333, -0.73333333, -0.3       , -0.29      ])
Coordinates:
  * time     (time) datetime64[ns] 1984-12-01 1985-12-01 ... 2013-12-01

解决方案 #

错误原因 #

xarray的spearman_r函数默认会对齐两个输入数据的维度索引，重采样后降水数据的time索引出现重复或错位，和指数数据的time索引不匹配，因此抛出错误。而我们需要忽略索引对齐，直接按数据的位置顺序计算相关性。

修正方案（保留xarray向量化优势） #

在重采样后重置降水数据的time坐标，让它和指数数据的time坐标完全一致，强制xarray按位置而非索引值匹配计算：

# Number of bootstrap samples
num_bootstrap_samples = 1000

# Function to calculate spearman correlation
def calculate_correlation(data1, data2):
    return xs.spearman_r(data1, data2, dim='time')

# Calculate the observed correlation
observed_correlation = calculate_correlation(cru_djf_uv, nino_djf)

# Initialize an array to store bootstrap samples
bootstrap_correlations = np.zeros((num_bootstrap_samples, cru_djf_uv.shape[1], cru_djf_uv.shape[2]))

# Perform bootstrap resampling
for i in range(num_bootstrap_samples):
    # Resample with replacement
    resampled_indices = np.random.choice(cru_djf_uv.shape[0], cru_djf_uv.shape[0], replace=True)
    resampled_cru = cru_djf_uv.isel(time=resampled_indices)
    
    # 重置time坐标为指数数据的time序列，跳过索引对齐检查
    resampled_cru = resampled_cru.assign_coords(time=nino_djf.time)
    
    # Calculate correlation for the resampled data
    bootstrap_correlations[i, :, :] = calculate_correlation(resampled_cru, nino_djf).values

备选方案（基于numpy数组，彻底避开索引问题） #

如果想完全脱离xarray的索引机制，可以直接提取数据的numpy数组进行计算：

import numpy as np

num_bootstrap_samples = 1000

# 提取numpy数组，避免索引干扰
cru_data = cru_djf_uv.values
nino_data = nino_djf.values

# 计算观测相关性（保留xarray结构用于后续分析）
observed_correlation = xs.spearman_r(cru_djf_uv, nino_djf, dim='time')

# Initialize an array to store bootstrap samples
bootstrap_correlations = np.zeros((num_bootstrap_samples, cru_djf_uv.shape[1], cru_djf_uv.shape[2]))

# Perform bootstrap resampling
for i in range(num_bootstrap_samples):
    # Resample with replacement
    resampled_indices = np.random.choice(cru_data.shape[0], cru_data.shape[0], replace=True)
    resampled_cru_data = cru_data[resampled_indices, :, :]
    
    # 遍历每个经纬度点计算相关性
    for lat_idx in range(cru_data.shape[1]):
        for lon_idx in range(cru_data.shape[2]):
            rho, _ = xs.spearman_r(resampled_cru_data[:, lat_idx, lon_idx], nino_data)
            bootstrap_correlations[i, lat_idx, lon_idx] = rho

# 可选：转换为xarray DataArray方便后续处理
bootstrap_correlations_da = xr.DataArray(
    bootstrap_correlations,
    dims=['sample', 'lat', 'lon'],
    coords={
        'lat': cru_djf_uv.lat,
        'lon': cru_djf_uv.lon,
        'sample': np.arange(num_bootstrap_samples)
    }
)

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