在R中使用Terra包滑动窗口计算多栅格偏相关系数的问题

一、terra::focal运行报错的修复 #

报错核心原因是函数返回值数量与默认输出设置不匹配：你的自定义函数rfun返回2个值（偏相关系数、P值），但focal默认要求返回值数量等于输入栅格的图层数（3层）。只需添加expand=TRUE参数，告知focal扩展输出为对应返回值数量的图层即可。

修正后的完整代码：

library(terra)
library(ppcor)

r <- rast(system.file("ex/logo.tif", package="terra"))  
names(r) <- c("y","x1","x2")

rfun <- function(r) {
  d <- na.omit(data.frame(r))
  if (nrow(d) < 25) {
    c(NA, NA)
  } else {
    pm <- ppcor::pcor.test(d[,1], d[,2], d[,3])
    c(pm$estimate, pm$p.value)
  }
}

# 添加expand=TRUE参数以支持多值返回
result <- terra::focal(r, w = c(5,5), fun = rfun, expand = TRUE)
names(result) <- c("pcor_estimate", "p_value")
result

二、三种方法结果不一致的原因 #

Method3错误使用了terra::focalReg函数：

• focalReg是专门用于滑动窗口回归分析的工具，要求传入的fun是回归模型函数（如lm），而非自定义的偏相关计算函数。强行传入rfun会导致输入逻辑混乱，计算结果完全偏离预期。
• Method1和Method2是直接提取目标窗口的正确数据进行计算，结果一致；用修正后的focal方法替代Method3，即可得到与前两种方法匹配的滑动窗口结果。

验证修正后的结果：

set.seed(20241001)
r1 <- rast(matrix(runif(100),10,10))
r2 <- rast(matrix(runif(100),10,10))
r3 <- rast(matrix(runif(100),10,10))
r <- c(r1,r2,r3) %>% setNames(c("y","x1","x2"))

# 用修正后的focal计算滑动窗口结果
result_focal <- terra::focal(r, w=5, fun=rfun, expand=TRUE)
# 提取(3,3)位置的结果，与Method1/2一致
as.matrix(result_focal$lyr1)[3,3] # 偏相关系数
as.matrix(result_focal$lyr2)[3,3] # P值

额外说明 #

• 若需要计算y与x2的偏相关（控制x1），只需将rfun中的pcor.test参数修改为ppcor::pcor.test(d[,1], d[,3], d[,2])。
• 边缘窗口（像素数不足25）会返回NA，符合你函数中设置的逻辑。

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