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

[Roger Bivand]

On Wed, 29 Jun 2011, Eelke Folmer wrote:

> 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


With byid=TRUE, you may end up with rgeos trying to re-assemble things in 
a different way than you intended. With reference to an earlier thread 
this month, if you do:

int <- gIntersects(grd.pol2, grd.pol, byid=TRUE)
vec <- vector(mode="list", length=dim(int)[2])
for (i in seq(along=vec)) vec[[i]] <- try(gIntersection(grd.pol2[i,], 
grd.pol[int[,i],], byid=TRUE))
out <- do.call("rbind", vec)
rn <- row.names(out)
nrn <- do.call("rbind", strsplit(rn, " "))
df <- data.frame(pol2=nrn[,1], pol=nrn[,2], area=sapply(slot(out,
   "polygons"), slot, "area"))
plot(grd.pol, border="red")
plot(grd.pol2, border="green", add=TRUE)
plot(out, col="#00FFFF64", add=TRUE, border="transparent")
text(coordinates(grd.pol), labels=row.names(grd.pol), cex=0.5)
head(df)

you can see that you are now getting a Polygons object for each 
intersection, which you should be able to tally.

Using a predicate first (and it is also possible to use STRtrees) also 
saves time, but we wanted to get the standard facilities of rgeos working 
before optimising, so the optimisation should be done by hand.

Hope this helps,

Roger

>
>
> 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
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no