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

Question

我正在处理大型栅格堆栈，因此需要重新采样和裁剪。我阅读了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中的过程如此缓慢，以及是否有一种加快过程的方法？（我使用了snow软件包进行并行处理，但这也无济于事）。这些是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 

Answer 1

我赞同@ JoshO'Brien的建议，直接使用GDAL，gdalUtils这很简单。

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

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个网格的32秒（使用4个并发进程）下，您需要大约20分钟才能处理365个文件。实际上，它应该比那快，因为两个线程最后可能什么都不做（10不是4的倍数）。