基于raster或terra实现栅格计算多结果输出

问题背景 #

我正在对一系列超大栅格运行线性拟合模型，需要从模型中返回系数、对应p值等两个及以上结果。在非栅格场景下，我会编写函数以列表形式返回所需结果，但在栅格数据处理中，使用raster包的calc或terra包的app时，当前方法需要重复执行相同计算提取不同结果，效率极低。

非栅格场景的示例代码：

#quick lm function to pull out coefficient and p.value
lm_fun <- function(x, y){
  lm <- lm(x ~ y)
  coef1 <- coef(summary(lm))[2,1]
  pval <- coef(summary(lm))[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)

当前低效的栅格实现：

# 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)

高效解决方案 #

使用terra包实现（推荐，性能更优） #

修改自定义函数，让它返回包含多个结果的向量，app函数会自动将这些向量转换为多图层的SpatRaster：

library(terra)

# 创建示例栅格栈
r1 <- r2 <- r3 <- rast(nrow=100, ncol=100)
values(r1) <- runif(ncell(r1), 2, 5)
values(r2) <- runif(ncell(r2), 4, 10)
values(r3) <- runif(ncell(r3), 6, 14)
s_rast <- c(r1, r2, r3)

timestep <- 1:nlyr(s_rast)

# 自定义函数：一次计算返回系数和p值
lm_multi <- function(x) {
  if (any(is.na(x))) {
    return(c(NA, NA))
  }
  model <- lm(x ~ timestep)
  sum_mod <- summary(model)
  coef_val <- sum_mod$coefficients[2, 1]
  p_val <- sum_mod$coefficients[2, 4]
  return(c(trend = coef_val, p_value = p_val))
}

# 应用函数，直接得到包含两个图层的SpatRaster
result_stack <- app(s_rast, lm_multi)

# 命名图层并查看结果
names(result_stack) <- c("trend_coef", "p_value")
plot(result_stack)

使用raster包实现 #

同样思路，让函数返回向量，calc会自动生成RasterStack：

library(raster)

# 创建示例栅格栈
r1 <- r2 <- r3 <- raster(nrow=100, ncol=100)
values(r1) <- runif(ncell(r1), 2, 5)
values(r2) <- runif(ncell(r2), 4, 10)
values(r3) <- runif(ncell(r3), 6, 14)
s <- stack(r1, r2, r3)

timestep <- 1:nlayers(s)

# 自定义函数
lm_multi <- function(x) {
  if (any(is.na(x))) {
    return(c(NA, NA))
  }
  model <- lm(x ~ timestep)
  sum_mod <- summary(model)
  c(sum_mod$coefficients[2,1], sum_mod$coefficients[2,4])
}

# 计算得到栅格栈
result_stack <- calc(s, lm_multi)
names(result_stack) <- c("trend_coef", "p_value")
plot(result_stack)

关键说明 #

• 核心优化：只执行一次线性拟合，同时提取所有需要的结果，避免重复建模带来的性能损耗，尤其适合超大栅格数据。
• 函数返回值必须是向量（而非列表），这样app/calc才能自动将每个单元格的多结果映射为不同的栅格图层。
• 增加了对NA值的处理，确保存在缺失值的单元格返回对应的NA结果，避免运行报错。

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