在R语言中如何计算两个全球网格数据集的空间相关性？

• 直接运行cor(data1, data2)返回NA：cor()仅支持二维矩阵输入，三维数组传入会被强制展平，且默认参数use = "everything"只要存在NA值就直接返回NA。
• 单月计算返回1080*1080矩阵：cor()默认对输入矩阵的列计算两两相关，你传入的单月数据为2160行（经度）、1080列（纬度）的矩阵，因此输出的是1080个纬度列的两两相关结果，和你要的两组数据相关性完全不符。

根据你的需求（定位不同区域的模拟效果优劣），推荐优先计算每个空间格点的时间维度相关系数，输出结果为2160*1080的矩阵，每个值对应对应经纬度位置的建模与观测匹配度。

# 计算每个经纬度格点的时间序列相关
spatial_cor <- apply(
  # 合并两个三维数组为四维数组，最后一维区分建模/观测
  X = array(c(data1, data2), dim = c(2160, 1080, 12, 2)),
  # 在前两个维度（经度、纬度）上迭代计算
  MARGIN = c(1, 2),
  FUN = function(grid_vals) {
    # 单个格点的12个月建模值和观测值计算相关，自动处理成对NA
    cor(grid_vals[,1], grid_vals[,2], use = "pairwise.complete.obs")
  }
)

如果需要看每个月的全局整体模拟效果，可以计算单月全格点的空间相关：

# 计算12个月各自的全局空间相关系数
monthly_global_cor <- sapply(1:12, function(month) {
  cor(
    as.vector(data1[,,month]), 
    as.vector(data2[,,month]), 
    use = "pairwise.complete.obs"
  )
})

• spatial_cor矩阵的每个值范围为-1到1，越接近1说明对应经纬度区域的模拟效果越好，你可以自行设置阈值划分优劣等级，比如>0.7为表现优秀、<0.3为表现较差，直接映射到全球网格即可可视化定位目标区域。
• 若对数据质量要求更高，可以把use参数改为"complete.obs"，仅当该格点12个月的建模、观测值都无NA时才计算相关，返回的NA会更多但结果更严谨。

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