使用R语言terra包合并大栅格时R崩溃的解决及内存估算问题

我有多个大尺寸TIFF栅格文件，需要合并为一个大栅格。尝试用R的terra包，先通过vrt()创建.vrt文件，再用writeRaster()写入最终栅格，但执行writeRaster()时R会话崩溃。

尝试过的代码 #

第一段代码 #

library(terra)
test1 <- rast('I:/DATA/output/MI_01.tif')
test2 <- rast('I:/DATA/output/MI_02.tif')
test3 <- rast('I:/DATA/output/MI_03.tif')
v <- vrt(c("I:/DATA/output/MI_01.tif", 
           "I:/DATA/output/MI_02.tif",
           "I:/DATA/output/MI_03.tif"),
         "test.vrt", overwrite=T)
plot(v)
writeRaster(v, "I:/DATA/output/MI_Merged.tif")

第二段代码 #

library(terra)
t.lst <- list.files("/lustre1/scratch/output/MI_tile/", pattern=".tif", full.names=TRUE)
MIs <- function(t.lst, fout="") {
  r <- vrt(t.lst)
  if (fout != "") {
    writeRaster(r, fout, overwrite=TRUE)
    fout
  } else {
    wrap(r)
  }
} 

MIs(t.lst, fout = "/lustre1/scratch/output/MI_final/MI_merge.tif")

待合并栅格信息 #

> test1
class       : SpatRaster 
dimensions  : 120000, 132818, 1  (nrow, ncol, nlyr)
resolution  : 25, 25  (x, y)
extent      : 2636075, 5956525, 1385925, 4385925  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source      : MI_01.tif 
name        : spat_o6xCi6omAQKQWfp_896553 
min value   :                           0 
max value   :                           1 
> test2
class       : SpatRaster 
dimensions  : 41203, 132818, 1  (nrow, ncol, nlyr)
resolution  : 25, 25  (x, y)
extent      : 2636075, 5956525, 4385925, 5416000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source      : MI_02.tif 
name        : spat_Qs2AebAaSIInfQI_2501483 
min value   :                            0 
max value   :                            1 
> test3
class       : SpatRaster 
dimensions  : 161203, 132818, 1  (nrow, ncol, nlyr)
resolution  : 25, 25  (x, y)
extent      : 2636075, 5956525, 1385925, 5416000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source      : MI_03.tif 
name        : spat_U51IfANH30JHuwp_503810 
min value   :                           0 
max value   :                           1 

运行环境与错误信息 #

在HPC（超级计算机）上运行，可用内存约750GiB，预估最终栅格约5GB，但返回段错误：

terra 1.7.55
|---------|---------|---------|---------|
 *** caught segfault ***
address 0x152a3df7ed20, cause 'memory not mapped'

Traceback:
 1: .External(list(name = "CppMethod__invoke_notvoid", address = <pointer: 0x2ee8190>,     dll = list(name = "Rcpp", path = "/x/miniconda3-new/envs/R_new/lib/R/library/Rcpp/libs/Rcpp.so",         dynamicLookup = TRUE, handle = <pointer: 0x43328d0>,         info = <pointer: 0x38933d0>), numParameters = -1L), <pointer: 0x7299890>,     <pointer: 0x3038160>, .pointer, ...)
 2: x@cpp$writeRaster(opt)
 3: .local(x, filename, ...)
 4: writeRaster(r, fout, overwrite = TRUE)
 5: writeRaster(r, fout, overwrite = TRUE)
 6: MIs(t.lst, fout = "/lustre1/scratch/output/MI_final/MI_final.tif")
An irrecoverable exception occurred. R is aborting now ...
/var/spool/slurmd/job61846166/slurm_script: line 19: 1794009 Segmentation fault      (core dumped) Rscript BuildVRTmap.R

具体问题 #

1. 如何在R中以低内存消耗合并这些栅格，避免R会话崩溃？
2. 使用vrt函数构建并写入大栅格时，如何估算所需内存？

问题1：低内存合并栅格的方法 #

方法1：直接用merge()分块处理 #

terra的merge()支持分块读取源栅格，无需加载全部数据到内存，跳过VRT步骤直接合并：

library(terra)
t.lst <- list.files("/lustre1/scratch/output/MI_tile/", pattern=".tif$", full.names=TRUE)
# 读取第一个栅格作为合并基准
base_rast <- rast(t.lst[1])
# 逐个合并剩余栅格，每次合并后释放临时内存
for (i in 2:length(t.lst)) {
  current_rast <- rast(t.lst[i])
  base_rast <- merge(base_rast, current_rast)
  rm(current_rast)
  gc()
}
writeRaster(base_rast, "/lustre1/scratch/output/MI_final/MI_merge.tif", overwrite=TRUE)

方法2：优化VRT写入参数 #

如果坚持用VRT，在writeRaster()中指定datatype和tilesize，强制分块写入，避免一次性申请大量内存：

v <- vrt(t.lst)
# 数据是0-1，用1字节无符号整型足够；设置分块大小控制单次写入内存
writeRaster(v, "/lustre1/scratch/output/MI_final/MI_merge.tif", 
            overwrite=TRUE, 
            datatype="INT1U",
            tilesize=c(1024,1024),
            progress="text")

方法3：调用GDAL命令行工具 #

HPC上直接用GDAL工具处理更稳定，完全绕过R的内存限制：

# 用gdal_merge.py合并，指定输出格式和数据类型
gdal_merge.py -o /lustre1/scratch/output/MI_final/MI_merge.tif -ot Byte -of GTiff /lustre1/scratch/output/MI_tile/*.tif

将这段命令写入脚本提交到HPC任务队列即可。

问题2：VRT写入时的内存估算方法 #

VRT本身是虚拟栅格，几乎不占内存，但writeRaster()时需要加载源栅格的对应块拼接，内存需求看以下几点：

1. 单块数据大小：tilesize宽×高×单像素字节数。比如你的数据用Byte类型（1字节），tilesize设为2048×2048，单块大小约4MB；如果默认全幅写入，未压缩的最终栅格总大小是161203×132818×1≈21GB（你之前预估的5GB是压缩后大小）。
2. 同时加载的块数：合并时某区域有N个源栅格重叠，内存需求≈单块大小×N。比如重叠3个栅格，就是4MB×3=12MB。
3. GDAL缓存设置：可以通过terraOptions()限制GDAL内存使用：

terraOptions(memfrac=0.2)  # 限制GDAL使用20%的可用内存

估算公式：最小内存需求 = 最大单块大小 × 重叠源栅格数，建议预留2-3倍冗余空间，避免内存碎片导致崩溃。

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