Help needed with extraction of raster statistics by polygons
Denys, The GIS functions in R have great utility, but speed is an issue; extract() in particular. So, if you make a 30m raster equivalent to the Grid30m but with the cell numbers as the values, you can use extract just once to determine which cells are to be extracted with each polygon. Then use regular indexing to find the mean for the raster cells covered by each polygon. This should be much faster. Dominik Dominik Schneider o 303.735.6296 | c 518.956.3978 On Tue, Aug 25, 2015 at 9:26 AM, Denys Dukhovnov via R-sig-Geo <
r-sig-geo at r-project.org> wrote:
Hello all! I need to spatially calculate 30-meter grid raster mean over US Census blocks in North East region (approx. 1.9 mln street blocks/polygons). I came up with the script that runs in parallel, but it still takes days to complete, even when the i7 CPU is used up to the max. My question: is there any other way to improve the processing speed given the set-up below (I included the backbone of the script to save space). I am new to R, so any help will be much appreciated. Thanks a lot! Regards, Denys
____________________________
library(sp)
library(rgdal)
library(raster)
library(foreach)
library(spatial.tools)
Grid30m <- raster("raster.tif") # loading main raster into R
NEfips <- c(09, 10, 11, ....) # list of state FIPS codes, whose
street block polygons are to be processed
ShapePath <- "T:\\..." # path to block shapefiles
sfQuickInit(cpus=8) # setting up and registering parallel
processing back-end
for (x in NEfips) {
State <- paste("block_", x , sep="")
Blocks <- readOGR(dsn=ShapePath, layer=State[1],
stringsAsFactors=F) # loading the block polygons
BlocksReproject <- spTransform(Blocks, crs(Grid30m))
GridMean <- foreach (i=1:length(BlocksReproject$GEOID10),
.combine='c', .packages="raster") %dopar% { # Parallel processing
loop for extracting mean raster value per street block polygon
extract(Grid30m, BlocksReproject[i], fun=mean)
}
write.table(...) # Generating output containing the mean raster
statistics for each block polygon.
}
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo