Return to Revisions

Let's benchmark some different approaches.

Here is some example data. Please note that we explode the MULTIPOLYGON geometry to simplify matters a bit.

library(sf)
library(raster)
library(terra)
library(exactextractr)

nc <- st_read(system.file("shape/nc.shp", package="sf"))
  nc <- st_cast(nc, "POLYGON")

r <- raster(extent(nc), res=0.05)
  r[] <- runif(ncell(r))
    proj4string(r) <- st_crs(nc)$proj4string

  plot(r)
    plot(nc[,1], add=TRUE)

raster::extract using lapply on resulting object and raster::extract using an internal function call

system.time({
  unlist(lapply(raster::extract(r, nc), mean))
})

system.time({
  raster::extract(r, nc, fun=mean)
})

terra::extract using tapply on resulting object and terra::extract using an internal function call. Note that terra is now returning a data.frame with a unique ID denoting the polygons (probably indicating rownames). As such, we use tapply rather than apply. This allows us to aggregate as stastic by a unique ID.

system.time({
  v.terra <- terra::extract( rast(r), vect(nc))
    tapply(v.terra[,"layer"], v.terra[,"ID"], mean) 
})

system.time({
  terra::extract( rast(r), vect(nc), fun=mean)
})

exactextractr::exact using lapply on resulting object and exactextractr::exact using an internal function call

system.time({
  v.exact <- exactextractr::exact_extract(r, nc) 
    unlist(lapply(v.exact, function(x) mean(x[,"value"]) ))
})

system.time({
  exactextractr::exact_extract(r, nc, "mean") 
})

All of these approaches result in a vector that can be joined back to the polygon data. You can see that the exactextractr::exact is the fastest approach. For both raster::extract and terra::extract you may observe that it can be faster to create an object and then use lapply, as opposed to calling a function internally. This is even more noticeable when it is a custom function (notably so). The exactextractr::exact function has several moment statistics available that are calculated in C++ rather than calling an R function. This speeds up things quite a bit compared to creating an object and using lapply however, this is not the case for custom function and I have noticed that even with this function with non-C++ function it is faster to create an object and use lapply rather than passing the function using the fun argument.

Here is a simple example of calculating the proportion of value >= to 0.25 and assigning results to the polygon data.

m <- function(x, p=0.25) {
  nrow(x[x[,1] >= p,]) / nrow(x) 
}
v.exact <- exactextractr::exact_extract(r, nc) 
  nc$p025 <- unlist(lapply(v.exact, m))

plot(nc[,"p025"])