portion of a raster cell covered by one or more polygons: is there a faster way to do this (in R)? -

pictures better words, please have @

what have

• a rasterlayer object (filled random values here illustration purposes only, actual values don't matter)
• a spatialpolygons object lots , lots of polygons in

you can re-create example data used image following code:

library(sp) library(raster) library(rgeos)  # create example raster r <- raster(nrows=10, ncol=15, xmn=0, ymn=0) values(r) <- sample(x=1:1000, size=150)  # create example (spatial) polygons p1 <- polygon(coords=matrix(c(50, 100, 100, 50, 50, 15, 15, 35, 35, 15), nrow=5, ncol=2), hole=false) p2 <- polygon(coords=matrix(c(77, 123, 111, 77, 43, 57, 66, 43), nrow=4, ncol=2), hole=false) p3 <- polygon(coords=matrix(c(110, 125, 125, 110, 67, 75, 80, 67), nrow=4, ncol=2), hole=false)  lots.of.polygons <- spatialpolygons(list(polygons(list(p1, p2, p3), 1))) crs(lots.of.polygons) <- crs(r) # copy crs raster polygons (please ignore potential problems related projections etc. now)   # plot both plot(r) #values in raster illustration purposes plot(lots.of.polygons, add=true) 

for each cell in raster, want know how of covered 1 or more polygons. or actually: area of polygons within raster cell, without outside cell in question. if there multiple polygons overlapping cell, need combined area.

the following code want, takes more week run actual data sets:

# empty example raster (don't need values): values(r) <- na  # copy of r hold results r.results <- r  (i in 1:ncell(r)){   r.cell <- r # fresh copy of empty raster   r.cell[i] <- 1 # set ith cell 1   p <- rastertopolygons(r.cell) # create polygon represents i-th raster cell   cropped.polygons <- gintersection(p, lots.of.polygons) # intersection of i-th raster cell , spatialpolygons    if (is.null(cropped.polygons)) {     r.results[i] <- na # if there's no polygon intersecting raster cell, return na ...   } else{     r.results[i] <- garea(cropped.polygons) # ... otherwise return area   } }  plot(r.results) plot(lots.of.polygons, add=true) 

i can squeeze out bit more speed using sapply instead of for-loop, bottleneck seems somewhere else. whole approach feels pretty awkward, , i'm wondering if i've missed obvious. @ first thought rasterize() should able easily, can't figure out put fun= argument. ideas?

[edited]

maybe gintersection(..., byid = t) gunaryunion(lots.of.polygons) (they enable treat cells @ once) faster loop (if gunaryunion() takes time, bad idea).

r <- raster(nrows=10, ncol=15, xmn=0, ymn=0) set.seed(1); values(r) <- sample(x=1:1000, size=150) rr <- rastertopolygons(r)  # joining intersecting polys , put polys single spatialpolygons lots.of.polygons <- gunaryunion(lots.of.polygons)   # in example, unnecessary  gi <- gintersection(rr, lots.of.polygons, byid = t)  ind <- as.numeric(do.call(rbind, strsplit(names(gi), " "))[,1])   # getting intersected rr's id r[] <- na r[ind] <- sapply(gi@polygons, function(x) slot(x, 'area'))  # bit faster garea(gi, byid = t)  plot(r) plot(lots.of.polygons, add=true)