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

[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:

> 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