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

[Eelke Folmer]

Hi Holly,
Could you provide some insight into your int object. Does it happen for 
all i?
My first guess is that your matrix 'int' might not be what you expect; 
try: grid[int[,i],] and see whether it is a geometry.
Eelke


Op 8:59 PM, Whiteley,Holly schreef:
> Dear list,
> I was wondering if someone could help me with polygon-polygon spatial
> queries?
>
> Similar to the below request in a previous thread, I have a
> SpatialPolygonDataFrame for an environmental factor (seabed sediment),
> and a gridded SpatialPolygonDataFrame. I am attempting to calculate the
> area of the different seabed sediment types within each cell of the
> polygon grid. I then hope to identify the dominant sediment type in
> terms of area for each cell and write it to a new dataframe against the
> IDs of the cells in the polygon grid.
>
> I started adapting the script that Roger Bivand provided in the below
> message (I hope I have understood correctly), but when I run this line:
>
>  > for (i in seq(along=vec)) vec[[i]] <- try(gIntersection(sediment[i,],
> grid[int[,i],], byid=TRUE))
>
> I get this error message, which I don't understand:
>
> Error in if (is.numeric(i) && i < 0) { :
> missing value where TRUE/FALSE needed
>
> I'd really appreciate any help you could offer me for getting round this
> issue.
>
> Many thanks for your time,
> Holly
>
> -----------------------------------
> Original message
> -----------------------------------
>
> Re: [R-sig-Geo] Calculate the surface of polygons within gridcells with
> gIntersection
> Roger Bivand
> Wed, 29 Jun 2011 06:43:29 -0700
>
> 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
>
>