[GIS] How to speed up extract(rstr,shp) ? – R

extractpolygonrrasterrasterization

I am trying to create a function that takes in a shapefile and a binary (1's and 0's) raster, with the same extent, and extracts the values from the raster. Once this occurs I desire to sum up those values and get a percentage of pixels within those boundaries that are 1. Then rasterize these percentages. Here is what I have so far. It is working just fine BUT IS EXTREMELY SLOW! 

Density2 <- function(rstr,shp,... ){

UrbanBuildup <- extract(rstr,shp)

summed.values <- lapply(UrbanBuildup, FUN = sum) 
#sum the values in    each respective bndry
number.values <- lapply(UrbanBuildup, FUN = length)

vector.of.summed.values <- matrix(summed.values,nrow=length(summed.values),ncol=1)
vector.of.number.values <- matrix(number.values,nrow=length(number.values),ncol=1)

vector.of.summed.values <- unlist(vector.of.summed.values)
vector.of.number.values <- unlist(vector.of.number.values)

sumUrbanBuildup <- vector.of.summed.values
sumUrbanBuildup[which(is.na(sumUrbanBuildup))] <- 0
numPixels <- vector.of.number.values

shp$den2 <- sumUrbanBuildup/numPixels
shp$den2[which(is.infinite(shp$den2))] <- 0
shp$den2[which(is.na(shp$den2))] <- 0

ext <- extent(shp)
r <- raster(ext, res=300) 
r <- rasterize(shp, r, field='den2')
return(r)
}

I read elsewhere it is faster to rasterize() the polygon first and use getValues() but I do not know how to apply this to my situation.

The project is for urban build-up/density using remote sensing data, hence the names.

Best Answer

Try mask(raster, shape). I did something similar a while back and it worked very fast for my purpose: Extract Raster from Raster using Polygon shapefile in R

From the package description:

Create a new Raster* object that has the same values as x, except for the cells that are NA (or other maskvalue) in a ’mask’. These cells become NA (or other updatevalue). The mask can be either another Raster* object of the same extent and resolution, or a Spatial* object (e.g. SpatialPolygons) in which case all cells that are not covered by the Spatial object are set to updatevalue. [...]

Simple code example from the documentation:

r <- raster(ncol=10, nrow=10)
m <- raster(ncol=10, nrow=10)
r[] <- runif(ncell(r)) * 10
m[] <- runif(ncell(r))
m[m < 0.5] <- NA
mr <- mask(r, m)
m2 <- m > .7
mr2 <- mask(r, m2, maskvalue=TRUE)

For more information just search for the mask function in the documentation of the raster package.

Related Question