1

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.

5
  • 2
    Have you tried the tabularaster package? It can be used to create a cell index before running the raster::extract function, and was a life saver for me. See vignette – Sandy AB Feb 3 at 16:20
  • Thanks a lot:)! Yesterday I saw it, but didn't know exactly what to do with it – Lenn Feb 4 at 7:58
  • Can you provide a bit of code to help me understand what you mean by "I can not append the aggregated value as a column at the vector" ? – dbaston Feb 4 at 15:39
  • Not sure what your exact issue is here. If you mean by "aggregate" that your values are returned as a list all you have to do is use lapply or pass the desired function to exact_extract eg., "mean" which is derived in C++ (see function help). This will FAR outperform the recommendation of raster::extract. – Jeffrey Evans Feb 4 at 19:21
  • 1
    Please clarify your language here. Your title is intimating that you are interested in focal functions (matrix operations) yet your question seems focused on extracting raster values associated with vector data. These are two very different things in raster algebra parlance. – Jeffrey Evans Feb 4 at 19:58
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"])

Hope this helps!

5
  • Why all the overhead when these functions can return several moments themselves? For example raster::extract(r, fun="mean", buffer=1000) will return the mean of all cells within a 1000m radius (does require a projected raster otherwise you need to figure out your buffer distance in DD. Another approach is that focal functions are fairly fast so, you can just calculate a statistic using a nXn window and then assign the straight pixel values. Your approach is the LONG away around and irrelevant for data extracted using polygon. There is also lapply. – Jeffrey Evans Feb 4 at 19:30
  • 1
    Thanks @JeffreyEvans - I'm sure you're right, as I am fairly novice at GIS in R, and perhaps my solution is not the most efficient. It did however do its job when I needed to extract raster values in a large polygon dataset (everything else I had tried was very slow), so I wouldn't say it is irrelevant. I'm not sure what you mean by recommending buffer in the extract function as I believe that is more the remit of focal statistics? (OP said focal but I believe meant zonal) I've learned a lot from your answers on here so will have to try your lapply solution. Cheers. – Sandy AB Feb 5 at 11:39
  • There are almost always a dozen ways to skin a cat in R so, whatever works for you is the right way. As you gain experience it becomes about code efficiency and reproducibility as well. The more verbose the code, the more can go wrong with a different problem, and the less reproducible the solution. Also, as you are using pipe approaches, start thinking about what is being done as each additional function call in the workflow. I would highly recommend reading the R Inferno (available free online) to start thinking about processing efficiency. – Jeffrey Evans Feb 5 at 15:13
  • Thank you both very much!! I don't have that much experience in R yet but really really like it. Sometimes it's just getting super confusing in comparison with what I learned in python for example. As @JeffreyEvans points out there are often mulitple (and it feels like way more as in other software) ways in R to achieve the same results. What can be an advantage might be a big problem for beginners and intermediates. Also in terms of code quality as something will most often end up working even though not knowing exactly why. Anyways, thanks a lot! – Lenn Feb 9 at 18:42
  • I though this would maybe be a process to be paralellized in order to be fast, but I think this is not the right approach to this issue – Lenn Feb 9 at 18:42
1

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"])
1
  • Wow, thanks a lot for the detailed answer! I wanted to take some minutes to read over it and so I just answered today... This definitely helped a lot in my understanding. Also of the basic processes behind raster::extract and similar functions – Lenn Feb 9 at 18:39

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.