Avoid a for loop in raster::extract(rst,shp)

Question

I am working on R to extract the mean and maximum value of a raster within a 3 meter buffer of some buildings.

For this, I have created a for loop that iterates iterates through each building to extract this two values. My current code looks as follows:

for (b in c(1:nrow(buildings_shp))){
    
    building <- buildings_shp[b,]
    
    buffered <- st_buffer(building, 3)
    
    raster_cropped <- crop(raster, extent(buffered))

    mean <- extract(depths_cropped, buffered, fun = mean, na.rm = TRUE)
    max <- extract(depths_cropped, buffered, fun = max, na.rm = TRUE)
    
    buildings_shp[b,"mean"] <- mean
    buildings_shp[b,"max"] <- max
    
  }

This loop, however, takes a considerable amount of time (~17 minutes for 1500 buildings), and the step that seems to take most time is the two extract lines. I would like to know if there are ways to speed up this process by:

a) avoiding the use of a loop - the reason for this loop is that I fear that if I use st_buffer on the entire dataset, then when buildings are closer than 3 meters I would generate overlapping geometries, which may cause an error. UPDATE - having tested this, I confirmed that the results produced are different when using overlapping geometries than when iterating through each one of them

b) parallelizing the for loop (i have tried the raster clustering feature, but it did not speed up the process, probably because it did not parallelize the loop itself but the extract function)

c) using other function than raster::extract. I have seen some posts recommending the velox package, but it seems like this package has been removed from CRAN.

Some dummy data (copied from the referenced question above)

library(raster)
library(sf)

raster <- raster(ncol=1000, nrow=1000, xmn=2001476, xmx=11519096, ymn=9087279, ymx=17080719)
raster []=rtruncnorm(n=ncell(raster ),a=0, b=10, mean=5, sd=2)
crs(raster ) <- "+proj=utm +zone=51 ellps=WGS84"
    
x1 <- runif(100,2001476,11519096)
y1 <- runif(100, 9087279,17080719)

buildings_shp <- st_buffer(st_sfc(st_point(c(x1[1],y1[1]), dim="XY"),crs=32651),200000)

Answer 1

I am not clear on what you are getting at in saying that the results differ between a loop and overlapping geometries. There is obviously something amiss there and I cannot reproduce it. By default, st_buffer does not dissolve buffer polygons like rgeos::gBuffer does. You do have to take care when operating on multipart geometries as you can inadvertently end up operating on the entire geometry at once (eg., taking the mean of all buffers and not individually).

Let's first add packages and create the example data. Please add the exactextractr package as it is essential here.

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

r <- raster(ncol=1000, nrow=1000, xmn=2001476, xmx=11519096, 
            ymn=9087279, ymx=17080719)
  r[] <- dnorm(runif(ncell(r)), 5, 3)
xy <- data.frame(runif(100,2001476,11519096),
                 runif(100, 9087279,17080719))
  xy <- xy %>% 
    sf::st_as_sf(coords = c(1,2),crs=32651)

Here we illustrate the difference between the polygons being dissolved or remaining independent in representing the original feature geometry.

First, with dissolving overlapping polygons

xyb <- st_buffer(xy, 200000) %>%
   st_union()
     plot(xyb)

Now, as a multipart feature geometry which, if I understand correctly is what you want to use. We also perform a sanity check to ensure that we have the correct number of feature buffers.

xyb <- st_buffer(xy, 200000)
  plot(xyb)
    nrow(xyb) == nrow(xy)

Now that we have the buffered points sorted, we can use the exact_extract function, from the exactextractr package, to return the max and mean. Please read the help document as, aside from speed, an advantage of the function is that it returns the fractional intersections of each cell. As such you can perform weighted means, sums, etc... or just filter out pixels with very small intersections. There are some weighted options in the statistics provided in the function or you can return a matrix of the "raw" results and operate on them yourself.

head( e <- exact_extract(r, xyb, c("mean", "max")) )
  nrow(xy)
  nrow(e)

Let's make sure that the results are the same when operating on the entire data or subsetting (as in a loop) and, remember that we are working with a normal distribution in the source raster.

p <- sample(1:nrow(xy),10)
e[p,]
exact_extract(r, xyb[p,], c("mean", "max")) 

Here is the raw extracted data. Each list element represents the extracted values of the given polygon, "value" column, and the fractional intersection of the associated pixel, "coverage_fraction" column.

  head( rd <- exact_extract(r, xyb) )

A quick example of doing something with the raw extracted values. We will still calculate max and mean but, will first remove any pixels with < 0.65 fractional intersections.

s <- function(x, p = 0.65) {
  idx <- which(x$coverage_fraction < p)
    if(length(idx) > 0) x <- x[-idx,]
    return(c(mean(x$value), max(x$value)))
  }
d <- lapply(rd, s)
  d <- as.data.frame(do.call(rbind, d))
    names(d) <- c("mean","max")
      head(d)

Join back to data

( xy <- cbind(xy, d) )

In the R spirit of "more than one way to skin a cat" here are two other approaches that collapse the entire workflow. Please note that the use of lapply is functionally exactly the same as a for loop.

do.call(rbind, lapply(1:nrow(xy), 
function(i) exact_extract(r, st_buffer(xy[i,],200000), c("mean", "max"))))
    
cbind(xy, exact_extract(r, st_buffer(xy, 200000), c("mean", "max")))

BTW, the system time for the last example with 1500 random observations was 1.61 seconds.