基于R语言terra包的掩膜+全局求和工作流提速咨询

问题背景 #

作为terra包新手，我面临一项计算提速难题：手头有50个尺寸为200×200×9000的SpatRaster格式瓦片，需为每个瓦片应用掩膜并计算各图层的全局求和，且要针对6种不同掩膜重复该操作。当前用lapply串行实现的代码可正常运行，但速度较慢，不确定是否值得并行化，恳请提供提速建议。

串行实现代码 #

library(terra)
r10 <- terra::rast(nrow = 10, ncol = 10, nlyrs = 5, vals = 10, xmin = 0, xmax = 1, ymin = 0, ymax = 1)
r20 <- terra::rast(nrow = 10, ncol = 10, nlyrs = 5, vals = 20, xmin = 0, xmax = 1, ymin = 1, ymax = 2)
r30 <- terra::rast(nrow = 10, ncol = 10, nlyrs = 5, vals = 30, xmin = 1, xmax = 2, ymin = 0, ymax = 1)
r40 <- terra::rast(nrow = 10, ncol = 10, nlyrs = 5, vals = 40, xmin = 1, xmax = 2, ymin = 1, ymax = 2)

large.mask <- terra::rast(nrow = 20, ncol = 20, nlyrs = 1, vals = rep(c(1,0), 200), xmin = 0, xmax = 2, ymin = 0, ymax = 2)

rast.list <- list(r10, r20, r30, r40)

compute.layer.sum <- function(x, large.mask, code) {
    a <- terra::crop(large.mask, x)
    b <- terra::match(a, code)
    d <- terra::mask(x, b)
    terra::global(d, sum, na.rm = TRUE)
}

out <- lapply(rast.list, function(x) compute.layer.sum(x, large.mask = large.mask, code = 1))

提速建议 #

1. 提前预处理掩膜，避免重复裁剪 #

每次循环里对大掩膜做crop是重复计算，可提前为每个瓦片、每个掩膜生成对应的裁剪后掩膜：

# 假设有6种掩膜存于masks.list
masks.cropped <- lapply(masks.list, function(mask) {
  lapply(rast.list, function(x) terra::crop(mask, x))
})

后续循环直接调用预处理好的masks.cropped，省去重复裁剪的时间开销。

2. 简化掩膜逻辑，替代match操作 #

原代码中terra::match(a, code)可直接用逻辑判断替代，效率更高：

# 替换原函数中的b <- terra::match(a, code)
b <- a == code

逻辑判断比向量匹配操作更轻量化，尤其针对大尺寸栅格时差异明显。

3. 并行化处理，大幅提升效率 #

完全值得并行化：每个瓦片的处理相互独立，且单瓦片9000图层的计算量不小，6种掩膜进一步放大了计算规模。推荐用future.apply包实现并行：

library(future.apply)
# 设置多进程并行（根据CPU核心数调整workers）
plan(multisession, workers = 4)

# 并行处理单个掩膜
out_parallel <- future_lapply(rast.list, function(x) {
  # 这里用预处理好的裁剪后掩膜
  mask_cropped <- masks.cropped[[1]][[which(rast.list == x)]]
  b <- mask_cropped == 1
  d <- x * b
  terra::global(d, sum, na.rm = TRUE)
})

# 针对6种掩膜，可嵌套并行或循环处理
all_out <- lapply(masks.cropped, function(mask_set) {
  future_lapply(seq_along(rast.list), function(i) {
    x <- rast.list[[i]]
    mask_cropped <- mask_set[[i]]
    b <- mask_cropped == 1
    d <- x * b
    terra::global(d, sum, na.rm = TRUE)
  })
})

注意：用multisession而非multicore，避免terra在多核心模式下的内存共享问题。

4. 替换mask操作，用乘法简化计算 #

原代码中terra::mask会将掩膜外的区域设为NA，可直接用栅格乘法替代，避免NA值处理的额外开销：

# 替换原函数中的mask步骤
d <- x * (a == code)

掩膜内区域为1，乘法后保留原数值；掩膜外为0，求和时不影响结果，比mask后再求和更高效。

5. 优化内存管理 #

处理大尺寸栅格时，及时清理临时对象并触发垃圾回收，避免内存累积：

compute.layer.sum <- function(x, mask_cropped, code) {
  b <- mask_cropped == code
  d <- x * b
  res <- terra::global(d, sum, na.rm = TRUE)
  # 清理临时对象
  rm(b, d)
  gc()
  return(res)
}

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