Converting GeoJSON into a Simple Feature in R

Question

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.

Answer 1

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.