在大型栅格上进行算术运算的最快方法

Question

我有大量大型栅格（全局范围，250 米分辨率；大约 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 支持...

在大型栅格上进行算术运算的最快方法

fastest way to conduct arithmetic on large rasters

r

raster

r-raster

一些示例方法：

一些虚拟数据以使其可重现