使用 terra 包计算 R 中两个多层栅格之间的像素相关性

Question

R terra 包中是否有一个函数相当于下面 raster 包的

corLocal()

语法

corLocal(x, y, method=c("spearman"), test=FALSE)

我想计算两个多层 SpatRaster 对象之间的像素相关性。 SpatRaster 具有相同的范围、分辨率和层数。我还需要找到每个像素的 p 值。我找到了

layerCor()

，但它只接受一个 SpatRaster 作为输入。

谢谢

Answer 1

对于 terra 1.7-69，您可以使用

xapp

目前这是 therra 的开发版本，您可以安装：

install.packages('terra', repos='https://rspatial.r-universe.dev')

示例数据

library(terra)
#terra 1.7.69
r <- rast(ncols=10, nrows=10, nlyr=5)
set.seed(1)
r <- init(r, runif)
s <- init(r, runif)

解决方案

x <- xapp(r, s, fun=cor)
#class       : SpatRaster 
#dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source(s)   : memory
#name        :      lyr.1 
#min value   : -0.9923990 
#max value   :  0.9808943 

同样，您可以使用

cor.test

，但您必须从该函数返回的值中提取所需的值。例如

xtest <- xapp(r, s, fun=\(x, y) cor.test(x, y)$p.value)

与

cor

不同，当输入为

cor.test

时，

NA

会失败。要解决这个问题，你可以这样做

naproof.cor.test <- function(x, y) {
    if (any(is.na(x)) || any(is.na(y))) return(NA)
    cor.test(x, y)$p.value
}

xtest <- xapp(r, s, naproof.cor.test)

您可以像这样一步获得相关系数和 p 值

my.cor <- function(x, y) {
    if (any(is.na(x)) || any(is.na(y))) return(c(NA, NA))
    unlist(cor.test(x, y)[c("estimate", "p.value")])
}

z <- xapp(r, s, my.cor)

z
#class       : SpatRaster 
#dimensions  : 10, 10, 2  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#source(s)   : memory
#names       : estimate.cor,      p.value 
#min values  :   -0.9923990, 0.0007945862 
#max values  :    0.9808943, 0.9992525680