[R-sig-Geo] Recombining polygon shapefiles using maptools
Roger Bivand
Roger.Bivand at nhh.no
Wed Apr 2 11:53:10 CEST 2008
On Wed, 2 Apr 2008, Roger Bivand wrote:
> 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)
Certainly untried - not spCbind() methods, but, of course, spRbind()
methods to bind the *rows*, sorry.
Roger
>
> 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
More information about the R-sig-Geo
mailing list