如何对两个同规格8万层栅格堆栈执行对应层逐元素相加？

问题背景 #

拥有两个terra包的SpatRaster栅格堆栈，各含约80000层，二者范围、CRS、分辨率及层数完全一致，需实现对应层的逐元素相加（即stack3[[n]] = stack1[[n]] + stack2[[n]]），尝试的三种方法均存在问题：

• 循环：逻辑正确但8万层耗时过长
• 直接+运算符：结果不符合预期
• terra::app()触发错误

错误原因分析 #

1. 直接+运算符的问题：
   terra中SpatRaster的算术运算符是对所有层的像素进行求和，而非对应层相加。例如stack1含层A1、A2，stack2含层B1、B2，直接stack1+stack2会计算每个像素的A1+B1+A2+B2，而非生成A1+B1、A2+B2的新堆栈。

2. terra::app()的错误：
   app()仅针对单个SpatRaster的所有层应用函数（如对每个像素的所有层求和），无法直接传入两个栅格对象处理对应层，需使用专门的多栅格对应层处理函数。

可行解决方案 #

方法1：使用terra::lapp()（推荐，高效） #

lapp()是terra中专门处理多个同参数SpatRaster对应层的函数，底层为C++批量实现，效率远高于循环。

代码示例：

library(terra)
# 传入两个栅格堆栈的列表，指定相加操作
stack3 <- lapp(list(stack1, stack2), fun = function(x, y) x + y)
# 更简洁的写法
stack3 <- lapp(list(stack1, stack2), `+`)

方法2：合并栅格后用terra::tapp()分组求和 #

若lapp()因特殊环境受限，可将两个堆栈合并后按对应层分组求和：

library(terra)
# 合并两个堆栈，stack1的层在前，stack2的层在后
combined_rast <- c(stack1, stack2)
# 生成分组索引：每2层为一组（stack1第n层与stack2第n层为一组）
group_indices <- rep(1:nlyr(stack1), 2)
# 按分组对层求和
stack3 <- tapp(combined_rast, group_indices, fun = sum)

方法3：优化循环（应急备选） #

若上述方法均不可用，可通过预分配内存、写入文件优化循环效率：

library(terra)
n_layers <- nlyr(stack1)
# 预先创建与stack1参数一致的空栅格
stack3 <- rast(stack1, nlyr = n_layers)
# 开启文件写入模式，避免内存溢出
writeStart(stack3, filename = "sum_stack.tif", overwrite = TRUE)
for (i in 1:n_layers) {
  layer_sum <- stack1[[i]] + stack2[[i]]
  writeValues(stack3, layer_sum, i)
}
writeStop(stack3)

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