验证Python实现时间序列与网格化数据Spearman相关性的代码逻辑

代码逻辑验证 #

你提供的代码核心逻辑是正确的，完全符合将网格化气候数据与自定义时间序列做Spearman秩相关分析的需求：

1. 时间序列与NetCDF网格化数据的读取逻辑无误；
2. 遍历每个网格点，提取该点的时间序列与输入气候指数时间序列，通过scipy.stats.spearmanr计算秩相关系数（rho）和显著性p值，这正是目标功能的核心；
3. 将p值≥0.1的区域设为NaN，只保留90%置信度的显著相关区域，显著性筛选逻辑正确；
4. 用热力图展示显著相关系数的空间分布，基础绘图逻辑合理。

关键前提校验 #

代码能输出有效结果的核心前提是：
输入的september.txt时间序列长度，必须与NetCDF数据中msl变量的时间维度长度完全一致。比如你提到时间序列有44个值，那么mslp_data['msl'].sizes['time']必须等于44，否则会因维度不匹配导致计算错误或结果无意义。建议在代码开头添加校验：

if len(time_series) != mslp_data['msl'].sizes['time']:
    raise ValueError("时间序列长度与网格化数据的时间维度不匹配，请检查数据！")

代码优化建议 #

1. 替换双重循环，提升计算效率 #

原代码的双重循环在网格点数量较多时（如全球格点）会非常缓慢，推荐用xarray的apply_ufunc实现向量化计算，大幅提升速度：

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from scipy.stats import spearmanr

# 读取数据
time_series = np.loadtxt("september.txt")
mslp_data = xr.open_dataset("area.nc")

# 时间维度校验
if len(time_series) != mslp_data['msl'].sizes['time']:
    raise ValueError("时间序列长度与网格化数据的时间维度不匹配，请检查数据！")

# 定义Spearman计算函数
def spearman_corr(grid_ts, index_ts):
    rho, p_value = spearmanr(grid_ts, index_ts, nan_policy='omit')
    return rho, p_value

# 向量化计算所有网格点的相关系数与p值
corr, p_val = xr.apply_ufunc(
    spearman_corr,
    mslp_data['msl'],
    time_series,
    input_core_dims=[['time'], []],  # 指定两个输入的核心维度
    vectorize=True,  # 自动对非核心维度（lat, lon）循环
    output_core_dims=[[], []]  # 两个输出对应lat, lon的二维数组
)

# 筛选显著相关区域
significant_corr = corr.where(p_val < 0.1)

# 用xarray自带绘图，自动匹配经纬度
significant_corr.plot(
    figsize=(10,5),
    cmap='coolwarm',
    cbar_kwargs={'label': 'Spearman Correlation Coefficient'},
    title='Temporal Correlation with Mean Sea Level Pressure (Significant at 90%)'
)
plt.show()

2. 缺失值处理 #

如果网格化数据中存在缺失值（NaN），在spearmanr中指定nan_policy='omit'，可自动忽略缺失的时间点对，避免计算报错。

3. 专业地理绘图 #

若需要更规范的地理投影（如显示海岸线、真实经纬度），可结合cartopy库，示例片段：

import cartopy.crs as ccrs

fig, ax = plt.subplots(figsize=(10,5), subplot_kw={'projection': ccrs.PlateCarree()})
significant_corr.plot(
    ax=ax,
    cmap='coolwarm',
    cbar_kwargs={'label': 'Spearman Correlation Coefficient'},
    transform=ccrs.PlateCarree()
)
ax.coastlines()
ax.set_title('Temporal Correlation with Mean Sea Level Pressure (Significant at 90%)')
plt.show()

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