处理大型栅格内存不足：Zonal计算优化及terra参数调整求助

问题描述 #

拥有大型栅格数据，尝试切片生成VRT后执行zonal计算仍触发std::bad_alloc内存错误，需找到可行的分块计算方法，并了解terraOptions参数调整方案。

用户提供的可复现代码：

library(terra)

f <- system.file("ex/lux.shp", package = "terra")
v <- vect(f)
r <- rast(v)
values(r) <- 1:ncell(r)
tester <- zonal(r, v, fun = "max") # 正常运行

r2 <- disagg(r, fact = c(4000, 5000)) # 生成超大栅格
tester <- zonal(r2, v, fun = "max") # 触发内存错误：Error: std::bad_alloc

filename <- "test_tile/test_.tif"
ff <- makeTiles(r2, 10000, filename) # 切片保存
vrtfile <- "test_tile/testing.vrt"
vrt_thing <- vrt(ff, vrtfile, overwrite=T)

tester <- zonal(vrt_thing, v, fun = "max") # 同样触发内存错误

# 当前terra配置
terraOptions()
# memfrac   : 0.6
# tolerance : 0.1
# verbose   : FALSE
# todisk    : FALSE
# tempdir   : C:/Users/vic039/AppData/Local/Temp/Rtmp6nthBI
# datatype  : FLT4S
# memmin    : 1
# progress  : 3

解决方案 #

1. 为什么VRT无法解决问题 #

VRT是虚拟栅格拼接文件，仅记录切片的位置信息，zonal函数处理时依然会尝试将整个虚拟栅格加载到内存计算，因此无法解决内存不足问题，必须手动实现分块计算与结果合并。

2. 手动分块计算Zonal统计 #

通过遍历每个切片，单独计算Zonal统计后合并结果，每次仅加载一小块数据到内存：

library(terra)

# 准备数据（复用示例代码）
f <- system.file("ex/lux.shp", package = "terra")
v <- vect(f)
r <- rast(v)
values(r) <- 1:ncell(r)
r2 <- disagg(r, fact = c(4000, 5000))

# 创建切片保存目录并生成切片
tile_dir <- "test_tile"
dir.create(tile_dir, showWarnings = FALSE)
tile_files <- makeTiles(r2, 10000, file.path(tile_dir, "test_.tif"))

# 初始化结果容器
final_result <- NULL

# 遍历每个切片计算并合并结果
for (tile_path in tile_files) {
  # 加载单个切片
  tile_rast <- rast(tile_path)
  # 裁剪切片到多边形范围，减少无效计算
  cropped_rast <- crop(tile_rast, v)
  # 计算当前切片的Zonal统计
  tile_zonal <- zonal(cropped_rast, v, fun = "max")
  
  # 合并结果：保留每个多边形的最大值
  if (is.null(final_result)) {
    final_result <- tile_zonal
  } else {
    # 按ID合并，取两个结果中的最大值
    merged <- merge(final_result, tile_zonal, by = "ID", suffixes = c("_prev", "_curr"))
    merged$max <- pmax(merged$max_prev, merged$max_curr, na.rm = TRUE)
    final_result <- merged[, c("ID", "max")]
  }
}

# 查看最终结果
print(final_result)

3. 调整terraOptions参数优化内存使用 #

修改以下参数可进一步降低内存压力：

• terraOptions(todisk = TRUE)：强制将计算过程中的中间结果写入磁盘，而非全部占用内存，是处理大型数据的核心设置。
• terraOptions(memfrac = 0.4)：限制terra使用的内存比例（默认0.6），系统内存紧张时可调低至0.3-0.5，避免内存耗尽。
• terraOptions(tempdir = "D:/temp/terra_temp")：指定临时目录到剩余空间充足的磁盘，避免临时文件占满系统盘。
• terraOptions(verbose = TRUE)：开启 verbose 模式，方便监控计算过程中的内存使用与进度。

4. 针对多次Zonal计算的优化 #

由于栅格是静态的，可提前完成切片保存，后续每次搭配不同多边形计算时，直接遍历已有的切片文件即可，无需重复切片，节省时间与资源。

内容的提问来源于stack exchange，提问作者we need a Mat. Stat.