29

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 
  • You might try this R code profiler (marcodvisser.github.io/aprof/Tutorial.html). It can tell you which lines take most of the time. The link also has guidelines for cutting down processing time in R. – J Kelly Jan 16 '15 at 14:54
  • Just my two cents here . . . but crop/getvalues method does not work when the number of pixels in the crop is very low. I'm not sure where the limit is, but I had issues on crops where there were just 1-5 pixels (I haven't determined the exact reason why (bit new still to spatial packages) but I bet the crop function depends on the pixel boundaries, so thus struggles to crop any individual pixels). Extract from the raster package has no such issue, but agreed is over twice the user time and a lot more than twice on system time. Just a warning to those who have low resolution rasters (and an in – Neal Barsch Apr 9 '17 at 14:13
  • 2
    There is a somewhat new package, velox, that has moved extract into C via the Rcpp package. It is giving ~10 fold increases in speed on extract operations using polygons. – Jeffrey Evans Apr 9 '17 at 14:53
  • @JeffreyEvans . Testing the answer to this question using velox now. Having issues with it allocating extremely large vectors however. – SeldomSeenSlim Oct 3 '17 at 15:41
23

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)
  • 3
    The suggestion that getValues was much faster than the extract does not seem valid because if you use extract you do not have to do crop and rasterize (or mask). The code in the original question does both, and that should about double processing time. – Robert Hijmans Jun 20 '16 at 16:48
  • the only way to know is by testing. – djas Jun 20 '16 at 18:09
  • What class is the polygonlist here, and what should the polygonlist[[i]][,j] do here (ELI5, please)? I'm newbie to parallel stuff, so I don't understand that very well. I couldn't get the function to return anything, until I changed to polygonlist[[i]][,j] to polygonlist[,j], which seems logical becaus [,j] is the jth element of a SpatialPolygonsDataframe, if that is the correct class? After changing that I got the process running and some outputs, but there is definitely still something wrong. ( I try to extract the median value inside n small polygons, so I changed a bit of the code too). – reima Nov 9 '17 at 14:04
  • @RobertH In my case, the cropping (and masking) makes it run about 3 times faster. I'm using a 100 million acre raster and the polygons are a tiny fraction of that. If I don't crop to the polygon, the process runs much slower. Here are my results: clip1 <- crop(rasterlayer, extent(single)) > system.time(ext<-extract(clip1,single)) #extracting from cropped raster user system elapsed 65.94 0.37 67.22 > system.time(ext<-extract(rasterlayer,single)) #extracting from a 100 million acre raster user system elapsed 175.00 4.92 181.10 – Luke Macaulay Jan 26 '18 at 21:10
4

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

  • I ran the two approaches and your method came out slightly slower in my use case. – Luke Macaulay Dec 14 '18 at 21:21
2

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.

2

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.

0

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()

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.