[R-sig-Geo] How to extract coordinates values from a shapefile?

[Paul Hiemstra]

On 06/09/2010 11:01 PM, Dylan Beaudette wrote:
> On Wednesday 09 June 2010, Hadley Wickham wrote:
>    
>> On Wed, Jun 9, 2010 at 3:24 PM, Edzer Pebesma
>>
>> <edzer.pebesma at uni-muenster.de>  wrote:
>>      
>>> The example provided by Matt assumes that each polygon consists of a
>>> single ring, and doesn't have islands, lakes etc. The function below
>>> pasts all coordinates to a single 2-column matrix. For clarity's sake, I
>>> avoided nested sapply's.
>>>
>>> library(maptools)
>>> xx<- readShapePoly(system.file("shapes/sids.shp",
>>> package="maptools")[1], IDvar="FIPSNO", proj4string=CRS("+proj=longlat
>>> +ellps=clrk66")) allcoordinates = function(x) {
>>>     ret = NULL
>>>     polys = x at polygons
>>>     for(i in 1:length(polys)) {
>>>         pp = polys[[i]]@Polygons
>>>         for (j in 1:length(pp))
>>>             ret = rbind(ret, coordinates(pp[[j]]))
>>>     }
>>>     ret
>>> }
>>> q = allcoordinates(xx)
>>>
>>> # check:
>>> plot(xx)
>>> lines(q, col='blue')
>>>        
>> And that's basically what fortify does, except it covers a few more cases.
>>
>> Hadley
>>      
> Not to be nit-picky, but wouldn't it be safer to pre-allocated "ret" instead
> of dynamically appending? This would allow the suggested function to scale to
> very large geometries.
>    
A combination of lapply and do.call("rbind" also works well in terms of 
performance.

cheers,
Paul
>
> Cheers,
> Dylan
>
>
>
>    


-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 274 3113 Mon-Tue
Phone:  +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul
http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770