Fastest way to perform focal operations using R

Question

I know there are already some questions out there discussing this issue, for example:

How parallelize the extract function for raster files in R?

Increasing speed of crop, mask, & extract raster by many polygons in R?

But I'm not really sure if there is any "best practice" way to approach this. The setup of problem I describe is that I have one (or more, but the number of rasters does not matter too much I think) raster and a vector-object with many polygons. Now I want to get some aggregated value for each polygon. What I found so for is that the normal raster::extract has the advantage of potentially returning a sp-object by setting sp=TRUE, however, it is raster slow.

The exact_extract-function seems to be way faster, but when using an aggregation function I can not append the aggregated value as a column at the vector (lets say sf-dataframe).

My question now would be if there is any way of making this even faster. Maybe using some kind of parallel approach or rasterizing the shape as proposed in the second link.

Answer 1

Following my comment I will try to explain the steps of using tabularaster to speed up the extraction process.

First let's get ourselves a vector data object and a raster:

library(raster)
library(tabularaster)
library(sf)
library(dplyr)

## Get some polygons

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


# Generate a raster that covers the extent of the sf object
r <- raster()
extent(r) <- extent(nc)
res(r) <- 0.05
r[] <- runif(ncell(r))

plot(r) # preview raster

We now have a vector object representing 100 counties in North Carolina, and a random raster spanning that area. We use tabularaster to create an index linking each polygon (from nc object) to each cell (from r object).

index <- cellnumbers(r, nc)

head(index, 10) 

# A tibble: 10 x 2
   object_ cell_
     <dbl> <int>
 1       1    54
 2       1    55
 3       1    56
 4       1    57
 5       1    58
 6       1    59
 7       1    60
 8       2    61
 9       2    62
10       2    63

What does that mean? It means that the first polygon in nc overlaps cells number 54 to 60 in the raster. The second polygon overlaps cells 61, 62, 63 etc.

You can check how many cells are contained within each polygon:

index %>% group_by(object_) %>% count()

# A tibble: 100 x 2
# Groups:   object_ [100]
   object_     n
     <dbl> <int>
 1       1    49
 2       2    26
 3       3    57
 4       4    26
 5       5    59
 6       6    36
 7       7    22
 8       8    34
 9       9    48
10      10    50
# ... with 90 more rows

Polygon 1 (nc[1, ]) overlaps 49 cells. Polygon 2 (nc[2, ]) overlaps 26 cells, and so on.

But what you really want to do is perform an extraction summarising the values of those cells for each polygon. Now that could be a sum, a mean, a modal value... I'll calculate the mean as an example.

The trick is to make use of the y argument in raster::extract() : "a numeric vector representing cell numbers". This tells the function exactly where to look and speeds things up massively. We perform the extraction at each cell that overlaps a polygon, then agregate based on those polygons and calculate our final value (in our case a mean, but you can use any function you want).

result <- index %>% 
   mutate(pixelvalue = raster::extract(r, cell_)) %>% # extract value of the raster at each pixel specified
   group_by(object_) %>%  # group by the polygon index
   summarise(finalvalue = mean(pixelvalue, na.rm = TRUE))  # agregate values

head(result)  # your results will differ as the raster is generated with random values!

# A tibble: 6 x 2
  object_ finalvalue
    <dbl>      <dbl>
1       1      0.506
2       2      0.489
3       3      0.544
4       4      0.510
5       5      0.497
6       6      0.434

You then end up with a dataframe that has a column representing the polygon index (from the original nc dataset), and the second column is the value you're after. As long as you didn't make any change to your nc object in the mean time, it's a one-to-one match and you can merge the two:

nc[result$object_, "myvalue"] <- result$finalvalue

plot(nc["myvalue"])

Answer 2

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"])