Extracting and summarize raster for a large polygon set

Question

I have a large multipolygon (5x5km grids) file (1.23 million polygons) that I would like to exact raster data to each polygon. Specifically I want to calculate the fractal dimensions (a user defined function), Morans I (from the raster package), and the mean NDVI value of all the pixels that fall into each polygon.

The raster covers the extent of Australia. Each band has an associated date (919 total bands in the raster). Importantly, all polygons need only one specific band.

However, this loop takes such a long time (roughly 2 seconds per polygon). The crop command is what makes the loop run slowly. What is the best way to speed this loop up?

Below is the basic for loop

#adding column for final values to be returned
#spdf3_buffed is a sf object

spdf3_buffed$FractID <- 0  
spdf3_buffed$meanNDVI$Moran <- 0
spdf3_buffed$meanNDVI <- 0

for (i in 1:nrow(spdf3_buffed)){
  print(i) 
  spSubset <- spdf3_buffed[i,]
  BandNumber <-spSubset$BandNumber    
  r <- raster(MODIS_raster_location,band=BandNumber)
  cropped_raster <- crop(r, st_bbox(spSubset))
  b <- as.numeric(fractalD(cropped_raster))
  c <- as.numeric(Moran(cropped_raster))
  layermeans <- cellStats(cropped_raster, stat='mean', na.rm=TRUE)
  u <- mean(layermeans)
  spdf3_buffed$FractID[[i]] <- b
  spdf3_buffed$Moran[[i]] <- c
  spdf3_buffed$meanNDVI[[i]] <- u
}

write.csv(spdf3_buffed,file="test.csv")

Answer 1

It seems that the slowest part of the code is the call to Moran. I am using raster::extent in lue of sf::st_bbox which seems to make a difference. I also moved the for loop into an lapply.

Let's add libraries and create some example data. I would recommend reading the data into R as a stack, and index out of the stack, rather than repeatedly reading a raster at each iteration. I am using the nc data that comes with sf and adding a "BandNumber" to emulate your code. This is what is being indexed from the stack using double brackets.

library(sf)
library(raster)

nc <- st_cast(st_read(system.file("shape/nc.shp", 
              package="sf")), "POLYGON")
 
i=500; j=500
r <- do.call(raster::stack, replicate(20, 
             raster::raster(matrix(sample(c(0,1), i*j, replace=TRUE), i, j)))) 
    extent(r) <- extent(nc)
      proj4string(r) <- st_crs(nc)$proj4string

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

nc$BandNumber <- sample(1:nlayers(r), nrow(nc), replace=TRUE)

Here is our call to lapply for implementing a loop. It looked like you wanted an echo so, I added one using cat, Sys.sleep and flush.console which is a trick to get the message to actually print. Since you do not tell us where fractalD comes from, I am using the lsm_c_frac_mn function from landscapemetrics. I would not use cellStats when you can simply pass a vector of raster values to mean with na.rm=TRUE.

s <- lapply(1:nrow(nc), function(i) {
  cat("processing", i, "of", nrow(nc), "\n")
    Sys.sleep(0.01)
      flush.console()
  rsub <- crop(r[[nc[i,]$BandNumber]], extent(nc[i,]))
    fd <- as.data.frame(landscapemetrics::lsm_c_frac_mn(rsub))[,6][2]
      c(fd, as.numeric(Moran(rsub)), mean(rsub[], na.rm=TRUE)) } )

The results are a list object so, we can collapse them into a data.frame and join back to the polygon data.

( s <- as.data.frame(do.call("rbind", s)) )
  names(s) <- c("FractID", "Moran", "meanNDVI")  
   s$CNTY_ID <- nc$CNTY_ID        
     nc <- dplyr::left_join(nc, s)

With over a million observations, this is still going to take quite some time. I would consider breaking the problem into multiple pieces. You could also try to multi-thread the problem but, breaking your polygon data into multiple pieces and running each subset in a difference instance of R will speed this up considerably, especially if you have the cores to also support multi-threading using something like future.apply::future_lapply (parallel version of lapply).