Recombining polygon shapefiles using maptools

On Thu, 27 Mar 2008, Don MacQueen wrote:

I have six polygon shapefiles. Two of them represent my area of
interest (call them A and B), and the other four (call them C, D, E,
F) represent holes in the first two.

I would like to create a single object that can be passed to
spsample() for spatial sampling, such that spsample will place
samples inside A and B, but not in C, D, E, or F.

I know how to do this by what might be called "brute force" (see below).

Sorry for not replying earlier. For this case, brute force may be the most
suitable.

The real question is, are there ways to do this more effectively
using higher level functions?

The higher level methods (spCbind() in maptools) are for Polygons objects
rather than for Polygon objects, because spCbind() expects to cbind two
lists of Polygons objects and two data frames.

The difficulty here is to find out how to pack and unpack your geometries
to use checkPolygonsHoles() in maptools. If you can put all your
geometries into a single Polygons object, checkPolygonsHoles() will return
a single Polygons object with the holes correctly identified, and that
will work with spsample. It will, however, treat the sammpled points as
lying within the same Polygons object, but maybe that doesn't matter.

This is untried:

Do spCbind() on the 5 SpatialPolygonDataFrames after having given the
constituent Polygons objects unique IDs (spChFIDs() method).

ALL <- spCbind(spCbind(spCbind(spCbind(A, B), C), D), E)

Add a constant vector to the output object, and use it as the IDs=
argument to unionSpatialPolygons()

ALL$all <- 1
out <- unionSpatialPolygons(as(ALL, "SpatialPolygons"), IDs=ALL$all)

Check out - it may be that the first pass through gpclib will be enough,
or

out1 <- sapply(slot(out, "polygons"), checkPolygonsHoles)

where out1 will be a list of Polygons object of length

length(slot(out, "polygons"))

If only one, just use spsample() on that (there is a sample.Polygons()
method), if more than one, build a SpatialPolygons object, and use
spsample() on that.

Have you considered using the spsurvey package - it is more targetted than
spsample() methods - or does spsample() meet your needs?

Hope this helps,

Roger

If there were, it might make for easier to understand scripts, for
example, or be easier to repeat using different sets of shapefiles
(the script below doesn't easily generalize, especially if any of the
shapefiles consist of multiple polygons).

Thanks
-Don



Here is my solution; I've run it and it works. I apologize for not
being able to supply the shapefiles and thus a reproducible example.

Each shapefile consists of a single polygon, and I don't need any of
the attribute information from the shapefiles.
This simplifies things, quite a lot, I think.

Extract the single polygon from each, into six separate two column matrices.

# A
vz1 <- readShapePoly('shapefiles-zones/Zone-TK')
tmp <- as(vz1 , 'SpatialPolygons')
tmp <- tmp at polygons[[1]]   ## since I know it has only one polygon
poly1 <- tmp at Polygons[[1]]@coords  ## a matrix of coordinates

## repeat for additional shapefiles 2 through 6


Combine the polygons following the example in ?overlay

## this example uses only the first three of my polygons
tmp <- Polygons(
                list(Polygon(poly1,hole=FALSE),  # A
                      Polygon(poly2,hole=FALSE),  # B
                      Polygon(poly3,hole=TRUE)),  # C
                ID=1)
sr <- SpatialPolygons(list(tmp))


plot(sr)
tmp <- spsample(sr, type='random', n=500)
points(tmp)  ## looks good!