使用 raster 或 terra 从栅格计算中返回多个输出

Question

我正在一系列非常大的栅格上运行线性拟合模型。我经常需要从模型返回两个或多个值 - 例如系数及其关联的 p 值。在非光栅情况下，我会创建一个函数，以列表的形式返回我想要的所有值，如下所示：

#quick lm function to pull out coefficient and p.value
lm_fun <- function(x, y){
  lm <- lm(x ~ y)
  coef1 <- coef(summary(test))[2,1]
  pval <- coef(summary(test))[2,4]
  lm_results <- list(coef1 = coef1, pval = pval)
}

#sample data
dat <- c(runif(100, 2, 5), runif(100, 4, 10), runif(100, 6, 14))
timestep <- 1:length(dat)

lm_results <- lm_fun(dat, timestep)

但是，在处理栅格数据时，自定义函数通过 calc （使用 raster 包）或 app （使用 terra 包）应用于栅格，并返回栅格或 SpatRaster 图层。无论如何，是否可以修改函数以生成栅格堆栈，并在不同层中具有不同的输出（例如系数、p 值）？现在，我正在使用以下方法，但这意味着我多次运行相同的计算来提取结果的不同部分，这在计算上效率极低。我愿意使用栅格或 terra。

# create three identical RasterLayer objects
r1 <- r2 <- r3 <- raster(nrow=100, ncol=100)
# Assign random cell values
values(r1) <- runif(ncell(r1), min=2, max=5)
values(r2) <- runif(ncell(r2), min=4, max=10) 
values(r3) <- runif(ncell(r3), min=6, max=14)
# combine three RasterLayer objects into a RasterStack
s <- stack(r1, r2, r3)
plot(s)

#calculate number of timesteps
nsteps <- dim(s)[3]
timestep <- 1:nsteps

#functions to calculate linear trend and associated p-value
trendfun <- function(x) { if (is.na(x[1])){ NA } else { lm(x ~ timestep)$coefficients[2]}}
pfun <- function(x) { if (is.na(x[1])){ NA } else { summary(lm(x ~ timestep))$coefficients[2,4]}}

# calculate trend and p-value using raster
library(raster)
s_trend <- calc(s, trendfun)
s_trend_pv <- calc(s, pfun)

#alternatively, calculate using terra package
library(terra)
s_rast <- rast(s)
s_trend <- app(s_rast, fun=trendfun)
s_pvalue <- app(s_rast, fun=pfun)

Answer 1

您可以使用

terra::app

。该函数应返回一个向量或矩阵（每个单元格对应一行）。这是一个插图。

set.seed(123)
r <- terra::rast(nrow=10, ncol=10, nlyr=3, vals=runif(300))
timestep <- 1:nlyr(r)

f <- function(d) {
    if (any(is.na(d))) {
        c(NA, NA)
    } else {
        m <- lm(d ~ timestep)
        cf <- summary(m)$coefficients
        cf[2, c(1,4)] 
    }
}

terra::app(r, f)
#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,    Pr(>|t|) 
#min values  : -0.4600143, 0.005280976 
#max values  :  0.4173450, 0.990162160 

另见

terra::regress

terra::regress(r, 1:nlyr(r))
#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       : (Intercept),          x 
#min values  :  -0.3435626, -0.4600143 
#max values  :   1.4499354,  0.4173450