[GIS] Data frame to SpatialPolygonsDataFrame with multiple polygons

polygon-creationr

I have a data frame with a bunch of polygons, several points per polygon. How can I convert it to a SpatialPolygonsDataFrame object?

Here's my data frame:

> buildings_df
        lon      lat               id
1  96.13597 16.83363 convocation_hall
2  96.13580 16.83315 convocation_hall
3  96.13631 16.83298 convocation_hall
4  96.13656 16.83351 convocation_hall
5  96.13400 16.83104    judson_church
6  96.13392 16.83073    judson_church
7  96.13432 16.83064    judson_church
8  96.13444 16.83085    judson_church
9  96.13596 16.83074      the_library
10 96.13586 16.83044      the_library
11 96.13615 16.83035      the_library
12 96.13626 16.83065      the_library
13 96.13599 16.83073      the_library

This can be reproduced with:

> dput(buildings_df)
structure(list(lon = c(96.135972, 96.1358, 96.136315, 96.136562, 
96.133997, 96.133922, 96.134319, 96.134437, 96.135961, 96.135864, 
96.136154, 96.136261, 96.135993), lat = c(16.833625, 16.833153, 
16.832978, 16.833512, 16.831038, 16.830729, 16.830637, 16.830853, 
16.83074, 16.830442, 16.83035, 16.830647, 16.830729), id = structure(c(1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L), .Label = c("convocation_hall", 
"judson_church", "the_library"), class = "factor")), .Names = c("lon", 
"lat", "id"), row.names = c(NA, -13L), class = "data.frame")

From this data.frame, I then convert it to a list, then to polygons, then to a SpatialPolygons object

# make a list
buildings_list <- split(buildings_df, buildings_df$id)
# only want lon-lats in the list, not the names
buildings_list <- lapply(buildings_list, function(x) { x["id"] <- NULL; x })

#  make data.frame into spatial polygon, cf. http://jwhollister.com/iale_open_science/2015/07/05/03-Spatial-Data-In-R/

# create SpatialPolygons Object, convert coords to polygon
ps <- sapply(buildings_list, Polygon)

# add id variable 
p1 <- Polygons(ps, ID = 1) 

# create SpatialPolygons object
my_spatial_polys <- SpatialPolygons(list(p1), proj4string = CRS("+proj=longlat +datum=WGS84") ) 

# let's see them
plot(my_spatial_polys)

So far so good:

And each of the polygons seems to be in there…

> str(my_spatial_polys)
Formal class 'SpatialPolygons' [package "sp"] with 4 slots
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 3
  .. .. .. .. ..$ convocation_hall:Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 96.1 16.8
  .. .. .. .. .. .. ..@ area   : num 3.08e-07
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
  .. .. .. .. ..$ judson_church   :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 96.1 43.3
  .. .. .. .. .. .. ..@ area   : num 0.0204
  .. .. .. .. .. .. ..@ hole   : logi TRUE
  .. .. .. .. .. .. ..@ ringDir: int -1
  .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
  .. .. .. .. ..$ the_library     :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 96.1 43.3
  .. .. .. .. .. .. ..@ area   : num 0.00512
  .. .. .. .. .. .. ..@ hole   : logi TRUE
  .. .. .. .. .. .. ..@ ringDir: int -1
  .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
  .. .. .. ..@ plotOrder: int [1:3] 1 2 3
  .. .. .. ..@ labpt    : num [1:2] 96.1 16.8
  .. .. .. ..@ ID       : chr "1"
  .. .. .. ..@ area     : num 0.0302
  ..@ plotOrder  : int 1
  ..@ bbox       : num [1:2, 1:2] 96.1 16.8 96.1 16.8
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

But they don't seem to be grouped in the way I'm expecting. When I try to set three ID variables, I get an error. How can I get the SpatialPolygons to have three polygons, each with its own ID?

Best Answer

One solution is to change p1 <- Polygons(ps, ID = 1) to an apply-like function so that the Polygons() function operates individually on each of the polygons in the ps object.

ps <- lapply(buildings_list, Polygon)

# add id variable
p1 <- lapply(seq_along(ps), function(i) Polygons(list(ps[[i]]), 
                                            ID = names(buildings_list)[i]  ))

# create SpatialPolygons object
my_spatial_polys <- SpatialPolygons(p1, proj4string = CRS("+proj=longlat +datum=WGS84") )

Then we can create the SpatialPolygonsDataFrame with the IDs as expected, one unique ID per polygon:

my_spatial_polys_df <- SpatialPolygonsDataFrame(my_spatial_polys, 
 data.frame(id = unique(buildings_df$id), 
            row.names = unique(buildings_df$id)))

> str(my_spatial_polys_df)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 3 obs. of  1 variable:
  .. ..$ id: Factor w/ 3 levels "convocation_hall",..: 1 2 3
  ..@ polygons   :List of 3
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 96.1 16.8
  .. .. .. .. .. .. ..@ area   : num 3.08e-07
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 96.1 16.8
  .. .. .. ..@ ID       : chr "convocation_hall"
  .. .. .. ..@ area     : num 0.0302
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 96.1 16.8
  .. .. .. .. .. .. ..@ area   : num 1.23e-07
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 96.1 16.8
  .. .. .. ..@ ID       : chr "judson_church"
  .. .. .. ..@ area     : num 0.0204
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 96.1 16.8
  .. .. .. .. .. .. ..@ area   : num 9.7e-08
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 96.1 16.8
  .. .. .. ..@ ID       : chr "the_library"
  .. .. .. ..@ area     : num 0.00512
  ..@ plotOrder  : int [1:3] 1 2 3
  ..@ bbox       : num [1:2, 1:2] 96.1 16.8 96.1 16.8
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

Related Solutions

[GIS] Extract all the polygon coordinates from a SpatialPolygonsDataframe

There are several ways to get at all of the coordinates for a SpatialPolygonsDataFrame or a SpatialLinesDataFrame.

Each of these methods records the object and the part within the object from the source, albeit with subtle differences and different names.

Say p is a SpatialPolygonsDataFrame.

#1. coerce to points, then to data frame
# names are Lines.NR Lines.ID Line.NR coords.x1 coords.x2
## start here for polys
lin <- as(p, "SpatialLinesDataFrame")  
## start here for lines
pts <- as.data.frame(as(lin, "SpatialPointsDataFrame"))


#2. convert to matrix 
# names are      object part cump hole         x        y
mpts <- raster::geom(p)

#3. convert to data frame
# names are        long      lat order  hole piece  id group
dpts <- ggplot2::fortify(p)

#4. convert to data frame
## names are   object_ branch_ island_ order_        x_       y_
spts <- spbabel::sptable(p)

Each of these outputs can be used interchangeably as long as you deal with the different name conventions, and for raster::geom remember that it's a matrix. (Actually the coercion one doesn't record the hole/island status so that is less useful depending on what's needed).

You can round trip the last one by

p2 <- spbabel::sp(spts, attr_tab = as.data.frame(p))

(There's an added bonus with spbabel in that you get a 'raster-like' print for the object, but with the data in 'tibble' print form. )

[GIS] Extract values from raster using point data with buffer in R

If buffer is in meters the code works, but if I use degrees, it won't.

Indeed; the buffer cannot be specified in degrees. ?extract, buffer says: "If the data are not projected (latitude/longitude), the unit should be meters. Otherwise it should be in map-units (typically also meters)."

Do you really want the buffer to be in degrees?

Is so, temporarily set the CRS to UTM, do the extract, and set it back.

cr <- crs(r)
crs(r) <- '+proj=utm +zone=10 +datum=WGS84'
e <- extract(r, xy, buffer=0.5)
crs(r) <- cr

Best Answer

Related Solutions

[GIS] Extract all the polygon coordinates from a SpatialPolygonsDataframe

[GIS] Extract values from raster using point data with buffer in R

Related Question