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

[Whiteley,Holly]

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