Terra处理超大栅格触发std::bad_alloc，求区域统计可行方案

问题根源：对VRT的误解 #

你之前的思路误区在于：VRT是虚拟拼接的栅格，它并没有实现"分块处理再自动合并"的逻辑。当调用zonal时，Terra依然会尝试将VRT对应的完整栅格数据加载到内存，因此不管是原大栅格还是VRT，都会触发内存溢出错误。

可行解决方案：分块计算+结果合并 #

核心思路是每次仅处理栅格的一小部分，计算该部分的区域统计值，再逐步合并得到全局结果。下面提供两种高效实现方式：

方式1：基于已生成的瓦片文件处理

如果你已经提前切割了瓦片，可以遍历每个瓦片计算局部统计，再汇总全局的range（需要保留每个区域的全局最小/最大值，最终计算差值）：

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))

# 生成瓦片（你的原有步骤）
dir.create("test_tile", showWarnings = FALSE)
filename <- "test_tile/test_.tif"
ff <- makeTiles(r2, 1000, filename)

# 初始化结果：记录每个区域的当前全局最小/最大值
result <- data.frame(ID = v$ID, min = Inf, max = -Inf)

# 遍历所有瓦片
for (tile_file in ff) {
  tile_rast <- rast(tile_file)
  # 计算当前瓦片内的区域最小/最大值
  tile_zonal <- zonal(tile_rast, v, fun = function(x) c(min(x), max(x)), na.rm = TRUE)
  colnames(tile_zonal) <- c("ID", "tile_min", "tile_max")
  
  # 合并并更新全局极值
  result <- merge(result, tile_zonal, by = "ID")
  result$min <- pmin(result$min, result$tile_min)
  result$max <- pmax(result$max, result$tile_max)
  result <- result[, c("ID", "min", "max")]
}

# 计算最终的range值
result$range <- result$max - result$min

方式2：直接分块读取原栅格（更高效，无需提前切瓦片）

利用Terra的blockSize函数自动计算合理的分块大小，直接在内存中分块处理原栅格，节省磁盘空间：

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))

# 计算合理的分块大小
bs <- blockSize(r2)
# 初始化结果
result <- data.frame(ID = v$ID, min = Inf, max = -Inf)

# 遍历每个数据块
for (i in 1:bs$n) {
  # 读取当前块
  block <- readStart(r2, bs$row[i], bs$nrows[i])
  # 计算块内的区域最小/最大值
  block_zonal <- zonal(block, v, fun = function(x) c(min(x), max(x)), na.rm = TRUE)
  colnames(block_zonal) <- c("ID", "tile_min", "tile_max")
  
  # 合并并更新全局极值（处理块中无对应区域像素的NA情况）
  result <- merge(result, block_zonal, by = "ID", all.x = TRUE)
  result$tile_min[is.na(result$tile_min)] <- Inf
  result$tile_max[is.na(result$tile_max)] <- -Inf
  
  result$min <- pmin(result$min, result$tile_min)
  result$max <- pmax(result$max, result$tile_max)
  result <- result[, c("ID", "min", "max")]
  
  # 释放当前块的内存
  readStop(r2)
}

# 计算最终range
result$range <- result$max - result$min

复用函数（适配多次不同多边形计算）

将分块逻辑封装成函数，后续只需传入不同的多边形矢量即可重复计算：

big_zonal_range <- function(big_rast, poly_vect) {
  bs <- blockSize(big_rast)
  result <- data.frame(ID = poly_vect$ID, min = Inf, max = -Inf)
  
  for (i in 1:bs$n) {
    block <- readStart(big_rast, bs$row[i], bs$nrows[i])
    block_zonal <- zonal(block, poly_vect, fun = function(x) c(min(x), max(x)), na.rm = TRUE)
    colnames(block_zonal) <- c("ID", "tile_min", "tile_max")
    
    result <- merge(result, block_zonal, by = "ID", all.x = TRUE)
    result$tile_min[is.na(result$tile_min)] <- Inf
    result$tile_max[is.na(result$tile_max)] <- -Inf
    
    result$min <- pmin(result$min, result$tile_min)
    result$max <- pmax(result$max, result$tile_max)
    result <- result[, c("ID", "min", "max")]
    
    readStop(big_rast)
  }
  
  result$range <- result$max - result$min
  return(result)
}

# 使用示例
tester1 <- big_zonal_range(r2, v)
# 传入其他多边形矢量v2
# tester2 <- big_zonal_range(r2, v2)

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