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

我正在处理大型光栅堆栈，我需要重新采样并剪辑它们。我读的TIFF文件列表，并创建堆栈：如何提高处理大型栅格堆栈的R处理速度？

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

来源

2015-11-09 Geo-sp

很难评论什么可能导致这种情况，而无需掌握光栅文件。也就是说，我偶尔使用由[** gdalUtils **]封装的GDAL函数来加速栅格操作（https://cran.r-project.org/web/packages/gdalUtils/index.html ）。在这里，我可能使用['gdal_translate（）']（http://www.gdal.org/gdal_translate.html），通过它的'-tr'参数设置分辨率，并通过'-r'参数设置所需的重采样算法。不知道它如何处理'RasterStack'，但它应该处理'RasterLayer's（或者，我猜，'* .tif *'磁盘上的文件）就好了。 –

确保您使用最新版本的光栅，因为'resample'的速度最近有了很大提高（可能还不够）。如果您显示节目和r以便我们可以看到正在发生的事情，我也会很有帮助。对于多核，你可以尝试beginCluster等 – RobertH

@ JoshO'Brien谢谢，我会尝试gdal，但光栅会更方便。 –

我第二@ 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

可以进一步与parallelise，例如，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)

在32秒，持续10个网格（使用4个同步过程），你正在寻找在约20分钟的时间365页的文件。实际上，它应该比这更快，因为两个线程最终可能什么都不做（10不是4的倍数）。

来源

2015-11-11 00:21:34 jbaums

谢谢！这工作正常，但如果我用它来读取所有文件，我得到内存错误;而光栅，你没有这个问题。你有什么建议吗？ –

@DNM哪部分导致内存错误？ – jbaums

parLapply/sapply会这样做。 –

为了便于比较，这是我得到：

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小时 - 但我没有尝试。

来源

2015-11-11 02:04:40 RobertH

令人惊讶的是它在我的机器上变慢了。 '用户系统已用完 20.45 5.22 25.75' –

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

回答

相关问题