使用 r 中的权重加速 raster::extract

Speed up raster::extract with weights in r

我想从 r 中的多边形定义的区域范围中提取栅格值的精确平均值。这可以使用 raster::extract 和选项 weights=TRUE。但是,对于大型栅格,此操作变得非常慢,并且该函数似乎没有并行化,因此 beginCluster() ... endCluster() 不会加快该过程。 我需要提取一系列栅格的值,此处示例为 r、r10 和 r100。有没有办法在 r 中加快速度,或者在 GDAL 中有其他方法吗?

r <- raster(nrow=1000, ncol=1000, vals=sample(seq(0,0.8,0.01),1000000,replace=TRUE))
r10 <- aggregate(r, fact=10)
r100 <- aggregate(r, fact=100)

v = Polygons(list(Polygon(cbind(c(-100,100,80,-120), c(-70,0,70,0)))), ID = "a")
v = SpatialPolygons(list(v))

plot(r)
plot(r10)
plot(r100)
plot(v, add=T)

system.time({
precise.mean <- raster::extract(r100, v, method="simple",weights=T, normalizeWeights=T, fun=mean)    
})
user  system elapsed 
0.251   0.000   0.253 
> precise.mean
      [,1]
[1,] 0.3994278    


system.time({
  precise.mean <- raster::extract(r10, v, method="simple",weights=T, normalizeWeights=T, fun=mean)    
})

user  system elapsed 
7.447   0.000   7.446 

precise.mean
      [,1]
[1,] 0.3995429

如果您首先调用 beginCluster(该函数然后处理并行化),它实际上应该 运行 更快。更好的方法是使用 2.7-14 版本,它的实现速度要快得多。它目前正在 CRAN 审查中,但您也可以在这里获得它:https://github.com/rspatial/raster

最后我使用 gdalUtils 直接在硬盘上工作解决了这个问题。

我使用命令 gdalwarp() 将光栅分辨率降低到 r10, 100。 然后 gdalwarp() 将生成的光栅的分辨率提高到 r 的原始分辨率。 然后 gdalwarp()cutline= "v.shp", crop_to_cutline =T 将栅格屏蔽到矢量 v。 然后 gdalinfo() 结合 subset(x(grep("Mean=",x))) 提取平均值。 所有这些都打包在一个 foreach() %dopar% 循环中以处理大量光栅和分辨率。

虽然复杂而且可能不如 extract::raster 精确,但它完成了工作。