Quantcast

How to extract coordinates values from a shapefile?

classic Classic list List threaded Threaded
16 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

How to extract coordinates values from a shapefile?

Thiago V. dos Santos
  Dear R colleagues,
  Does anyone know if it's possible to create a vector with coordinate values extracted from a shape loaded with readShapePoly command of "maptools" package?
  For example, the following code loads the shapefile containing Brazilian states division:
library(maptools)br<-readShapePoly("BR-map/3.shp")
  What I need to do is create a vector which contains all lat and lon values inside this shapefile...
  Thanks in advance,
  Thiago Veloso.


     
        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Henrique Dallazuanna
Try:

coordinates(br)

On Wed, Jun 9, 2010 at 3:45 PM, Thiago Veloso <[hidden email]>wrote:

>   Dear R colleagues,
>   Does anyone know if it's possible to create a vector with coordinate
> values extracted from a shape loaded with readShapePoly command of
> "maptools" package?
>   For example, the following code loads the shapefile containing Brazilian
> states division:
> library(maptools)br<-readShapePoly("BR-map/3.shp")
>   What I need to do is create a vector which contains all lat and lon
> values inside this shapefile...
>   Thanks in advance,
>   Thiago Veloso.
>
>
>
>        [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>

--
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O

        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Nikhil Kaza-2
In reply to this post by Thiago V. dos Santos
Look at this. Not sure if this is indeed what you want or if you  
actually have to unproject them.

http://www.mail-archive.com/r-sig-geo@.../msg04715.html

require(ggplot2)
fortify(br)


Nikhil Kaza
Asst. Professor,
City and Regional Planning
University of North Carolina

[hidden email]

On Jun 9, 2010, at 2:45 PM, Thiago Veloso wrote:

>   Dear R colleagues,
>   Does anyone know if it's possible to create a vector with  
> coordinate values extracted from a shape loaded with readShapePoly  
> command of "maptools" package?
>   For example, the following code loads the shapefile containing  
> Brazilian states division:
> library(maptools)br<-readShapePoly("BR-map/3.shp")
>   What I need to do is create a vector which contains all lat and  
> lon values inside this shapefile...
>   Thanks in advance,
>   Thiago Veloso.
>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Thiago V. dos Santos
In reply to this post by Henrique Dallazuanna
  Thank you, Henrique.
  This command returned a short list of coordinates (maybe the limits?):
>coordinates(br)          [,1]      [,2][1,] -53.31971 -29.70562[2,] -50.48913 -27.24445> 
  Any other suggestion?
  Cheers,
  Thiago.

--- On Wed, 9/6/10, Henrique Dallazuanna <[hidden email]> wrote:

From: Henrique Dallazuanna <[hidden email]>
Subject: Re: [R-sig-Geo] How to extract coordinates values from a shapefile?
To: "Thiago Veloso" <[hidden email]>
Cc: [hidden email]
Date: Wednesday, 9 June, 2010, 15:53

Try:

coordinates(br)

On Wed, Jun 9, 2010 at 3:45 PM, Thiago Veloso <[hidden email]> wrote:


  Dear R colleagues,

  Does anyone know if it's possible to create a vector with coordinate values extracted from a shape loaded with readShapePoly command of "maptools" package?

  For example, the following code loads the shapefile containing Brazilian states division:

library(maptools)br<-readShapePoly("BR-map/3.shp")

  What I need to do is create a vector which contains all lat and lon values inside this shapefile...

  Thanks in advance,

  Thiago Veloso.







        [[alternative HTML version deleted]]




_______________________________________________

R-sig-Geo mailing list

[hidden email]

https://stat.ethz.ch/mailman/listinfo/r-sig-geo





--
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O




     
        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Thiago V. dos Santos
In reply to this post by Nikhil Kaza-2
  Dear Nikhil,
  Thank you for your input. I tried your suggestion and received the following error:
> fortify(br)Using UF_05 to define regions.Error: isTRUE(gpclibPermitStatus()) is not TRUE>
  Do you know what else could be done?
  Best wishes,
  Thiago.
--- On Wed, 9/6/10, Nikhil Kaza <[hidden email]> wrote:

From: Nikhil Kaza <[hidden email]>
Subject: Re: [R-sig-Geo] How to extract coordinates values from a shapefile?
To: "Thiago Veloso" <[hidden email]>
Cc: [hidden email]
Date: Wednesday, 9 June, 2010, 16:00

Look at this. Not sure if this is indeed what you want or if you actually have to unproject them.

http://www.mail-archive.com/r-sig-geo@.../msg04715.html

require(ggplot2)
fortify(br)


Nikhil Kaza
Asst. Professor,
City and Regional Planning
University of North Carolina

[hidden email]

On Jun 9, 2010, at 2:45 PM, Thiago Veloso wrote:

>   Dear R colleagues,
>   Does anyone know if it's possible to create a vector with coordinate values extracted from a shape loaded with readShapePoly command of "maptools" package?
>   For example, the following code loads the shapefile containing Brazilian states division:
> library(maptools)br<-readShapePoly("BR-map/3.shp")
>   What I need to do is create a vector which contains all lat and lon values inside this shapefile...
>   Thanks in advance,
>   Thiago Veloso.
>
>
>
>     [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



     
        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Nikhil Kaza-2
You need to execute gpclibPermit()  to enable gpclib.

library(maptools) should have issued a warning to that effect.


Nikhil Kaza
Asst. Professor,
City and Regional Planning
University of North Carolina

[hidden email]

On Jun 9, 2010, at 3:34 PM, Thiago Veloso wrote:

>   Dear Nikhil,
>
>   Thank you for your input. I tried your suggestion and received the  
> following error:
>
> > fortify(br)
> Using UF_05 to define regions.
> Error: isTRUE(gpclibPermitStatus()) is not TRUE
> >
>
>   Do you know what else could be done?
>
>   Best wishes,
>
>   Thiago.
>
> --- On Wed, 9/6/10, Nikhil Kaza <[hidden email]> wrote:
>
> From: Nikhil Kaza <[hidden email]>
> Subject: Re: [R-sig-Geo] How to extract coordinates values from a  
> shapefile?
> To: "Thiago Veloso" <[hidden email]>
> Cc: [hidden email]
> Date: Wednesday, 9 June, 2010, 16:00
>
> Look at this. Not sure if this is indeed what you want or if you  
> actually have to unproject them.
>
> http://www.mail-archive.com/r-sig-geo@.../msg04715.html
>
> require(ggplot2)
> fortify(br)
>
>
> Nikhil Kaza
> Asst. Professor,
> City and Regional Planning
> University of North Carolina
>
> [hidden email]
>
> On Jun 9, 2010, at 2:45 PM, Thiago Veloso wrote:
>
> >   Dear R colleagues,
> >   Does anyone know if it's possible to create a vector with  
> coordinate values extracted from a shape loaded with readShapePoly  
> command of "maptools" package?
> >   For example, the following code loads the shapefile containing  
> Brazilian states division:
> > library(maptools)br<-readShapePoly("BR-map/3.shp")
> >   What I need to do is create a vector which contains all lat and  
> lon values inside this shapefile...
> >   Thanks in advance,
> >   Thiago Veloso.
> >
> >
> >
> >     [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Thiago V. dos Santos
  Prof. Kaza,
  Thank you very much. It worked like a charm!
  Best wishes,
  Thiago.

--- On Wed, 9/6/10, Nikhil Kaza <[hidden email]> wrote:

From: Nikhil Kaza <[hidden email]>
Subject: Re: [R-sig-Geo] How to extract coordinates values from a shapefile?
To: "Thiago Veloso" <[hidden email]>
Cc: [hidden email]
Date: Wednesday, 9 June, 2010, 16:37

You need to execute gpclibPermit()  to enable gpclib. 
library(maptools) should have issued a warning to that effect.

 Nikhil KazaAsst. Professor,City and Regional PlanningUniversity of North Carolina
[hidden email]
On Jun 9, 2010, at 3:34 PM, Thiago Veloso wrote:
  Dear Nikhil,
  Thank you for your input. I tried your suggestion and received the following error:
> fortify(br)Using UF_05 to define regions.Error: isTRUE(gpclibPermitStatus()) is not TRUE>
  Do you know what else could be done?
  Best wishes,
  Thiago.
--- On Wed, 9/6/10, Nikhil Kaza <[hidden email]> wrote:

From: Nikhil Kaza <[hidden email]>
Subject: Re: [R-sig-Geo] How to extract coordinates values from a shapefile?
To: "Thiago Veloso" <[hidden email]>
Cc: [hidden email]
Date: Wednesday, 9 June, 2010, 16:00

Look at this. Not sure if this is indeed what you want or if you actually have to unproject them.

http://www.mail-archive.com/r-sig-geo@.../msg04715.html

require(ggplot2)
fortify(br)


Nikhil Kaza
Asst. Professor,
City and Regional Planning
University of North Carolina

[hidden email]

On Jun 9, 2010, at 2:45 PM, Thiago Veloso wrote:

>   Dear R colleagues,
>   Does anyone know if it's possible to create a vector with coordinate values extracted from a shape loaded with readShapePoly command of "maptools" package?
>   For example, the following code loads the shapefile containing Brazilian states division:
> library(maptools)br<-readShapePoly("BR-map/3.shp")
>   What I need to do is create a vector which contains all lat and lon values inside this shapefile...
>   Thanks in advance,
>   Thiago Veloso.
>
>
>
>     [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

       



     
        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Matt Beard
>
> library(maptools)
>
> # A simple way to print out a list of coordinates for each polygon in your
> shapefile:
>
> # Path and filename of polygon shapefile
> testfile <- '/media/PKBACK# 001/FNR210/QGISLab/habitats/habitats.shp'
>
> # Read in polygon shapefile using handy maptools function
> test <- readShapePoly(testfile)
>
> # Extract the list of Polygons objects
> polys <- slot(test,"polygons")
>
> # Within each Polygons object
> #    Extract the Polygon list (assuming just one per Polygons)
> #    And Extract the coordinates from each Polygon
> for (i in 1:length(polys)) {
>    print(paste("Polygon #",i))
>    print(slot(slot(polys[[i]],"Polygons")[[1]],"coords"  ))
> }
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

edzer
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@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')


On 06/09/2010 09:58 PM, Matt Beard wrote:

>>
>> library(maptools)
>>
>> # A simple way to print out a list of coordinates for each polygon in your
>> shapefile:
>>
>> # Path and filename of polygon shapefile
>> testfile <- '/media/PKBACK# 001/FNR210/QGISLab/habitats/habitats.shp'
>>
>> # Read in polygon shapefile using handy maptools function
>> test <- readShapePoly(testfile)
>>
>> # Extract the list of Polygons objects
>> polys <- slot(test,"polygons")
>>
>> # Within each Polygons object
>> #    Extract the Polygon list (assuming just one per Polygons)
>> #    And Extract the coordinates from each Polygon
>> for (i in 1:length(polys)) {
>>    print(paste("Polygon #",i))
>>    print(slot(slot(polys[[i]],"Polygons")[[1]],"coords"  ))
>> }
>>
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Hadley Wickham-2
On Wed, Jun 9, 2010 at 3:24 PM, Edzer Pebesma
<[hidden email]> 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@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

--
Assistant Professor / Dobelman Family Junior Chair
Department of Statistics / Rice University
http://had.co.nz/

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Dylan Beaudette
On Wednesday 09 June 2010, Hadley Wickham wrote:

> On Wed, Jun 9, 2010 at 3:24 PM, Edzer Pebesma
>
> <[hidden email]> 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@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.


Cheers,
Dylan



--
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: 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
>>
>> <[hidden email]>  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@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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Paul Hiemstra
In reply to this post by Dylan Beaudette
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
>>
>> <[hidden email]>  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@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.
>    
In addition to my previous suggestion, I tried using lapply + do.call
instead of the for loop based implementation of edzer:

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@polygons
     for(i in 1:length(polys)) {
         pp = polys[[i]]@Polygons
         for (j in 1:length(pp))
             ret = rbind(ret, coordinates(pp[[j]]))
     }
     ret
}

allcoordinates_lapply = function(x) {
     polys = x@polygons
     return(do.call("rbind", lapply(polys, function(pp) {
     do.call("rbind", lapply(pp@Polygons, coordinates))
     })))
}

q = allcoordinates(xx)
z = allcoordinates_lapply(xx)
all.equal(q,z)

So far it produces the same output. And now comes the nice part, doing a
test with a large polygon set and timing the performance of both functions:

 > system.time(x <- allcoordinates(nuts))
    user  system elapsed
  22.781   8.572  31.422
 > system.time(y <- allcoordinates_lapply(nuts))
    user  system elapsed
   0.248   0.004   0.250
 > all.equal(x,y)
[1] TRUE

So apart from the imo nicer looking function in terms of syntax, it is
also several orders of magnitude faster. This effect will only become
stronger if the polygon set becomes larger. So, lapply + do.call is your
friend :).

cheers,
Paul

 > sessionInfo()
R version 2.11.0 (2010-04-22)
i486-pc-linux-gnu

locale:
  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
  [5] LC_MONETARY=C              LC_MESSAGES=en_US.UTF-8
  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
  [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] maptools_0.7-34 lattice_0.18-5  sp_0.9-62       foreign_0.8-40

loaded via a namespace (and not attached):
[1] grid_2.11.0

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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Paul Hiemstra
On 06/11/2010 11:43 AM, Paul Hiemstra wrote:

> 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
>>>
>>> <[hidden email]>  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@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.
> In addition to my previous suggestion, I tried using lapply + do.call
> instead of the for loop based implementation of edzer:
>
> 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@polygons
>     for(i in 1:length(polys)) {
>         pp = polys[[i]]@Polygons
>         for (j in 1:length(pp))
>             ret = rbind(ret, coordinates(pp[[j]]))
>     }
>     ret
> }
>
> allcoordinates_lapply = function(x) {
>     polys = x@polygons
>     return(do.call("rbind", lapply(polys, function(pp) {
>     do.call("rbind", lapply(pp@Polygons, coordinates))
>     })))
> }
More information on why this is so much faster can be found in the R
Inferno Chapter 2, this gives a good description of what goes wrong in
the background when growing an R object.

http://www.burns-stat.com/pages/Tutor/R_inferno.pdf

The following basic R code also provides some info:

# Just letting the object grow
# without preallocation
meth1 = function(n) {
   vec <- numeric(0);
   for(i in 1:n) vec <- c(vec, i)
}

# preallocation
meth2 = function(n) {
   vec <- numeric(n)
   for(i in 1:n) vec[i] <- i
}

# Direct creation
# Not relevant for the
# above case
meth3 = function(n) {
   vec <- 1:n
}

# The suggestion I did
# with do.call and lapply
meth4 = function(n) {
   do.call("c", lapply(1:n, function(x) x))
}

n = 100000
# Speed test
system.time(meth1(n)) # Slooooow
system.time(meth2(n))
system.time(meth3(n))
system.time(meth4(n))

cheers,
Paul

>
> q = allcoordinates(xx)
> z = allcoordinates_lapply(xx)
> all.equal(q,z)
>
> So far it produces the same output. And now comes the nice part, doing
> a test with a large polygon set and timing the performance of both
> functions:
>
> > system.time(x <- allcoordinates(nuts))
>    user  system elapsed
>  22.781   8.572  31.422
> > system.time(y <- allcoordinates_lapply(nuts))
>    user  system elapsed
>   0.248   0.004   0.250
> > all.equal(x,y)
> [1] TRUE
>
> So apart from the imo nicer looking function in terms of syntax, it is
> also several orders of magnitude faster. This effect will only become
> stronger if the polygon set becomes larger. So, lapply + do.call is
> your friend :).
>
> cheers,
> Paul
>
> > sessionInfo()
> R version 2.11.0 (2010-04-22)
> i486-pc-linux-gnu
>
> locale:
>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=C              LC_MESSAGES=en_US.UTF-8
>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] maptools_0.7-34 lattice_0.18-5  sp_0.9-62       foreign_0.8-40
>
> loaded via a namespace (and not attached):
> [1] grid_2.11.0
>
>>
>> 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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

Thiago V. dos Santos
In reply to this post by edzer
  This function worked like a charm, but I can't individually invoke the columns which contain the coordinates of variable "q".
> names(q)NULL
  For example, I need to plot the coordinates using polygon function, whose help tells us that "‘polygon’ draws the polygons whose vertices are given in ‘x’ and ‘y’".
  So, I'd like to be able to refer to lon and lat columns of variable "q" as "q$lon" and "q$lat".
  Is it possible??
  Best wishes,
  Thiago.

--- On Wed, 9/6/10, Edzer Pebesma <[hidden email]> wrote:

From: Edzer Pebesma <[hidden email]>
Subject: Re: [R-sig-Geo] How to extract coordinates values from a shapefile?
To: [hidden email]
Date: Wednesday, 9 June, 2010, 17:24

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@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')


On 06/09/2010 09:58 PM, Matt Beard wrote:

>>
>> library(maptools)
>>
>> # A simple way to print out a list of coordinates for each polygon in your
>> shapefile:
>>
>> # Path and filename of polygon shapefile
>> testfile <- '/media/PKBACK# 001/FNR210/QGISLab/habitats/habitats.shp'
>>
>> # Read in polygon shapefile using handy maptools function
>> test <- readShapePoly(testfile)
>>
>> # Extract the list of Polygons objects
>> polys <- slot(test,"polygons")
>>
>> # Within each Polygons object
>> #    Extract the Polygon list (assuming just one per Polygons)
>> #    And Extract the coordinates from each Polygon
>> for (i in 1:length(polys)) {
>>    print(paste("Polygon #",i))
>>    print(slot(slot(polys[[i]],"Polygons")[[1]],"coords"  ))
>> }
>>
>>
>
>     [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



     
        [[alternative HTML version deleted]]


_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: How to extract coordinates values from a shapefile?

edzer
long = q[,1]
lat  = q[,2]

On 06/14/2010 05:05 AM, Thiago Veloso wrote:

>   This function worked like a charm, but I can't individually invoke the columns which contain the coordinates of variable "q".
>> names(q)NULL
>   For example, I need to plot the coordinates using polygon function, whose help tells us that "‘polygon’ draws the polygons whose vertices are given in ‘x’ and ‘y’".
>   So, I'd like to be able to refer to lon and lat columns of variable "q" as "q$lon" and "q$lat".
>   Is it possible??
>   Best wishes,
>   Thiago.
>
> --- On Wed, 9/6/10, Edzer Pebesma <[hidden email]> wrote:
>
> From: Edzer Pebesma <[hidden email]>
> Subject: Re: [R-sig-Geo] How to extract coordinates values from a shapefile?
> To: [hidden email]
> Date: Wednesday, 9 June, 2010, 17:24
>
> 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@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')
>
>
> On 06/09/2010 09:58 PM, Matt Beard wrote:
>>>
>>> library(maptools)
>>>
>>> # A simple way to print out a list of coordinates for each polygon in your
>>> shapefile:
>>>
>>> # Path and filename of polygon shapefile
>>> testfile <- '/media/PKBACK# 001/FNR210/QGISLab/habitats/habitats.shp'
>>>
>>> # Read in polygon shapefile using handy maptools function
>>> test <- readShapePoly(testfile)
>>>
>>> # Extract the list of Polygons objects
>>> polys <- slot(test,"polygons")
>>>
>>> # Within each Polygons object
>>> #    Extract the Polygon list (assuming just one per Polygons)
>>> #    And Extract the coordinates from each Polygon
>>> for (i in 1:length(polys)) {
>>>     print(paste("Polygon #",i))
>>>     print(slot(slot(polys[[i]],"Polygons")[[1]],"coords"  ))
>>> }
>>>
>>>
>>
>>     [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

--
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      [hidden email]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Loading...