Skip to content

unpacking polygon geometry with holes

4 messages · Colin Robertson, Colin Rundel, Roger Bivand

Original
Colin Robertson
# 
Dear List,

I am trying to manipulate polygons that are the component pieces that
result from the gDifference function.
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
Colin Rundel
# 
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).


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:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Sun, 22 Apr 2012, Colin Rundel wrote:

        

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      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: