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.
extract
is the bottleneck of your process, have a look at the velox package, which does efficient data extraction