使用移动窗口计算R中多个栅格图层之间的成对偏相关

Question

现在我有三个栅格数据，分别命名为

y

、

x1

和

x2

。我想使用移动窗口创建 5X5 窗口来计算偏相关系数。最终结果将返回

y

和

x1

的偏相关系数和对应的P值，以及

y

和

x2

的偏相关系数和对应的P值。我在Terra扩展包中找到了

focal

函数，但不知道如何使用它来完成偏相关分析。非常感谢您的帮助。

> library(terra)
> r <- rast(system.file("ex/logo.tif", package="terra"))  
> names(r) <- c("y","x1","x2")
> r
class       : SpatRaster 
dimensions  : 77, 101, 3  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
coord. ref. : Cartesian (Meter) 
source      : logo.tif 
colors RGB  : 1, 2, 3 
names       :   y,  x1,  x2 
min values  :   0,   0,   0 
max values  : 255, 255, 255

> 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)
+   }
+ }
> terra::focal(r,w = c(5,5),fun = rfun)
error: [focal] test failed

> 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"))
> res <- focal(x = r,w = 5,fun = function(x) {
+   d <- na.omit(data.frame(r))
+   if (nrow(d) < 25) {
+     c(NA,NA)
+   }else{
+     pm <- lm(y~x1+x2,data = d)
+     c(broom::tidy(pm)$estimate[1],broom::tidy(pm)$estimate[2])
+   }
+ })
Warning message:
[readStart] source already open for reading 
> res
class       : SpatRaster 
dimensions  : 10, 10, 6  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 10, 0, 10  (xmin, xmax, ymin, ymax)
coord. ref. :  
source(s)   : memory
names       :      lyr1,        lyr2,      lyr3,        lyr4,      lyr5,        lyr6 
min values  : 0.4240551, -0.07680683, 0.4240551, -0.07680683, 0.4240551, -0.07680683 
max values  : 0.4240551, -0.07680683, 0.4240551, -0.07680683, 0.4240551, -0.07680683 

我自己模拟生成了几个网格数据，都是10行10列。我使用焦点函数成功地做到了

OLS

，但是我的系数结果显示每个网格的系数值是相同的。这哪里出了问题？

Answer 1

您可以使用

terra::focalReg

为此

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

rfun <- function(y, x, ...) {
   d <- na.omit(data.frame(y, x))
   # the ncol==2 test is to avoid a bug in testing to find the output number of layers 
   if ((ncol(d) == 2) | (nrow(d) < 25)) {  
     c(NA, NA)
   }else{
     pm <- ppcor::pcor.test(d[,1], d[,2], d[,3])
     c(pm$estimate,pm$p.value)
   }
}

f <- focalReg(r, w=5, fun=rfun)