Manipulating polygon data imported from Google BigQuery GIS within R

Question

I have some polygon shapefiles (representing UK postcodes) which I have uploaded to BigQuery, using this guide, hoping to take advantage of the inbuilt GIS support described here.

I am trying to import these into R, and from there plot these in an R Shiny app (or use for other analysis).

Whilst I have managed to import the polygon data into R, using the "bigrquery" package, the data extracted is in an unfamiliar format to me - a table where the first column is a GeoJSON string for each shape (containing only the geography), and the other columns are attributes:

# A tibble: 5 x 4
  geom                                                            POSTCODE   UPP PC_AREA
  <chr>                                                           <chr>    <int> <chr>  
1 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.095015148797922,5… AB10 1AL     3 AB     
2 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.095031775835887,5… AB10 1AN     4 AB     
3 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.096682747887519,5… AB10 1AP     5 AB     
4 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.09786728657654,57… AB10 1AS     6 AB     
5 "{\"type\":\"Polygon\",\"coordinates\":[[[-2.095637704208373,5… AB10 1AU     7 AB     

I am trying to convert this within R to something more familiar e.g. an sf object or a SpatialPolygonsDataFrame, which will be easier to further manipulate/plot etc.

Is there a function I can use to do this?

Answer 1

If you row-bind together the list that you get when reading GeoJSON feature strings with read_sf together you get an sf data frame:

> do.call(rbind, lapply(d$geom,read_sf))
Simple feature collection with 3 features and 0 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
                        geometry
1 MULTIPOLYGON (((-81.47276 3...
2 MULTIPOLYGON (((-81.47276 3...
3 MULTIPOLYGON (((-81.47276 3...

That's the crux of the process, the fact that st_read and read_sf can eat a GeoJSON feature string and return an sf object. Here's one way of adding that to the original data frame and making it spatial:

> ds = d # work on a copy

Get the geometry of all the geom elements:

> ds$geom = do.call(rbind, lapply(d$geom,read_sf))$geometry

This isn't spatial yet:

> ds
                            geom ID         Z
1 MULTIPOLYGON (((-81.47276 3...  A 0.3324811
2 MULTIPOLYGON (((-81.47276 3...  B 0.8053311
3 MULTIPOLYGON (((-81.47276 3...  C 0.7607013
> class(ds)
[1] "data.frame"

So tell it where its geometry is:

> st_geometry(ds)=ds$geom
> ds
Simple feature collection with 3 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  ID         Z                           geom
1  A 0.3324811 MULTIPOLYGON (((-81.47276 3...
2  B 0.8053311 MULTIPOLYGON (((-81.47276 3...
3  C 0.7607013 MULTIPOLYGON (((-81.47276 3...