r - 对大型栅格进行算术运算的最快方法

Question

我有大量的大型栅格(全局范围，250 m 分辨率；大约 1e10 个浮点单元格)--文件名在向量中 deltaX.files .我想将其中的每一个添加到另一个光栅，文件名 X.tif .由于此操作可能需要数天才能完成，我想知道添加栅格以使其尽可能快的最快方法是什么。
我可以想到一些方法，但我不确定哪个是最有效的选择，或者是否有比其中任何一个更好的选择。
所以，我的问题是是否有一种方法可以优化或显着加快大型栅格的运算速度。请注意，我有一个支持 CUDA 的 NVidia GPU，因此非常欢迎可以在 GPU 上并行化的解决方案。 请注意，我使用的是 Linux ystsem。

一些示例方法:
注意下面的代码块要插入在每个之前，以确定默认输出文件压缩、内存分配和启动并行集群

rasterOptions(chunksize = 1e10, maxmemory = 4e10)
f.opt = '-co COMPRESS=ZSTD -co PREDICTOR=2'
f.type = 'FLT4S'
beginCluster()

选项1)

for (f in deltaX.files) {
  s = stack('X.tif', f)
  calc(s, sum, filename = paste0('new_', f), datatype = f.type, options = f.opt)
}

选项 (2)

X = raster('X.tif')
for (f in deltaX.files) {
  dX = raster(f)
  overlay(X, dX, fun=sum, filename = paste0('new_', f), datatype = f.type, options = f.opt)
}

选项 (3)

X = raster('X.tif')
for (f in deltaX.files) {
  dX = raster(f)
  Y = X + dX
  writeRaster(Y, filename = paste0('new_', f), datatype = f.type, options = f.opt)
}

选项(4):使用 gdal_calc.py 而不是 R

for (f in deltaX.files) {
  system(cmd)
  cmd = paste0("gdal_calc.py -A X.tif ", "-B ", f, " --outfile=", 'temp.tif', ' --calc="A+B"')
  system(cmd)
  system(paste('gdal_translate -ot Float32', f.opt, 'temp.tif', paste0('new_', f)))
  system('rm temp.tif')
}

请注意，我在成功生成完全压缩的输出文件时遇到了麻烦，因此还需要在每个文件上使用 gdal_translate 进行压缩的附加步骤。然而，在一些测试运行中，它似乎产生了损坏的值，所以我真的对 R 解决方案最感兴趣，而不是使用 gdal_calc.py .

一些虚拟数据使这个可重现

X = raster(vals = rnorm(65000 * 160000), ncol = 160000, nrow = 65000)
writeRaster(X, 'X.tif',  datatype = f.type, options = f.opt)
for (i in 1:10) {
  dX = raster(vals = rnorm(65000 * 160000), ncol = 160000, nrow = 65000)
  writeRaster(X, paste0('dX', i, '.tif'),  datatype = f.type, options = f.opt)
}
deltaX.files = paste0('dX', 1:10, '.tif')

Answer 1

我建议使用 terra (一个旨在取代 raster 的新包——它更简单、更快)。它现在可以从 CRAN 获得，但对于尖端技术，您可以 install from github
可能最好的方法是

library(terra)
r <- rast(c('X.tif')
for (f in deltaX.files) {
    s <- rast(f)
    x <- c(r, s)
    y <- app(x, sum, filename=paste0('new_', f), datatype="INT2S",
           wopt=list(gdal="COMPRESS=LZW") )
 }

也许下面的要快一点；但问题是它没有文件名参数。但你可以解决这个问题

library(terra)
r <- rast(c('X.tif')
for (f in deltaX.files) {
    s <- rast(f)
    x <- r + s
    tempfile <- sources(x)$source[1]
    file.rename(tempfile, paste0('new_', f))
 }

或者，在一个步骤中(这会创建一个巨大的文件——可能不需要):

r <- rast(c('X.tif')
s <- rast(deltaX.files)
# combine them as separate sub-datasets
x <- sds(r, s)
y <- sum(x, filename="file.tif")

或者像这样(很快，但它会转到临时文件，完成后可以重命名，但不能设置所有写入选项)

 z <- r + s 

尚无GPU支持...