如何提高处理大型栅格堆栈的 R 处理速度？答案

【问题标题】：How to increase R processing speed dealing with large raster stacks?如何提高处理大型栅格堆栈的 R 处理速度？
【发布时间】：2015-11-09 20:21:06
【问题描述】：

我正在处理大型光栅堆栈，我需要重新采样和剪辑它们。我阅读了 Tiff 文件列表并创建了堆栈：

files <- list.files(path=".", pattern="tif", all.files=FALSE, full.names=TRUE) 
s <- stack(files)
r <- raster("raster.tif")
s_re <- resample(s, r,method='bilinear')
e <- extent(-180, 180, -60, 90)
s_crop <- crop(s_re, e)

这个过程需要几天才能完成！但是，使用 ArcPy 和 python 会快得多。我的问题是：为什么 R 中的过程如此缓慢，是否有办法加快该过程？（我使用雪包进行并行处理，但这也无济于事）。这些是r 和s 层：

> r
class       : RasterLayer 
dimensions  : 3000, 7200, 21600000  (nrow, ncol, ncell)
resolution  : 0.05, 0.05  (x, y)
extent      : -180, 180, -59.99999, 90.00001  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

> s
class       : RasterStack 
dimensions  : 2160, 4320, 9331200, 365  (nrow, ncol, ncell, nlayers)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

【问题讨论】：

如果没有您的光栅文件，很难评论可能导致此问题的原因。也就是说，我偶尔会通过使用gdalUtils 包装的 GDAL 函数来加速光栅操作。在这里，我可能会使用gdal_translate()，通过其-tr 参数设置分辨率，并通过其-r 参数设置所需的重采样算法。不确定它如何处理RasterStacks，但它应该可以处理RasterLayers（或者，我猜，*.tif* 磁盘上的文件）就好了。
确保您使用最新版本的栅格，因为resample 的速度最近已经提高了很多（可能还不够）。如果您显示 show(s) 和 r，我也会很有帮助，以便我们可以看到发生了什么。对于多核，您可以尝试 beginCluster 等
@JoshO'Brien 谢谢，我会尝试 gdal，但光栅会更方便。
@RobertH 我使用的是光栅版本 2.4-20。我添加了 r 和 s 层的信息。谢谢！

标签： r parallel-processing raster snow

【解决方案1】：

我赞同 @JoshO'Brien 直接使用 GDAL 的建议，gdalUtils 让这一点变得简单。

这是一个使用与您的尺寸相同的双精度网格的示例。对于 10 个文件，在我的系统上大约需要 55 秒。它是线性缩放的，因此 365 个文件大约需要 33 分钟。

library(gdalUtils)
library(raster)

# Create 10 rasters with random values, and write to temp files
ff <- replicate(10, tempfile(fileext='.tif'))
writeRaster(stack(replicate(10, {
  raster(matrix(runif(2160*4320), 2160), 
         xmn=-180, xmx=180, ymn=-90, ymx=90) 
})), ff, bylayer=TRUE)

# Define clipping extent and res
e <- bbox(extent(-180, 180, -60, 90))
tr <- c(0.05, 0.05)

# Use gdalwarp to resample and clip 
# Here, ff is my vector of tempfiles, but you'll use your `files` vector
# The clipped files are written to the same file names as your `files`
#  vector, but with the suffix "_clipped". Change as desired.
system.time(sapply(ff, function(f) {
  gdalwarp(f, sub('\\.tif', '_clipped.tif', f), tr=tr, r='bilinear', 
           te=c(e), multi=TRUE)
}))

##    user  system elapsed 
##    0.34    0.64   54.45

您可以进一步并行化，例如，parLapply:

library(parallel)
cl <- makeCluster(4)
clusterEvalQ(cl, library(gdalUtils))
clusterExport(cl, c('tr', 'e'))

system.time(parLapply(cl, ff, function(f) {
  gdalwarp(f, sub('\\.tif', '_clipped.tif', f), tr=tr, 
           r='bilinear', te=c(e), multi=TRUE)
}))

##    user  system elapsed 
##    0.05    0.00   31.92

stopCluster(cl)

对于 10 个网格（使用 4 个同时处理），时间为 32 秒，365 个文件大约需要 20 分钟。实际上，它应该比这更快，因为两个线程最后可能什么都不做（10 不是 4 的倍数）。

【讨论】：

谢谢！这工作正常，但如果我用它来读取所有文件，我会得到内存错误；而使用光栅你没有这个问题。你有什么建议吗？
parLapply/sapply 这样做。
@DNM - 我无法想象为什么......在我的系统上，GDAL 以块的形式处理网格。当我运行上述任何一个时，我的 RAM 使用量几乎没有变化。每个 GDAL 进程对我来说似乎使用了大约 100-200MB 的 RAM（因此具有 4 个进程的并行版本使用可能高达大约 1GB）。您的可用 RAM 已经很少了吗？
这很奇怪！现在我有 26 GB 的可用内存。
Error: cannot allocate vector of size 12.7 Gb

【解决方案2】：

作为比较，这是我得到的：

library(raster)
r <- raster(nrow=3000, ncol=7200, ymn=-60, ymx=90) 
s <- raster(nrow=2160, ncol=4320) 
values(s) <- 1:ncell(s)
s <- writeRaster(s, 'test.tif')

x <- system.time(resample(s, r, method='bilinear'))
#   user  system elapsed 
#  15.26    2.56   17.83

10 个文件需要 150 秒。所以我预计 365 个文件需要 1.5 小时 --- 但我没有尝试。

【讨论】：

令人惊讶的是它在我的机器上速度较慢。 user system elapsed 20.45 5.22 25.75