如何使用R语言raster包计算降水栅格的空间相关性？

完全可以借助R语言raster包实现你需要的空间相关性分析，最终输出符合要求的相关性栅格图层。

实现逻辑 #

你的需求核心是逐像元计算指定年段内，目标研究点位的年降水序列和周边每个栅格像元年降水序列的关联指标，最后掩膜掉阈值范围内的像元即可，和参考成果的生成逻辑完全匹配。

注：常规皮尔逊相关系数取值范围为[-1,1]，参考成果设置-3 < 相关值 < 3的区域不显示，推测原作者输出的是相关系数显著性检验的t统计量（t值绝对值>3大致对应p<0.01的显著性水平），你可以根据自己的指标需求调整计算逻辑。

参考实现代码 #

# 加载raster包
library(raster)

# ----------------------
# 步骤1：导入基础数据
# ----------------------
# 替换为你的年降水栅格存放路径，栅格文件建议按年份排序，每个图层对应1年的年降水量
prcp_path <- "你的年降水栅格文件夹路径"
prcp_file_list <- list.files(prcp_path, pattern = "\\.tif$", full.names = TRUE)
prcp_stack <- stack(prcp_file_list)

# 替换为你的目标研究点位坐标，坐标系要和降水栅格保持一致
# 示例坐标为阿克恰莱的大致WGS84经纬度
target_loc <- SpatialPoints(cbind(38.91, 36.71), proj4string = CRS(proj4string(prcp_stack)))

# 提取目标点位的年降水时间序列
target_prcp_ts <- extract(prcp_stack, target_loc)[1, ]

# ----------------------
# 步骤2：逐像元计算关联指标
# ----------------------
# 自定义计算函数，以下提供两种计算逻辑，按需选择
calc_corr <- function(cell_ts) {
  # 缺失值占比超过10%的像元直接返回NA
  if (sum(is.na(cell_ts)) > length(cell_ts)*0.1) return(NA)
  # 逻辑1：计算皮尔逊相关系数
  # cor_val <- cor(cell_ts, target_prcp_ts, use = "complete.obs")
  # 逻辑2：计算相关显著性t统计量（匹配参考图的±3阈值逻辑）
  cor_test_res <- cor.test(cell_ts, target_prcp_ts, use = "complete.obs")
  cor_val <- cor_test_res$statistic
  return(cor_val)
}

# 栅格计算生成关联指标图层
corr_raster <- calc(prcp_stack, fun = calc_corr)

# ----------------------
# 步骤3：按规则掩膜+输出结果
# ----------------------
# 将-3~3之间的像元设为NA，即不显示
corr_raster_final <- mask(corr_raster, corr_raster > -3 & corr_raster < 3, maskvalue = TRUE)

# 输出栅格文件，替换为你自己的输出路径
writeRaster(corr_raster_final, "降水相关性结果.tif", overwrite = TRUE)

参考成果示例 #

补充说明 #

• 如果你的降水数据是NC等其他格式，raster包也支持直接读取构建栅格栈，无需提前转成tiff格式。
• 如果研究区范围大、栅格分辨率高，计算速度慢，可以替换为terra包处理，语法和raster基本一致，计算效率提升明显。

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