Skip to content

Retrieving Vertices Coordinates from SpatialPolygons

5 messages · Enrico R. Crema, Barry Rowlingson, Hadley Wickham +1 more

David Winsemius
# 
If it is in an S4 slot then there should be an extraction function,  
possibly with an obvious name like coords(). Failing that you could  
offer the name of the package that you are using. Perhaps that would  
contain an example to be used for illustration.
Barry Rowlingson
#
On Sat, Mar 21, 2009 at 12:33 PM, Enrico R. Crema <e.crema at ucl.ac.uk> wrote:
There's not really such thing as a SpatialPolygon - there's only
SpatialPolygons - so maybe you've got a SpatialPolygons object with
only one feature in it..

 str(thing) will help here. This is what I get if I str() something
that is a SpatialPolygons object with only one feature in it:
Formal class 'SpatialPolygons' [package "sp"] with 4 slots
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 1 10
  .. .. .. .. .. .. ..@ area   : num 1
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 0.5 0.5 1.5 1.5 0.5
9.5 10.5 10.5 9.5 9.5
  .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. .. ..$ : chr [1:5] "s1" "s1" "s1" "s1" ...
  .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 1 10
  .. .. .. ..@ ID       : chr "g1"
  .. .. .. ..@ area     : num 1
  ..@ plotOrder  : int 1
  ..@ bbox       : num [1:2, 1:2] 0.5 9.5 1.5 10.5
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "r1" "r2"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr NA

 I'll explain, stepping through the structure breadth first, and backwards:

 @proj4string is the coordinate reference system slot
 @bbox is the bounding box slot
 @plotOrder is the order to plot the polygons
 @polygons is the list of Polygons objects.

 @polygons[[1]] is the first (and in this case, only) feature. It is
an object of class 'Polygons' (non-spatial, since there's no
@proj4string in this part of the structure).

 @polygons[[1]]@Polygons is the 'Polygons' slot of class 'Polygons',
and is a list of rings that make up the feature.

 @polygons[[1]]@Polygons[[1]] is an object of class 'Polygon'.

 @polygons[[1]]@Polygons[[1]]@coords is the coordinates of the Polygon:
x    y
s1 0.5  9.5
s1 0.5 10.5
s1 1.5 10.5
s1 1.5  9.5
s1 0.5  9.5

 It's complex because the SpatialPolygons class can hold multiple
features (such as administrative subdivisions), and each feature can
be made from multiple rings (islands, holes). It seems overkill when
all you really want is to store a single simple closed polygon though!
But now you know how to get it. There may be methods for some of these
@slot calls.

 Personally, I forget this everytime, but str() is immensely useful.
As is a copy of the Bivand/Pebesma/Gomez-Rubio book!

Barry
Hadley Wickham
# 
This came up on R-sig-geo two days ago and this is what I said:

I have the following code in ggplot2 for turning a SpatialPolygon into
a regular data frame of coordinates.  You'll need to load ggplot2, and
then run fortify(yoursp).

fortify.SpatialPolygonsDataFrame <- function(shape, region = NULL) {
 attr <- as.data.frame(shape)
 # If not specified, split into regions based on first variable in attributes
 if (is.null(region)) {
   region <- names(attr)[1]
   message("Using ", region, " to define regions.")
 }

 # Figure out how polygons should be split up into the region of interest
 polys <- split(as.numeric(row.names(attr)), attr[, region])
 cp <- polygons(shape)

 # Union together all polygons that make up a region
 try_require(c("gpclib", "maptools"))
 unioned <- unionSpatialPolygons(cp, invert(polys))

 coords <- fortify(unioned)
 coords$order <- 1:nrow(coords)
 coords
}

fortify.SpatialPolygons <- function(model, data, ...) {
 ldply(model at polygons, fortify)
}

fortify.Polygons <- function(model, data, ...) {
 subpolys <- model at Polygons
 pieces <- ldply(seq_along(subpolys), function(i) {
   df <- fortify(subpolys[[model at plotOrder[i]]])
   df$piece <- i
   df
 })

 within(pieces,{
   order <- 1:nrow(pieces)
   id <- model at ID
   piece <- factor(piece)
   group <- interaction(id, piece)
 })
}

fortify.Polygon <- function(model, data, ...) {
 df <- as.data.frame(model at coords)
 names(df) <- c("long", "lat")
 df$order <- 1:nrow(df)
 df$hole <- model at hole
 df
}

fortify.SpatialLinesDataFrame <- function(model, data, ...) {
 ldply(model at lines, fortify)
}

fortify.Lines <- function(model, data, ...) {
 lines <- model at Lines
 pieces <- ldply(seq_along(lines), function(i) {
   df <- fortify(lines[[i]])
   df$piece <- i
   df
 })

 within(pieces,{
   order <- 1:nrow(pieces)
   id <- model at ID
   piece <- factor(piece)
   group <- interaction(id, piece)
 })
}

fortify.Line <- function(model, data, ...) {
 df <- as.data.frame(model at coords)
 names(df) <- c("long", "lat")
 df$order <- 1:nrow(df)
 df
}


Hadley
On Sat, Mar 21, 2009 at 7:33 AM, Enrico R. Crema <e.crema at ucl.ac.uk> wrote:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      David Winsemius
    

    

      
      #
    

  


  

      
On Mar 21, 2009, at 10:31 AM, Barry Rowlingson wrote:

        

            

      
      
      
    

        
The last example on the following page suggests to me that there is a  
coordinates function for this class, although I do not see it  
mentioned in the body of the class description:

http://finzi.psych.upenn.edu/R/library/sp/html/SpatialPolygons-class.html

            

      
      
      
    

        
David Winsemius, MD
Heritage Laboratories
West Hartford, CT