如何用R语言从匹配层名的多层栅格列表计算对应层均值栅格栈？

我有一组范围、图层名称完全匹配的多层栅格（SpatRaster栈），需要对对应图层分别计算均值，最终得到一个保留原图层结构的单个栅格栈。手动执行terra::mean(stack.1, stack.2, na.rm = T)可以实现，但现在需要批量处理大量文件，急需通用方案。

尝试过用terra::rast()将列表转为单个栈后取均值，结果会把所有图层平均成单层，没法保留对应层；用lapply()和do.call()只能对列表元素单独操作，完成不了整体的对应层均值计算，求解决办法。

代码示例 #

# 栅格栈结构
stack.1
class       : SpatRaster 
dimensions  : 776, 720, 6  (nrow, ncol, nlyr)
resolution  : 0.6, 0.6  (x, y)
extent      : 278110.2, 278542.2, 4752927, 4753393  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 16N (EPSG:32616) 
source      : pepr.2022.100.tif 
names       : preds~_temp, preds~_temp, preds~_temp, preds~vapor, preds~vapor, preds~vapor 
min values  :   -2.074884,    29.46503,    7.464921,   -118.1709,    2325.729,    984.7896 
max values  :    1.352281,    43.14127,   11.517185,    128.4751,    2929.794,   1289.9479

# 将多个栅格栈存入列表
stack.list <- list(stack.1, stack.2, stack.3)

# 期望结果：保留原图层结构的均值栅格栈
mean(stack.1, stack.2, stack.3, na.rm = TRUE)
class       : SpatRaster 
dimensions  : 776, 720, 6  (nrow, ncol, nlyr)
resolution  : 0.6, 0.6  (x, y)
extent      : 278110.2, 278542.2, 4752927, 4753393  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 16N (EPSG:32616) 
source(s)   : memory
names       : preds~_temp, preds~_temp, preds~_temp, preds~vapor, preds~vapor, preds~vapor 
min values  :   -1.928547,    29.14394,    5.955294,   -115.4321,    2044.692,    941.3278 
max values  :    1.215450,    38.19327,    9.235113,    127.4777,    2525.456,   1166.1364

方法1：用do.call()直接传递列表参数 #

terra::mean()本身支持接收多个SpatRaster对象作为参数，用do.call()可以把列表中的元素拆分为单独参数传入，完全等价于手动输入多个栈的操作，是最简洁的方案：

mean_stack <- do.call(terra::mean, c(stack.list, list(na.rm = TRUE)))

该方法会自动对应所有栈的同名/同位置图层计算均值，直接保留原图层的名称和结构。

方法2：terra::app()结合图层索引遍历 #

如果需要更灵活的自定义操作，可以按图层索引批量提取对应图层计算均值：

# 获取单个栈的图层数量
n_layers <- nlyr(stack.list[[1]])

# 遍历每个图层索引，提取所有栈的对应图层并计算均值
mean_stack <- rast(lapply(1:n_layers, function(i) {
  layers_to_avg <- lapply(stack.list, function(x) x[[i]])
  app(rast(layers_to_avg), mean, na.rm = TRUE)
}))

# 还原原图层名称
names(mean_stack) <- names(stack.list[[1]])

方法3：利用terra::sds()处理多维度栅格 #

将列表转为SpatRasterDataset后，对数据集维度（即各个栅格栈）计算均值：

# 转为SpatRasterDataset
sds_obj <- sds(stack.list)
# 对第三个维度（多组栅格栈）计算均值
mean_stack <- app(sds_obj, mean, na.rm = TRUE)

此方法适合需要对多组栅格进行维度级操作的场景，结果同样保留原图层结构。

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