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

[Colin Robertson]

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