7

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.

  • 4
    If extract is the bottleneck of your process, have a look at the velox package, which does efficient data extraction – Loïc Dutrieux Mar 8 '17 at 21:09
  • I believe this will certainly do the trick! Thanks. – Christian Sprague Mar 8 '17 at 22:46
4

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.

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.