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

Question

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

3

我正在处理大型光栅堆栈，并需要重新采样和裁剪它们。我读取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

- Geo-sp

3

如果没有手头上的栅格文件，很难评论可能导致这种情况的原因。尽管如此，我有时会使用由gdalUtils包装的GDAL函数加快栅格操作速度。在这里，我可能会使用gdal_translate()，通过其-tr参数设置分辨率，并通过其-r参数设置所需的重采样算法。不确定它如何处理RasterStack，但应该可以很好地处理RasterLayer（或猜测是磁盘上的*.tif*文件）。 - Josh O'Brien

1

请确保您使用最新版本的光栅化程序，因为“重新采样”速度最近有了很大的提升（可能还不够）。如果您展示show(s)和r，我们可以看到发生了什么，这也会非常有帮助。对于多核心计算，您可以尝试beginCluster等函数。 - Robert Hijmans

@JoshO'Brien 谢谢，我会尝试使用gdal，但raster会更方便。 - Geo-sp

@RobertH 我正在使用栅格版本2.4-20。我已经添加了r和s图层的信息。谢谢！ - Geo-sp

2个回答

1

对比一下，这是我得到的内容：

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

- Robert Hijmans

令人惊讶的是，在我的机器上速度较慢。 用户系统经过的时间 20.45 5.22 25.75 - Geo-sp

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jbaums · Accepted Answer

我支持@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)

使用4个并行进程，在10个网格的情况下需要32秒，对于365个文件来说大约需要20分钟。实际上，速度应该会更快，因为在最后两个线程可能没有做任何事情（10不是4的倍数）。