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

Question

I'm extracting the area and percent cover of different land use types from a raster based on several thousand polygon boundaries. I've found that the extract function works much faster if I iterate through each individual polygon and crop then mask the raster down to the size of the particular polygon. Nonetheless, it's pretty slow, and I'm wondering if anyone has any suggestions for improving the efficiency and speed of my code.

The only thing I've found related to this is this response by Roger Bivand who suggested using GDAL.open() and GDAL.close() as well as getRasterTable() and getRasterData(). I looked into those, but have had trouble with gdal in the past and don't know it well enough to know how to implement it.

Reproducible Example:

library(maptools)  ## For wrld_simpl
library(raster)

## Example SpatialPolygonsDataFrame
data(wrld_simpl) #polygon of world countries
bound <- wrld_simpl[1:25,] #name it this to subset to 25 countries and because my loop is set up with that variable  

## Example RasterLayer
c <- raster(nrow=2e3, ncol=2e3, crs=proj4string(wrld_simpl), xmn=-180, xmx=180, ymn=-90, ymx=90)
c[] <- 1:length(c)

#plot, so you can see it
plot(c)    
plot(bound, add=TRUE) 

Fastest Method so far

result <- data.frame() #empty result dataframe 

system.time(
     for (i in 1:nrow(bound)) { #this is the number of polygons to iterate through
      single <- bound[i,] #selects a single polygon
      clip1 <- crop(c, extent(single)) #crops the raster to the extent of the polygon, I do this first because it speeds the mask up
      clip2 <- mask(clip1,single) #crops the raster to the polygon boundary

      ext<-extract(clip2,single) #extracts data from the raster based on the polygon bound
      tab<-lapply(ext,table) #makes a table of the extract output
      s<-sum(tab[[1]])  #sums the table for percentage calculation
      mat<- as.data.frame(tab) 
      mat2<- as.data.frame(tab[[1]]/s) #calculates percent
      final<-cbind(single@data$NAME,mat,mat2$Freq) #combines into single dataframe
      result<-rbind(final,result)
      })

   user  system elapsed 
  39.39    0.11   39.52 

Parallel Processing

Parallel processing cut the user time by half, but negated the benefit by doubling the system time. Raster uses this for the extract function, but unfortunately not for the crop or mask function. Unfortunately, this leaves a slighly larger amount of total elapsed time due to "waiting around" by the "IO."

beginCluster( detectCores() -1) #use all but one core

run code on multiple cores:

  user  system elapsed 
  23.31    0.68   42.01 

then end the cluster

endCluster()

Slow Method: The alternative method of doing an extract directly from the raster function takes a lot lot longer, and I'm not sure about the data management to get it into the form I want:

system.time(ext<-extract(c,bound))
   user  system elapsed 
1170.64   14.41 1186.14 

Answer 1

I have finally gotten around to improving this function. I found that for my purposes, it was fastest to rasterize() the polygon first and use getValues() instead of extract(). The rasterizing isn't much faster than the original code for tabulating raster values in small polygons, but it shines when it came to large polygon areas that had large rasters to be cropped and the values extracted. I also found getValues() was much faster than the extract() function.

I also figured out the multi-core processing using foreach().

I hope this is useful for other people who want an R solution for extracting raster values from a polygon overlay. This is similar to the "Tabulate Intersection" of ArcGIS, which did not work well for me, returning empty outputs after hours of processing, like this user.

#initiate multicore cluster and load packages
library(foreach)
library(doParallel)
library(tcltk)
library(sp)
library(raster)

cores<- 7
cl <- makeCluster(cores, output="") #output should make it spit errors
registerDoParallel(cl)

Here's the function:

multicore.tabulate.intersect<- function(cores, polygonlist, rasterlayer){ 
  foreach(i=1:cores, .packages= c("raster","tcltk","foreach"), .combine = rbind) %dopar% {

    mypb <- tkProgressBar(title = "R progress bar", label = "", min = 0, max = length(polygonlist[[i]]), initial = 0, width = 300) 

    foreach(j = 1:length(polygonlist[[i]]), .combine = rbind) %do% {
      final<-data.frame()
      tryCatch({ #not sure if this is necessary now that I'm using foreach, but it is useful for loops.

        single <- polygonlist[[i]][j,] #pull out individual polygon to be tabulated

        dir.create (file.path("c:/rtemp",i,j,single@data$OWNER), showWarnings = FALSE) #creates unique filepath for temp directory
        rasterOptions(tmpdir=file.path("c:/rtemp",i,j, single@data$OWNER))  #sets temp directory - this is important b/c it can fill up a hard drive if you're doing a lot of polygons

        clip1 <- crop(rasterlayer, extent(single)) #crop to extent of polygon
        clip2 <- rasterize(single, clip1, mask=TRUE) #crops to polygon edge & converts to raster
        ext <- getValues(clip2) #much faster than extract
        tab<-table(ext) #tabulates the values of the raster in the polygon

        mat<- as.data.frame(tab)
        final<-cbind(single@data$OWNER,mat) #combines it with the name of the polygon
        unlink(file.path("c:/rtemp",i,j,single@data$OWNER), recursive = TRUE,force = TRUE) #delete temporary files
        setTkProgressBar(mypb, j, title = "number complete", label = j)

      }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) #trycatch error so it doesn't kill the loop

      return(final)
    }  
    #close(mypb) #not sure why but closing the pb while operating causes it to return an empty final dataset... dunno why. 
  }
}

So to use it, adjust the single@data$OWNER to fit with the column name of your identifying polygon (guess that could have been built into the function...) and put in:

myoutput <- multicore.tabulate.intersect(cores, polygonlist, rasterlayer)

Answer 2

This question is very old so my answer is just for posterities. Since 2020-11-20 there is the new exactextractr package which implement the exact_extract() function, which is the new benchmark for raster extraction. The function take a Raster*object and an sf object as input and return the same output of raster::extract() plus the fraction of pixels covered by polygons. As raster::extract(), exact_extract() can take a fun argument to summarize data within polygons. Taking your original example (extracting area and percent cover of land uses within polygons), I've made a timing benchmark:

#generate data
data(wrld_simpl) #polygon of world countries
bound <- wrld_simpl[1:25,] #name it this to subset to 25 countries and because my loop 
is set up with that variable  

## Example RasterLayer
c <- raster(nrow=2e3, ncol=2e3, crs=proj4string(wrld_simpl), xmn=-180, xmx=180, 
ymn=-90, ymx=90)
c[] <- sample(1:8,length(c),replace = T)#I prefer to use a sample with replacement to better simulate land uses analysis

#plot, so you can see it
plot(c)    
plot(bound, add=TRUE) 

On my machin the original for loop took 18 sec to run

result <- data.frame() #empty result dataframe 

system.time(
for (i in 1:nrow(bound)) { #this is the number of polygons to iterate through
single <- bound[i,] #selects a single polygon
clip1 <- crop(c, extent(single)) #crops the raster to the extent of the polygon, I do 
this first because it speeds the mask up
clip2 <- mask(clip1,single) #crops the raster to the polygon boundary

ext<-extract(clip2,single) #extracts data from the raster based on the polygon bound
tab<-lapply(ext,table) #makes a table of the extract output
s<-sum(tab[[1]])  #sums the table for percentage calculation
mat<- as.data.frame(tab) 
mat2<- as.data.frame(tab[[1]]/s) #calculates percent
final<-cbind(single@data$NAME,mat,mat2$Freq) #combines into single dataframe
result<-rbind(final,result)
})#18.61 seconds

just by replacing extract() with exact_extract() (and doing some minor adjustment to data format) the loop run in 2 seconds!

library(exactextractr)
library(sf)
result <- data.frame() #empty result dataframe 

system.time(
for (i in 1:nrow(bound)) { #this is the number of polygons to iterate through
single <- bound[i,] #selects a single polygon
single_sf <- st_as_sf(single)
clip1 <- crop(c, extent(single)) #crops the raster to the extent of the polygon, I do 
this first because it speeds the mask up
clip2 <- mask(clip1,single) #crops the raster to the polygon boundary

ext<-as.data.frame(exact_extract(clip2,single_sf)) #extracts data from the raster 
based on the polygon bound
ext <- ext[complete.cases(ext),]
ext <- list(ext$value)
tab<-lapply(ext,table) #makes a table of the extract output
s<-sum(tab[[1]])  #sums the table for percentage calculation
mat<- as.data.frame(tab) 
mat2<- as.data.frame(tab[[1]]/s) #calculates percent
final<-cbind(single@data$NAME,mat,mat2$Freq) #combines into single dataframe
result<-rbind(final,result)
})#2.09

But looping is not the best way to deal with this problem. It is way better to do the extraction within all polygons and then manipulate the output to obtain the requested output (which is more concise and elegant):

library(tidyverse)
o <- Sys.time()
b <- st_as_sf(bound)
ex <- exact_extract(c,b,include_cols="NAME")
ex <- do.call(rbind,ex)
result <- ex %>% group_by(NAME) %>% summarise(
soil.use = sort(unique(value)),
count = table(value),
freq = table(value) / length(value)
)
Sys.time()-o #0.74 seconds

To generalize the problem just to the extraction routine, the extraction time only takes 0.11 seconds on my machine !!! This is a huge improvement.

For better details on the function and the package see https://github.com/isciences/exactextractr and ?exact_extract

EDIT

Since version 0.9 of exactextractr, the function exact_extract can take a spatRaster object as an argument, making the extraction even faster

Answer 3

Speed up extracting raster (raster stack) from point, XY or Polygon

Great answer Luke. You must be a R wizard! Here is a very minor tweak to simplify your code (may improve performance slightly in some cases). You can avoid some operations by using cellFromPolygon (or cellFromXY for points) and then clip and getValues.

Extract polygon or points data from raster stacks ------------------------

 library(raster)  
 library(sp)   

  # create polygon for extraction
  xys= c(76.27797,28.39791,
        76.30543,28.39761,
        76.30548,28.40236,
        76.27668,28.40489)
  pt <- matrix(xys, ncol=2, byrow=TRUE)
  pt <- SpatialPolygons(list(Polygons(list(Polygon(pt)), ID="a")));
  proj4string(pt) <-"+proj=longlat +datum=WGS84 +ellps=WGS84"
  pt <- spTransform(pt, CRS("+proj=sinu +a=6371007.181 +b=6371007.181 +units=m"))
  ## Create a matrix with random data & use image()
  xy <- matrix(rnorm(4448*4448),4448,4448)
  plot(xy)

  # Turn the matrix into a raster
  NDVI_stack_h24v06 <- raster(xy)
  # Give it lat/lon coords for 36-37°E, 3-2°S
  extent(NDVI_stack_h24v06) <- c(6671703,7783703,2223852,3335852)
  # ... and assign a projection
  projection(NDVI_stack_h24v06) <- CRS("+proj=sinu +a=6371007.181 +b=6371007.181 +units=m")
  plot(NDVI_stack_h24v06)
  # create a stack of the same raster
  NDVI_stack_h24v06 = stack( mget( rep( "NDVI_stack_h24v06" , 500 ) ) )


  # Run functions on list of points
  registerDoParallel(16)
  ptm <- proc.time()
  # grab cell number
  cell = cellFromPolygon(NDVI_stack_h24v06, pt, weights=FALSE)
  # create a raster with only those cells
  r = rasterFromCells(NDVI_stack_h24v06, cell[[1]],values=F)
  result = foreach(i = 1:dim(NDVI_stack_h24v06)[3],.packages='raster',.combine=rbind,.inorder=T) %dopar% {
     #get value and store
     getValues(crop(NDVI_stack_h24v06[[i]],r))
  }
  proc.time() - ptm
  endCluster()

user system elapsed 16.682 2.610 2.530

  registerDoParallel(16)
  ptm <- proc.time()
  result = foreach(i = 1:dim(NDVI_stack_h24v06)[3],.packages='raster',.inorder=T,.combine=rbind) %dopar% {
        clip1 <- crop(NDVI_stack_h24v06[[i]], extent(pt)) #crop to extent of polygon
        clip2 <- rasterize(pt, clip1, mask=TRUE) #crops to polygon edge & converts to raster
         getValues(clip2) #much faster than extract
  }
  proc.time() - ptm
  endCluster()

user system elapsed 33.038 3.511 3.288

Answer 4

I also struggled with this for some time, trying to calculate the area share of land cover classes of a ~300mx300m grid map in a ~1kmx1km grid. The latter was a polygon file. I tried the multicore solution but this was still too slow for the number of grid cells I had. Instead I:

1. Rasterized the 1kmx1km grid giving all grid cells a unique number
2. Used the allign_rasters (or gdalwarp directly) from the gdalUtils package with the r="near" option to increase the resolution of the 1kmx1km grid to 300mx300m, same projection etc.
3. Stack the 300mx300m land cover map and the 300mx300m grid from step 2, using the raster package: stack_file <- stack(lc, grid).
4. Create a data.frame to combine the maps: df <- as.data.frame(rasterToPoints(stack_file)), which maps the grid cell numbers of the 1kmx1km map to the 300mx300m land cover map
5. Use dplyr to calculate the share of land cover class cells in the 1kmx1km cells.
6. Create a new raster on the basis of step 5 by linking it to the original 1kmx1km grid.

This procedure runs pretty quick and without memory issues on my pc, when I tried it on a land cover map with > 15 mill grid cells at 300mx300m.

I assume the approach above will also work if one wants to combine a polygon file with irregular shapes with raster data. Perhaps, one could combine step 1&2 by directly rasterizing the polygon file to a 300mx300 grid using rasterize (raster probably slow) or gdal_rasterize. In my case I needed to reproject as well so I used gdalwarp to both reproject and disaggregate at the same time.

Answer 5

If a loss in the precision of the overlay is not terribly important - assuming it is precise to begin with - one can typically achieve much greater zonal calculation speeds by first converting the polygons to a raster. The raster package contains the zonal() function, which should work well for the intended task. However, you can always subset two matrices of the same dimension using standard indexing. If you must maintain polygons and you don't mind GIS software, QGIS is actually must faster at zonal statistics than either ArcGIS or ENVI-IDL.

Answer 6

I have to face this same problem to extract values inside polygon from a big mosaic (50k x 50k). My polygon are only have 4 points. The fastest method I found is to crop mosaic into bound of polygon, triangulate polygon into 2 triangles, then check whether points in the triangle (The fastest algorithm I found). Compare with extract function, the run time is reduced from 20 s into 0.5 s. However, the function crop still requires about 2 s for each polygon.

Sorry I cannot provide the full reproducible example. The R codes below don't include the inputs files.

This method is only working for simple polygons.

par_dsm <- function(i, image_tif_name, field_plys2) {
    library(raster)
    image_tif <- raster(image_tif_name)
    coor <- field_plys2@polygons[[i]]@Polygons[[1]]@coords
    ext <- extent(c(min(coor[,1]), max(coor[,1]), min(coor[,2]), max(coor[,2])))

    extract2 <- function(u, v, us, vs) {
        u1 <- us[2]  - us[1]
        u2 <- us[3]  - us[2]
        u3 <- us[1]  - us[3]
        v1 <- vs[1]  - vs[2]
        v2 <- vs[2]  - vs[3]
        v3 <- vs[3]  - vs[1]
        uv1 <- vs[2] * us[1] - vs[1] * us[2]
        uv2 <- vs[3] * us[2] - vs[2] * us[3]
        uv3 <- vs[1] * us[3] - vs[3] * us[1]

        s1 <- v * u1 + u * v1 + uv1
        s2 <- v * u2 + u * v2 + uv2
        s3 <- v * u3 + u * v3 + uv3
        pos <- s1 * s2 > 0 & s2 * s3 > 0
        pos 
    }

    system.time({
        plot_rect <- crop(image_tif, ext, snap ='out')
        system.time({
        cell_idx <- cellFromXY(plot_rect, coor[seq(1,4),])
        row_idx <- rowFromCell(plot_rect, cell_idx)
        col_idx <- colFromCell(plot_rect, cell_idx)

        rect_idx <- expand.grid(lapply(rev(dim(plot_rect)[1:2]), function(x) seq(length.out = x)))

        pixel_idx1 <- extract2(
            rect_idx[,2], rect_idx[,1], 
            row_idx[c(1,2,3)], col_idx[c(1,2,3)])
        pixel_idx2 <- extract2(
            rect_idx[,2], rect_idx[,1], 
            row_idx[c(1,4,3)], col_idx[c(1,4,3)])
        pixel_idx <- pixel_idx1 | pixel_idx2
        })
    })
    mean(values(plot_rect)[pixel_idx])
}

# field_plys2: An object of polygons
# image_taf_name: file name of mosaic file
library(snowfall)
sfInit(cpus = 14, parallel = TRUE)
system.time(plot_dsm <- sfLapply(
    seq(along = field_plys2), par_dsm, image_tif_name, field_plys2))
sfStop()