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

Posted 2023-04-18

【中文标题】对大型栅格进行算术运算的最快方法

【英文标题】：fastest way to conduct arithmetic on large rasters

【发布时间】：2020-07-10 15:56:31

【问题描述】：

我有大量大型栅格（全局范围，250 m 分辨率；大约 1e10 个浮点单元格）- 文件名在矢量 deltaX.files 中。我想将这些中的每一个添加到另一个栅格，文件名X.tif。由于此操作可能需要几天时间才能完成，我想知道添加栅格以使其尽可能快的最快方法。

我可以想到一些方法，但我不确定哪种方法最有效，或者是否还有比这些方法更好的选择。

所以，我的问题是是否有办法优化或显着加快大型栅格的算术运算。请注意，我有一个启用 CUDA 的 NVidia GPU，因此非常欢迎可以在 GPU 上并行处理的解决方案。 请注意，我使用的是 Linux 系统。

---

一些示例方法：

注意下面要插入的代码块，以确定默认的输出文件压缩、内存分配和启动并行集群

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

【参考方案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 支持...

【讨论】：

谢谢，我正在试用y <- sum(x... 版本。似乎比光栅包快几倍。但我找不到输出文件。不在当前工作目录中，搜索文件系统一无所获。我使用的确切命令是y <- sum(x, filename = new.f, na.rm=TRUE, datatype = "FLT4S", wopt = list(gdal = c("COMPRESS=ZSTD", "PREDICTOR=2")))

非常抱歉！ sum 没有文件名参数。我需要看看我是否可以解决这个问题；我建议上面的解决方法。

优秀。它现在工作（使用app）。我注意到它只使用一个核心。 terra 中是否有任何并行支持（例如栅格中的 beginCluster）？否则，我可以同时进行几个 R 会话，并在每个会话中并行处理不同的文件子集。

app 有一个 nodes 参数，但这仅适用于您使用自己的函数（而不是 sum、mean、min ...）时。在这种情况下，我会进行手动并行化（或并行化 for 循环----但是您必须在该循环中创建所有 SpatRaster 对象（它们不会序列化，您不能将它们发送到节点，它们必须在那里创建） )。