如何基于xarray计算非洲区域降水数据集逐季Spearman秩相关并空间可视化

Spearman秩相关的核心逻辑是先对原始变量在目标维度上做秩转换，再对转换后的秩序列计算Pearson相关系数，你原本的Pearson计算函数不需要做任何修改，仅需要新增一步秩转换的预处理步骤即可实现需求。

核心实现代码 #

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# 原有的Pearson相关计算函数可直接保留使用
def covariance(x, y, dims=None):
    return xr.dot(x - x.mean(dims), y - y.mean(dims), dims=dims) / x.count(dims)

def correlation(x, y, dims=None):
    return covariance(x, y, dims) / (x.std(dims) * y.std(dims))

# 新增：沿指定维度的秩转换函数
def rank_transform(da, dim='time'):
    # 沿指定维度做秩转换，缺失值保持原有状态不参与计算
    return da.rank(dim=dim, na_option='keep')

# 单季节计算示例（对应你原有DJF季节的计算逻辑）
# 先对两个数据集的DJF时段数据做秩转换
da1_djf_rank = rank_transform(da1_sum_season.sel(time=da1_sum_season.time.dt.month.isin([12,1,2])), dim='time')
da2_djf_rank = rank_transform(da2_sum_season.sel(time=da2_sum_season.time.dt.month.isin([12,1,2])), dim='time')
# 对秩转换后的数据计算Pearson相关，得到的结果即为Spearman秩相关系数
ds2_djf_spearman = correlation(da1_djf_rank, da2_djf_rank, dims='time')

# 更高效的批量计算方案：一次性计算四个季节的Spearman相关
# 先对全时序数据做秩转换
da1_rank = rank_transform(da1_sum_season, dim='time')
da2_rank = rank_transform(da2_sum_season, dim='time')
# 按季节分组逐组计算相关
seasonal_spearman = xr.Dataset()
for season, da1_season in da1_rank.groupby('time.season'):
    da2_season = da2_rank.sel(time=da2_rank.time.dt.season == season)
    seasonal_spearman[season] = correlation(da1_season, da2_season, dims='time')

空间可视化示例 #

以DJF季节的相关系数结果为例绘制空间分布图：

fig = plt.figure(figsize=(10,8))
ax = plt.axes(projection=ccrs.PlateCarree())
# 绘制相关系数热力图
cs = ds2_djf_spearman.plot(
    ax=ax,
    transform=ccrs.PlateCarree(),
    cmap='RdBu_r',
    vmin=-1, vmax=1,
    cbar_kwargs={'shrink':0.8, 'label':'Spearman秩相关系数'}
)
# 添加非洲区域的海岸线、网格线
ax.coastlines(resolution='50m', linewidth=0.8)
ax.gridlines(draw_labels=True, linestyle='--', alpha=0.7)
ax.set_title('两个数据集DJF季节降水Spearman秩相关空间分布', fontsize=12)
# 限定绘图范围与你的数据坐标匹配
ax.set_extent([-20.25, 60.25, -40.25, 40.25], crs=ccrs.PlateCarree())
plt.show()

注意事项 #

• 气象领域常用的DJF季节对应12、1、2三个月，你原有代码仅筛选了12月，可根据实际需求调整时间筛选逻辑
• 秩转换时na_option='keep'参数会自动跳过缺失值，与你原有计算逻辑兼容，不会影响无缺测区域的结果
• 计算前需确保两个数据集的时间序列长度、时间戳完全对齐，否则需要先做时间维度的匹配处理

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