[R-sig-Geo] Calculate the surface of polygons within gridcells with gIntersection

[Eelke Folmer]

Hi all,

I have a polygon shapefile with a column in the attribute table denoting 
class (in my case "habitat"). I construct a gridded 
SpatialPolygonsDataFrame. For each gridcell I'd like to be able to 
obtain the surface of each of the classes from the polygon shapefile 
resulting in something like this:

cell	classA	classB	classC
1	10	25	31
2	53	23	0
3	13	25	..
4	..	..	..	

I tried by means of gIntersection (rgeos) to make a new geometry but I 
run into difficulties (see below for a reproducable example). If I were 
to obtain a new geometry I would probably be able to match the grid with 
the area (perhaps looping through the classes?).
I'd suspect that this is a relatively common operation and I wonder if 
there is a more direct way (within R).

Example:
require(raster)
require(rgeos)

ext = extent(c(0, 10, 0, 10))
grd     <- raster(ext, nrows=10, ncols=10, crs=NA )
grd.pol <- rasterToPolygons(grd, fun=NULL, n=4, na.rm=TRUE, digits=12)

grd2     <- raster(ext-1.5, nrows=5, ncols=5, crs=NA )
grd.pol2 <- rasterToPolygons(grd2, fun=NULL, n=4, na.rm=TRUE, digits=12)

plot(grd.pol)
plot(grd.pol2, add=T, lty=2)

int <- gIntersection(grd.pol2, grd.pol, byid=T)

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, 
"rgeos_intersection") : Geometry collections may not contain other 
geometry collections


I have tried with various other geometries (polygon shapefiles) and run 
into the same problem. I have also searched the the mailing list for a 
hint how to carry on but with no luck.

Thanks for a powerful library anyway!

Eelke