Speed up raster::extract with weights in r

user_2257981 0
Viewed 812

I want to extract the precise mean value of raster values from an area extent defined by a polygon in r. This works using raster::extract with the option weights=TRUE. However, this operation becomes prohibitively slow with large rasters and the function doesn't seem to be parallelized, thus beginCluster() ... endCluster() does not speed up the process. I need to extract the values for a range of rasters, exemplified here as r, r10 and r100. Is there a way to speed this up in r, or is there an alternative way of doing this in 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

r1000 and the area to be extracted

raster extract vector parallel-processing r
2 Answers

In the end I resorted the problem using gdalUtils working directly on the harddisk.

I used the command gdalwarp() to reduce the raster resolution to r10, 100. Then gdalwarp() to increase the resolution of the resulting raster to the original resolution of r. Then gdalwarp() with cutline= "v.shp", crop_to_cutline =T to mask the raster to the vector v. And then gdalinfo() combined with subset(x(grep("Mean=",x))) to extract the mean values. All of this was packed in a foreach() %dopar% loop to process a number of rasters and resolution.

While complicated and probably not as precise as extract::raster, it did the job.

user_22579810

It should actually run faster if you first call beginCluster (the function then deals with the parallelization). Even better would be to use version 2.7-14 which has a much faster implementation. It is currently under review at CRAN, but you can also get it here: https://github.com/rspatial/raster

user_6352450
Related