基于文件名匹配对多栅格执行栅格减法运算

有两组栅格数据，想要通过栅格运算将其中一组减去另一组，但无法正确匹配对应的栅格对。已知用terra库可以轻松实现单个栅格对的减法：

library(terra)
s <- rast(ncols = 22, nrows = 25, nlyrs = 1, xmin = 0, vals = 10)
p <- rast(ncols = 22, nrows = 25, nlyrs = 1, xmin = 0, vals = 8)
r <- s - p

但文件夹中有大量栅格文件，需要根据文件名末尾的匹配字符串配对执行减法（比如用above_b.tif减去below_b.tif），文件命名格式如下：

above_a.tif
above_b.tif
above_c.tif

below_a.tif
below_b.tif
below_c.tif

已经通过以下代码获取两类文件列表，但无法实现按文件名匹配的栅格减法：

list_above <- list.files(pattern = 'above', all.files = T, full.names = F)
list_below <- list.files(pattern = 'below', all.files = T, full.names = F) 

尝试的减法代码无法满足需求：

temp_sub <- lapply(list_above, function(x) {
  to_subtract <- list_below
  subtracted <- x - to_subtract
  subtracted
})

请问如何实现基于文件名匹配的批量栅格减法？是否需要将栅格整理为有序的栅格栈？

方法1：按文件名标识匹配处理 #

核心是提取每个文件名的匹配标识（如_a、_b），为每个above文件找到对应的below文件后执行减法：

library(terra)

# 获取带完整路径的文件列表，方便直接读取
list_above <- list.files(pattern = 'above', full.names = TRUE)
list_below <- list.files(pattern = 'below', full.names = TRUE)

# 定义函数提取文件名中的匹配标识
get_match_id <- function(file_path) {
  # 从文件名中去除前缀和后缀，得到核心匹配标识
  sub("(above_|below_|\\.tif)", "", basename(file_path))
}

# 遍历每个above文件，匹配对应below文件并计算
for (above_file in list_above) {
  current_id <- get_match_id(above_file)
  # 找到对应标识的below文件
  matched_below <- list_below[grepl(current_id, list_below)]
  
  # 读取栅格并执行减法
  r_above <- rast(above_file)
  r_below <- rast(matched_below)
  r_diff <- r_above - r_below
  
  # 保存结果，命名格式为diff_标识.tif
  writeRaster(r_diff, filename = paste0("diff_", current_id, ".tif"), overwrite = TRUE)
}

方法2：栅格栈批量运算（高效版） #

如果所有栅格的空间范围、分辨率、投影完全一致，可以将两类栅格分别整理为有序栅格栈，直接执行栈间减法：

library(terra)

# 对文件列表排序，确保两类栅格的顺序严格对应
list_above <- sort(list.files(pattern = 'above', full.names = TRUE))
list_below <- sort(list.files(pattern = 'below', full.names = TRUE))

# 创建栅格栈
stack_above <- rast(list_above)
stack_below <- rast(list_below)

# 栈间减法：每个对应图层自动相减
stack_diff <- stack_above - stack_below

# 按图层保存结果，每个图层对应一个文件
writeRaster(stack_diff, 
            filename = paste0("diff_", 1:nlyr(stack_diff), ".tif"), 
            overwrite = TRUE, 
            bylayer = TRUE)

关键注意点 #

• 方法1不依赖文件排序，只要文件名标识匹配即可，适合文件顺序混乱的场景；
• 方法2要求两类栅格数量一致，且排序后顺序完全对应，运算效率更高；
• 执行前务必确认所有栅格的空间参考、分辨率、范围完全一致，否则terra会自动重采样（可能不符合预期）或报错。

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