[R-sig-Geo] Recombining polygon shapefiles using maptools

[Roger Bivand]

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!
>
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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