在R的terra包中计算两个多层栅格的逐像元相关性

目前terra包没有直接对应raster包corLocal()的内置函数，但可以通过app()函数自定义实现逐像元的相关性及p值计算，完全满足你对两个同范围、同分辨率、同层数SpatRaster对象的分析需求。

实现步骤与代码示例 #

假设你有两个SpatRaster对象x和y，两者层数一致（均为k层）：

1. 合并两个栅格对象
   将x和y合并为一个多层SpatRaster，方便app()逐像元提取所有层的值：

library(terra)
# 构造示例数据（可替换为你的实际数据）
x <- rast(ncols=10, nrows=10, nlyr=5)
values(x) <- rnorm(ncell(x)*nlyr(x))
y <- rast(ncols=10, nrows=10, nlyr=5)
values(y) <- values(x) + rnorm(ncell(x)*nlyr(x), sd=0.2)

# 合并栅格
xy <- c(x, y)
k <- nlyr(x) # 获取层数

2. 定义自定义计算函数
   编写函数处理单个像元的层值对，计算Spearman相关系数与p值，同时处理缺失值和样本量不足的情况：

pixel_cor <- function(v, k) {
  # 拆分两个栅格的层值
  x_vals <- v[1:k]
  y_vals <- v[(k+1):(2*k)]
  
  # 过滤非缺失值对
  valid_pairs <- !is.na(x_vals) & !is.na(y_vals)
  if (sum(valid_pairs) < 2) {
    # 有效样本量不足时返回NA
    return(c(cor=NA_real_, p=NA_real_))
  }
  
  # 计算Spearman相关性与p值
  cor_result <- cor.test(x_vals[valid_pairs], y_vals[valid_pairs], 
                         method="spearman", exact=FALSE)
  return(c(cor=unname(cor_result$estimate), p=cor_result$p.value))
}

3. 逐像元应用函数
   使用app()函数将自定义逻辑应用到每个像元，得到包含相关系数和p值的两层SpatRaster：

# 执行逐像元计算
cor_p_raster <- app(xy, pixel_cor, k=k)

# 重命名结果层，方便识别
names(cor_p_raster) <- c("spearman_cor", "p_value")

结果说明 #

• cor_p_raster[[1]]：存储每个像元的Spearman相关系数
• cor_p_raster[[2]]：存储对应像元的相关性检验p值

你可以根据需求替换cor.test()中的method参数（如"pearson"），实现不同类型的相关性计算。

内容的提问来源于stack exchange，提问作者San-Diego-WX