[R-sig-Geo] unpacking polygon geometry with holes

[Roger Bivand]

On Sun, 22 Apr 2012, Colin Rundel wrote:

> Hi Colin,
>
> You can try the following code to separate the resulting polygons and 
> their holes from the other polygons. Basically you can find the parent / 
> hole relationship by checking the comment attribute of the Polygon 
> object. Currently it is a kludgy hack since sp objects were not 
> originally aware of this relationship, but hopefully it is something 
> that will be improved on in future versions of sp and related packages 
> (see Roger's recent work on readOGR for reading in polygons).

Hi Colins,

I'm afraid the recent work is unlikely to help here. The underlying 
problem is that the data are returning nested geometry collections, which 
isn't predictable ahead of time. I tend then to look for a logical matrix 
to iterate over by row to disambiguate, and then stick the output together 
again afterwards - the row.names could be improved:

load(url("http://dl.dropbox.com/u/8776628/polys.rdata"))
library(rgeos)
oo <- gOverlaps(p1, p2, byid=TRUE)
res <- vector(mode="list", length=nrow(oo))
for (i in seq(along=res)) {
   if(any(oo[i,])) {
     res[[i]] <- gDifference(p1[oo[i,]], p2)
     row.names(res[[i]]) <- paste(i, row.names(res[[i]]), sep="_")
   }
}
res1 <- res[!sapply(res, is.null)]
out <- do.call("rbind", res1)
plot(p1, col="grey")
plot(p2, border="red", add=TRUE, angle=45, density=10, col="red")
plot(out, border="green", add=TRUE, angle=135, density=10, col="green")

Your reconstruction:

plot(pList, add=TRUE, border="blue", angle=215, density=15, col="blue")

isn't quite the same, because p1 geometries that don't overlap with p2 
geometries are included, and I'm not sure what the actual one-sided 
difference should be. I'm also uncertain whether this was intended to be a 
one-sided difference, or a symmetric difference. I'm wondering whether 
there is a way to disambiguate the nested geometry collections internally.

res <- vector(mode="list", length=length(slot(p1, "polygons")))
for (i in seq(along=res)) {
   res[[i]] <- gDifference(p1[i], p2)
   if (!is.null(res[[i]]))
     row.names(res[[i]]) <- paste(i, row.names(res[[i]]), sep="_")
}
res1 <- res[!sapply(res, is.null)]
out <- do.call("rbind", res1)

is the same as pList in total area, so may do the same. I can see two 
cases in which p1/p2 intersections are in the output object.

Roger


>
>
> p_diff = gDifference(p1, p2, byid=F)
>
> x = (slot(p_diff, 'polygons')[[1]])
>
> comm = as.numeric(strsplit(comment(x)," ")[[1]])
>
> parents = which(comm == 0)
> holes = which(comm != 0)
>
> for(i in 1:length(parents)) {
>    p = parents[i]
>    polys = c(p, which(comm == p))
>    xx[[i]] = Polygons(x at Polygons[polys], i)
> }
>
> pList = as.SpatialPolygons.PolygonsList(xx)
>
> -Colin
>
>
> On Apr 21, 2012, at 11:39 AM, Colin Robertson wrote:
>
>> Dear List,
>>
>> I am trying to manipulate polygons that are the component pieces that
>> result from the gDifference function.
>>
>>> spdf1 = readShapePoly(paste(pathToDropbox, "spatiallab/hschange/hs99.shp", sep="")) #rgdal wouldnt insstall
>>> spdf2 = readShapePoly(paste(pathToDropbox, "spatiallab/hschange/hs00.shp", sep=""))
>>> p1 = polygons(spdf1)
>>> p2 = polygons(spdf2)
>>
>> The above data is in a workspace here:
>> http://dl.dropbox.com/u/8776628/polys.rdata
>>
>> If I use byid=T in gDifference I get an error:
>>
>> p_diff = gDifference(p1, p2, byid=T)
>>
>> Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_difference") :
>>  Geometry collections may not contain other geometry collections
>>
>> Even though both are SpatialPolygons:
>> class(p1)
>> [1] "SpatialPolygons"
>> attr(,"package")
>> [1] "sp"
>> class(p2)
>> [1] "SpatialPolygons"
>> attr(,"package")
>> [1] "sp"
>>
>> So I am trying to work with the result of byid=F and unpacking it using:
>>
>> #unpack the polygons to make them useful (find better way to do this)
>> x <- (slot(p_diff, 'polygons')[[1]])
>> xx <- slot(x, "Polygons")
>> for(i in 1:length(xx)) {
>> 	xx[[i]] <- Polygons(x at Polygons[i], i)
>> }
>>
>> pList <- as.SpatialPolygons.PolygonsList(xx)
>> p_diff2 <- pList
>> length(p_diff2)
>> [1] 39
>>
>> Which works except for interior rings are not preserved (there are two)
>>
>> which(unlist(lapply(slot(slot(p_diff, 'polygons')[[1]], 'Polygons'),
>> function(P) P at hole)) == TRUE)
>> [1] 14 29
>>
>> Attempting to use checkPolygonsHoles has not worked because each of
>> the polygons are an individual Polygons object at this point, but I
>> cant figure out how to keep interior rings associated with their
>> parent polygons.
>>
>> polysList <- as.SpatialPolygons.PolygonsList(xx)
>> chkList <- lapply(polysList at polygons, checkPolygonsHoles)
>> polysList <- as.SpatialPolygons.PolygonsList(chkList)
>>
>> plot(polysList,col="red", pbg="white") #interior polys are red not white
>>
>>
>> Any assistance much appreciated -
>>
>> Best,
>>
>> Colin
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> _______________________________________________
> 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