[R-sig-Geo] How to efficiently clip large grids/raster with polygons?

[Ariel Ortiz-Bobea]

Hello everyone,

I've been exploring ways to "clip a grid/raster with a polygon/shape file"
in R but what I have seen so far seems more suited for relatively fine grids
and large polygons (e.g. many grids fall within each polygon). 

I'm confronted with a relatively coarse data grid for the US (0.3 deg ~ 32
km) and would like to clip these grids using county boundaries and calculate
mean values for each polygon/county. Basically, many counties will fall "in
between" several grids and so, ideally, I need to calculate a weighted mean
for each polygon by weighting each grid value by the fraction of the polygon
it is covering.

I found a way to do this with the "extract" function in the raster package
using the "weight" option (see example below and ?extract), but it takes
about 60 minutes for one grid/raster layer for all US counties (~ 1
second/county/layer). My problem is that I need to do this for about 30
years * 365 days * 8 daily snapshots = 87,600 layers / variable. Enough to
keep me 8+ years waiting for each variable, and I have about 25 variables to
process...

Do you know a (hopefully way) more efficient way to do this? Thank you in
advance for any guidance!

Ariel

PS: I'm running on a 2.66 GHz core 2 duo iMac / 8 GB RAM / Mac OS X 10.7
(Lion) / R 2.12.2

FILES:
-------
shape file: http://terpconnect.umd.edu/~ortiz/R/US48_rr-fixed.grb (~500KB)
GRIB grid file: http://terpconnect.umd.edu/~ortiz/R/US48_rr-fixed.grb
(~12MB)

CODE:
-------
#used libraries
library(rgdal) 
library(raster)

#setting up data
grid <- readGDAL("US48_rr-fixed.grb")
proj <- CRS(proj4string(fixed))
countyshape0	<- readShapePoly("co99_d00.shp", proj4string =
CRS("+proj=longlat"))
countyshape 	<- spTransform(countyshape0, proj) 
removestate	<- countyshape$STATE %in% c("02","15","72") # leave out Alaska,
Hawai and Puerto Rico
countyshape48	<- countyshape[!removestate,]
rasteredgrid <-raster(grid, 3) # chooses one layer

#visual glimpse:
image(rasteredgrid)
plot(countyshape48, add=TRUE)

#cookie cutting/raster clipping
timer <- proc.time() 	# start the clock
test <- extract(flxraster, countyshape48, weights=TRUE, fun=mean) 
proc.time() - timer 	# stop the clock

-----
Ariel Ortiz-Bobea
PhD Candidate in Agricultural & Resource Economics
University of Maryland - College Park
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/How-to-efficiently-clip-large-grids-raster-with-polygons-tp6804923p6804923.html
Sent from the R-sig-geo mailing list archive at Nabble.com.