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
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
