1

I'm wondering if I have maximized the speed at which a mean of an area buffered around a point in a raster can be extracted.

Fastest way I have ever found to extract a raster, with the pre-cropping suggestion by @dbaston:

If you have the velox raster already (even if you have to buffer the shape dynamically), this is lightning:

test7_LIGHTNING <- system.time(mclapply(seq_along(x1), function(x){
  q <- vras$crop(poly);vras$extract(poly, fun=function(t) mean(t,na.rm=T))
}))

> test7_LIGHTNING
      user  system elapsed 
     0.001   0.005   0.355 

These are still quick if you need to dynamically load velox raster (simulating loading different rasters across a set, but here loading the same raster repetitively from sample reproducible data below):

test8 <- system.time(mclapply(seq_along(x1), function(x){ ras<-velox("testras_so.tif");ras$crop(poly);ras$extract(poly, fun=function(t) mean(t,na.rm=T)) }))

test9 <- system.time(mclapply(seq_along(x1), function(x){ ras<-velox("testras_so.tif");ras$crop(st_buffer(st_sfc(st_point(c(x1[x],y1[x]), dim="XY"),crs=32651),200000));ras$extract(poly, fun=function(t) mean(t,na.rm=T)) }))


> test8
   user  system elapsed 
  0.011   0.016   4.450 
> test9
   user  system elapsed 
  0.006   0.012   4.333 

ORIGINAL QUESTION: Can performance be improved any further on these LOCALLY? I use parallel mclapply already, and I know I could get further gains by setting up and running this on a cluster (use a cluster or get more cpu's is not the answer I'm looking for).

Replicate some data:

library(raster)
library(parallel)
library(truncnorm)
library(gdalUtils)
library(velox)
library(sf)
ras <- raster(ncol=1000, nrow=1000, xmn=2001476, xmx=11519096, ymn=9087279, ymx=17080719)
ras[]=rtruncnorm(n=ncell(ras),a=0, b=10, mean=5, sd=2)
crs(ras) <- "+proj=utm +zone=51 ellps=WGS84"

writeRaster(ras,"testras_so.tif", format="GTiff")

gdalbuildvrt(gdalfile = "testras_so.tif", 
             output.vrt = "testvrt_so.vrt")

x1 <- runif(100,2001476,11519096)
y1 <- runif(100, 9087279,17080719)

poly <- st_buffer(st_sfc(st_point(c(x1[1],y1[1]), dim="XY"),crs=32651),200000)
vras <- velox("testvrt_so.vrt")
###########

Tests (test1: if have poly and velox raster, test2: if have to generate buffer but have velox raster, test3: if have to generate velox from VR (simulating having different rasters) but having the buffer, test4: have to generate both (from VR):test5: generate velox from tif have buffer, test6: generate both (tif version).

#Test time if have poly and velox raster
test1 <- system.time(mclapply(seq_along(x1), function(x){
  vras$extract(poly, fun=function(t) mean(t,na.rm=T))
}))

#Test time if have to generate buffer but have velox raster
test2 <- system.time(mclapply(seq_along(x1), function(x){
  vras$extract(st_buffer(st_sfc(st_point(c(x1[x],y1[x]), dim="XY"),crs=32651),200000), fun=function(t) mean(t,na.rm=T))
}))

#Test time if have to generate velox from VR (simulating having different rasters) but having the buffer
test3 <- system.time(mclapply(seq_along(x1), function(x){
  velox("testvrt_so.vrt")$extract(poly, fun=function(t) mean(t,na.rm=T))
}))

#Test time if have to generate velox from VR AND generate buffer (simulating a list of rasters with different buffers each)
test4 <- system.time(mclapply(seq_along(x1), function(x){
  velox("testvrt_so.vrt")$extract(st_buffer(st_sfc(st_point(c(x1[x],y1[x]), dim="XY"),crs=32651),200000), fun=function(t) mean(t,na.rm=T))
}))

#Test time if have to generate velox from TIF (simulating having different rasters) but having the buffer
test5 <- system.time(mclapply(seq_along(x1), function(x){
  velox("testras_so.tif")$extract(poly, fun=function(t) mean(t,na.rm=T))
}))

#Test time if have to generate velox from TIF AND generate buffer (simulating a list of rasters with different buffers each)
test6 <- system.time(mclapply(seq_along(x1), function(x){
  velox("testras_so.tif")$extract(st_buffer(st_sfc(st_point(c(x1[x],y1[x]), dim="XY"),crs=32651),200000), fun=function(t) mean(t,na.rm=T))
}))

My results (yours will vary with cores due to mclapply running parallel):

   > test1
   user  system elapsed 
  0.007   0.022   3.501 
> test2
   user  system elapsed 
  0.008   0.025   3.851 
> test3
   user  system elapsed 
  0.018   0.036  10.546 
> test4
   user  system elapsed 
  0.020   0.042  10.754 
> test5
   user  system elapsed 
  0.015   0.034   9.143 
> test6
   user  system elapsed 
  0.015   0.033   9.364

Can anybody make any suggestions to make this faster or have I maxed local speed here?

7

I get much faster results with velox if I crop the raster before running extract, e.g.:

r <- velox("testras_so.tif")
r$crop(poly)
r$extract(poly)

I've also been working on a package with an optimized extract function that may be of interest:

devtools::install_github('isciences/exactextractr')
library(exactextractr)

exact_extract(ras, poly)                # get a matrix with weights and values,
                                        # as with raster::extract(, weights=TRUE)
exact_extract(ras, poly, weighted.mean) # pass weights and values to an R function
exact_extract(ras, poly, 'mean')        # use a built-in stat

For this case, runtimes look comparable to velox but take into account partially covered cells (comparable to weights=TRUE with raster::extract).

library(microbenchmark)

v <- velox(ras)
microbenchmark(
  v$extract(poly),
  {q <- v$copy(); q$crop(poly); q$extract(poly)},
  exact_extract(ras, poly)
)

# Unit: milliseconds
#                     expr      min       lq     mean   median       uq       max neval
# v$extract(poly)          38.937180 39.749301 41.703686 40.247412 41.217893 56.372611   100
# { q <- v$copy()...        2.510894  2.639588  2.769510  2.725196  2.844109  4.850091   100
# exact_extract(ras, poly)  1.818466  2.039748  3.193071  2.255309  2.382927 35.236612   100
  • 1
    Why are you distributing this under an Apache License? You will never be able to distribute this on CRAN and it is not exactly a GNU open source. – Jeffrey Evans Jun 15 '18 at 15:38
  • Huh. Do you know why this is? Everything I've read tells me that velox isn't actually loading any values, just the metadata and based on that I shouldn't be getting these performance increases by pre-cropping, but you are right, tests show pre cropping is faster! – Neal Barsch Jun 15 '18 at 22:44
  • test6 <- system.time(mclapply(seq_along(x1), function(x){ ras<-velox("testras_so.tif");ras$crop(poly);ras$extract(poly, fun=function(t) mean(t,na.rm=T)) })) test7 <- system.time(mclapply(seq_along(x1), function(x){ velox("testras_so.tif")$extract(poly, fun=function(t) mean(t,na.rm=T)) })) test6(elapsed):4.400 test7(elapsed):7.791 – Neal Barsch Jun 15 '18 at 22:45
  • Also thanks for this, with that this is lightning in parallel, velox, I'm doing 30x extraction performance I was a couple years ago with sp and base raster packages – Neal Barsch Jun 15 '18 at 22:53
  • @JeffreyEvans there are a number of Apache-licensed packages on CRAN; see length(grep('Apache', available.packages()[, "License"])). I'm open to something else, but I'd need to understand what problems this license poses (it is more permissive than GPL and recommended by FSF for projects that don't want to use copyleft) – dbaston Jun 18 '18 at 13:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.