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

[Colin Robertson]

Roger and Colin,

Thanks very much for this.

Roger you second method seems to get me just about there, although I
did notice one instance of one 'multi-part' polygon in the output.
Where two parts of one input polygon were the 'difference'. I am think
a little more about whether single-parts are really needed.

In response to your question, I actually am after the symmetric
difference, but I found that using gSymmdifference would dissolve
boundaries between coincident polygons in the output, and I need them
as separate geometries. So instead I am running gDifference twice and
combining difference1, difference2, and the intersection to give me
all component geometries of a union of p1 and p2 with boundaries in
tact, and each geometry 'single part' with its own id.

Thanks again -

Colin


On Sun, Apr 22, 2012 at 10:36 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>