I would like to read a spatial dataset from a Socrata repository into R and then convert it into a simple feature object.
The dataset consists of capital improvement projects represented as polygons:
The data used in the app above is accessible via a Socrata Open Data API (SODA):
library(tibble)
library(dplyr)
library(purrr)
library(sf)
library(RSocrata)
obj <- read.socrata("https://data.seattle.gov/resource/pdbw-sw7q.json") %>% as_tibble()
And the spatial data appears to be in the the_geom.coordinates
column:
# Inspect the object
glimpse(obj)
#> Observations: 113
#> Variables: 13
#> $ creationdate <chr> "2017-06-02T13:51:50.854Z", "2017-06-02T1...
#> $ creator <chr> "Transportation_SeattleCityGIS", "Transpo...
#> $ editdate <chr> "2017-06-02T13:51:50.854Z", "2017-06-02T1...
#> $ editor <chr> "Transportation_SeattleCityGIS", "Transpo...
#> $ globalid <chr> "4a78a16a-9ea4-4a81-8011-ad974c80b357", "...
#> $ objectid <chr> "967", "646", "968", "11862", "11521", "1...
#> $ project_id <chr> "TC36717008", "TC367240", "TC36659003", "...
#> $ projectname <chr> "Safe Routes to School - S Fisher Place S...
#> $ shape_area <chr> "4.53868971176276E-8", "0.000002901627518...
#> $ shape_length <chr> "0.000918371270091608", "0.02024322483978...
#> $ status <chr> "ACTIVE", "ACTIVE", "ACTIVE", "ACTIVE", "...
#> $ the_geom.type <chr> "Polygon", "Polygon", "Polygon", "Polygon...
#> $ the_geom.coordinates <list> [<-122.26983, -122.26984, -122.27015, -1... <-- here
Taking a look at the last column reveals that each polygon is stored as an array (or a list of arrays for multipolygons):
# Inspect the spatial data
obj %>% select(the_geom.type, the_geom.coordinates) %>%
mutate(class = map_chr(the_geom.coordinates, class))
#> # A tibble: 113 x 3
#> the_geom.type the_geom.coordinates class
#> <chr> <list> <chr>
#> 1 Polygon <dbl [1 x 5 x 2]> array
#> 2 Polygon <dbl [1 x 16 x 2]> array
#> 3 Polygon <dbl [1 x 5 x 2]> array
#> 4 Polygon <dbl [1 x 35 x 2]> array
#> 5 Polygon <dbl [1 x 24 x 2]> array
#> 6 Polygon <dbl [1 x 15 x 2]> array
#> 7 Polygon <list [2]> list
#> 8 <NA> <NULL> NULL
#> 9 Polygon <list [2]> list
#> 10 Polygon <dbl [1 x 10 x 2]> array
#> # ... with 103 more rows
obj %>% slice(1) %>% pull
#> [[1]]
#> , , 1
#>
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] -122.2698 -122.2698 -122.2702 -122.2702 -122.2698
#>
#> , , 2
#>
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 47.52145 47.5213 47.52131 47.52145 47.52145
I have been unable to transform these arrays into polygons using the tools provided by the sf
package:
# Try to convert one row from the `the_geom.coordinates` column
# into a POYLGON or MULTIPOLYGON
obj[1, "the_geom.coordinates"] %>% st_polygon
#> Error in vapply(x, ncol, 0L): values must be length 1,
#> but FUN(X[[1]]) result is length 0
obj[1, "the_geom.coordinates"] %>% st_multipolygon
#> Error in MtrxSetSet(x, dim, type = "MULTIPOLYGON", needClosed = TRUE):
#> polygons not (all) closed
Any advice on how to convert obj
into an sf
object would be much appreciated.
You've not got geoJSON. That URL seems to get a JSON download with the geometry encoded in some way.
With the latest sf
package, you can do:
> d = read_sf("https://data.seattle.gov/resource/pdbw-sw7q.geojson")
>
However if you need to add RSocrata API keys etc, or you are requesting a large number and need to batch it (ie get 1000 at a time) then you'll have to do that manually. RSocrata will attempt to get batched data with multiple requests.
There is an open request on the RSocrata github site for geoJSON functionality: https://github.com/Chicago/RSocrata/issues/43 but it didn't seem very well received.