Sf: easy conversion for GeoJSON as R list to sf

Created on 20 Jan 2017  路  41Comments  路  Source: r-spatial/sf

As discussed in https://github.com/edzer/sfr/issues/39#issuecomment-273908034, easily converting GeoJSON as an R list to sf would be really helpful. @ateucher has already done all the hard work for the reverse sf -> GeoJSON https://github.com/ropensci/geojsonio/pull/95.

Most helpful comment

We now have

> geo <- c("{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68152563269095,36.43764870908927]}",
+          "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67408758213843,36.43366018922779]}",
+          "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67708346361097,36.44208638659282]}",
+          "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67886661944996,36.44110273135671]}",
+          "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68089232041565,36.44173155205561]}")
> st_as_sfc(geo, GeoJSON=TRUE)
Geometry set for 5 features 
geometry type:  POINT
dimension:      XY
bbox:           xmin: -118.6815 ymin: 36.43366 xmax: -118.6741 ymax: 36.44209
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POINT (-118.681525632691 36.4376487090893)
POINT (-118.674087582138 36.4336601892278)
POINT (-118.677083463611 36.4420863865928)
POINT (-118.67886661945 36.4411027313567)
POINT (-118.680892320416 36.4417315520556)

All 41 comments

If someone can provide me with examples of each of the seven types, as R data, I'll take a look.

Would it be safe to assume that character data is GeoJSON? I know geojsonio::geojson_json uses classes json and geo_json.

@edzer, fortunately with the work done by @ateucher, this should be easy. I'll try to do at least a couple quick examples unless someone else wants to :)

sf assumes that character data is WKT; I may hope the R lists with GeoJSON have a class?

@edzer, this will provide a pretty thorough start using your examples in sf.

library(sf)
# use newest geojsonio from Github
# devtools::install_github("ropensci/geojsonio")
library(geojsonio)
library(magrittr)

# use sf examples
p1 = st_point(c(1,2))
geojsonio::geojson_list(p1)

(p2 = st_point(c(1,2,3)))
geojsonio::geojson_list(p2)

p3 = st_point(c(1,2,3), "XYM")
geojsonio::geojson_list(p3)

pts = matrix(1:10, , 2)
(mp1 = st_multipoint(pts))
geojsonio::geojson_list(mp1)

pts = matrix(1:15, , 3)
(mp2 = st_multipoint(pts))
geojsonio::geojson_list(mp2)

(mp3 = st_multipoint(pts, "XYM"))
geojsonio::geojson_list(mp3)

pts = matrix(1:20, , 4)
(mp4 = st_multipoint(pts))
geojsonio::geojson_list(mp4)

pts = matrix(1:10, , 2)
(ls1 = st_linestring(pts))
geojsonio::geojson_list(ls1)

pts = matrix(1:15, , 3)
(ls2 = st_linestring(pts))
geojsonio::geojson_list(ls2)

(ls3 = st_linestring(pts, "XYM"))
pts = matrix(1:20, , 4)
geojsonio::geojson_list(ls3)

(ls4 = st_linestring(pts))
geojsonio::geojson_list(ls4)


outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE)
hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE)
hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE)
pts = list(outer, hole1, hole2)
(ml1 = st_multilinestring(pts))
geojsonio::geojson_list(ml1)

pts3 = lapply(pts, function(x) cbind(x, 0))
(ml2 = st_multilinestring(pts3))
geojsonio::geojson_list(ml2)

(ml3 = st_multilinestring(pts3, "XYM"))
geojsonio::geojson_list(ml3)

pts4 = lapply(pts3, function(x) cbind(x, 0))
(ml4 = st_multilinestring(pts4))
geojsonio::geojson_list(ml4)

outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE)
hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE)
hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE)
pts = list(outer, hole1, hole2)
(pl1 = st_polygon(pts))
geojsonio::geojson_list(pl1)

pts3 = lapply(pts, function(x) cbind(x, 0))
(pl2 = st_polygon(pts3))
geojsonio::geojson_list(pl2)

(pl3 = st_polygon(pts3, "XYM"))
geojsonio::geojson_list(pl3)

pts4 = lapply(pts3, function(x) cbind(x, 0))
(pl4 = st_polygon(pts4))
geojsonio::geojson_list(pl4)

pol1 = list(outer, hole1, hole2)
pol2 = list(outer + 12, hole1 + 12)
pol3 = list(outer + 24)
mp = list(pol1,pol2,pol3)
(mp1 = st_multipolygon(mp))
geojsonio::geojson_list(mp1)

pts3 = lapply(mp, function(x) lapply(x, function(y) cbind(y, 0)))
(mp2 = st_multipolygon(pts3))
geojsonio::geojson_list(mp2)

(mp3 = st_multipolygon(pts3, "XYM"))
geojsonio::geojson_list(mp3)

pts4 = lapply(mp2, function(x) lapply(x, function(y) cbind(y, 0)))
(mp4 = st_multipolygon(pts4))
geojsonio::geojson_list(mp4)

(gc = st_geometrycollection(list(p1, ls1, pl1, mp1)))
geojsonio::geojson_list(gc)
st_geometrycollection() # empty geometry
c(st_point(1:2), st_point(5:6)) %>%
  geojsonio::geojson_list(mp4)

c(st_point(1:2), st_multipoint(matrix(5:8,2))) %>%
  geojsonio::geojson_list()

c(st_multipoint(matrix(1:4,2)), st_multipoint(matrix(5:8,2))) %>%
  geojsonio::geojson_list()
c(st_linestring(matrix(1:6,3)), st_linestring(matrix(11:16,3))) %>%
  geojsonio::geojson_list()

c(st_multilinestring(list(matrix(1:6,3))), st_multilinestring(list(matrix(11:16,3)))) %>%
  geojsonio::geojson_list()

pl = list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0)))
c(st_polygon(pl), st_polygon(pl)) %>%
  geojsonio::geojson_list()
c(st_polygon(pl), st_multipolygon(list(pl))) %>%
  geojsonio::geojson_list()

c(st_linestring(matrix(1:6,3)), st_point(1:2)) %>%
  geojsonio::geojson_list()

c(st_geometrycollection(list(st_point(1:2), st_linestring(matrix(1:6,3)))),
  st_geometrycollection(list(st_multilinestring(list(matrix(11:16,3)))))) %>%
  geojsonio::geojson_list()

c(st_geometrycollection(list(st_point(1:2), st_linestring(matrix(1:6,3)))),
  st_multilinestring(list(matrix(11:16,3))), st_point(5:6), 
  st_geometrycollection(list(st_point(10:11)))) %>%
  geojsonio::geojson_list()

Here's a a first shot at it. It returns sfc with one feature, as these have a crs. Alternatively, one could return single simple feature geometries, without crs (which must be 4326 anyway).

library(sf)
# use newest geojsonio from Github
# devtools::install_github("ropensci/geojsonio")
library(geojsonio)
library(magrittr)
library(testthat)

st_as_sfc.geo_list = function(x, ...) {
    x = switch(x$type,
        Point = st_point(x$coordinates),
        MultiPoint = st_multipoint(x$coordinates),
        LineString = st_linestring(x$coordinates),
        MultiLineString = st_multilinestring(x$coordinates),
        Polygon = st_polygon(x$coordinates),
        MultiPolygon = st_multipolygon(x$coordinates),
        GeometryCollection = st_geometrycollection(
            lapply(x$geometries, function(y) st_as_sfc.geo_list(y)[[1]])),
        stop("unkown class")
    )
    st_sfc(x, crs = st_crs(4326))
}

# use sf examples
p1 = st_point(c(1,2))
expect_equal(st_as_sfc(geojsonio::geojson_list(p1))[[1]], p1)

(p2 = st_point(c(1,2,3)))
expect_equal(st_as_sfc(geojsonio::geojson_list(p2))[[1]], p2)

pts = matrix(1:10, , 2)
(mp1 = st_multipoint(pts))
expect_equal(st_as_sfc(geojsonio::geojson_list(mp1))[[1]], mp1)

pts = matrix(1:15, , 3)
(mp2 = st_multipoint(pts))
expect_equal(st_as_sfc(geojsonio::geojson_list(mp2))[[1]], mp2)

pts = matrix(1:10, , 2)
(ls1 = st_linestring(pts))
expect_equal(st_as_sfc(geojsonio::geojson_list(ls1))[[1]], ls1)

pts = matrix(1:15, , 3)
(ls2 = st_linestring(pts))
expect_equal(st_as_sfc(geojsonio::geojson_list(ls2))[[1]], ls2)

outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE)
hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE)
hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE)
pts = list(outer, hole1, hole2)
(ml1 = st_multilinestring(pts))
expect_equal(st_as_sfc(geojsonio::geojson_list(ml1))[[1]], ml1)

pts3 = lapply(pts, function(x) cbind(x, 0))
(ml2 = st_multilinestring(pts3))
expect_equal(st_as_sfc(geojsonio::geojson_list(ml2))[[1]], ml2)

outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE)
hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE)
hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE)
pts = list(outer, hole1, hole2)
(pl1 = st_polygon(pts))
geojsonio::geojson_list(pl1)
st_as_sfc(geojsonio::geojson_list(pl1))[[1]]
expect_equal(st_as_sfc(geojsonio::geojson_list(pl1))[[1]], pl1)

pts3 = lapply(pts, function(x) cbind(x, 0))
(pl2 = st_polygon(pts3))
expect_equal(st_as_sfc(geojsonio::geojson_list(pl2))[[1]], pl2)

pol1 = list(outer, hole1, hole2)
pol2 = list(outer + 12, hole1 + 12)
pol3 = list(outer + 24)
mp = list(pol1,pol2,pol3)
(mp1 = st_multipolygon(mp))
expect_equal(st_as_sfc(geojsonio::geojson_list(mp1))[[1]], mp1)

pts3 = lapply(mp, function(x) lapply(x, function(y) cbind(y, 0)))
(mp2 = st_multipolygon(pts3))
expect_equal(st_as_sfc(geojsonio::geojson_list(mp2))[[1]], mp2)

(gc = st_geometrycollection(list(p1, ls1, pl1, mp1)))
# x = geojsonio::geojson_list(gc)
# lapply(x$geometries, st_as_sfc.geo_list)
expect_equal(st_as_sfc(geojsonio::geojson_list(gc))[[1]], gc)

gc = st_geometrycollection() # empty geometry
expect_equal(st_as_sfc(geojsonio::geojson_list(gc))[[1]], gc)

gc = st_geometrycollection(list(st_point(1:2), st_multipoint(matrix(5:8,2))) )
expect_equal(st_as_sfc(geojsonio::geojson_list(gc))[[1]], gc)

gc = st_geometrycollection(
    list(st_multilinestring(list(matrix(1:6,3))), 
    st_multilinestring(list(matrix(11:16,3)))))
expect_equal(st_as_sfc(geojsonio::geojson_list(gc))[[1]], gc)

Did this work out, and should we add it to sf?

Did anyone look at this?

@edzer I'll have a look at these soon.
@timelyportfolio did you check these? I think these would be crucial for mapedit.

I should be able to test tomorrow. Sorry for the delay. Yes, these will be very important for mapedit.

  • [ ] What about FEATURE and FEATURECOLLECTION? I'll try to work through some code to handle these.

  • [ ] In some geojson, I get coordinates as a list which causes the conversion to fail.

    structure(list(type = "Point", coordinates = list(-78.0080795288, 
        39.9357821415, 0L)), .Names = c("type", "coordinates"), class = "geo_list")
    
  • [ ] Work through geojsonio tests in reverse https://github.com/ropensci/geojsonio/blob/master/tests/testthat/test-sf_classes.R

I guess a FEATURECOLLECTION would map into an sf object, where non-geometry properties go into the non-geometry columns, and a FEATURE is a single row (record) of such a table.

AFAICS, coordinates in geojson should be a vector or array. Who created this geo_list?

Thanks @edzer. By no means did I intend to imply that coordinates as a list are correct (I have no idea :)). I just discovered as I randomly pulled every geojson example I could find. Here is the code.

library(geojsonio)

file <- system.file("examples", "norway_maple.kml", package = "geojsonio")
nwy_map <- geojson_read(as.location(file), what = "list")
listviewer::jsonedit(nwy_map)

pt <- nwy_map$features[[1]]$geometry
class(pt) <- "geo_list"

# test geojson conversion functions
library(sf)
st_as_sfc.geo_list = function(x, ...) {
  x = switch(x$type,
             Point = st_point(x$coordinates),
             MultiPoint = st_multipoint(x$coordinates),
             LineString = st_linestring(x$coordinates),
             MultiLineString = st_multilinestring(x$coordinates),
             Polygon = st_polygon(x$coordinates),
             MultiPolygon = st_multipolygon(x$coordinates),
             GeometryCollection = st_geometrycollection(
               lapply(x$geometries, function(y) st_as_sfc.geo_list(y)[[1]])),
             stop("unkown class")
  )
  st_sfc(x, crs = st_crs(4326))
}

I think this geojsonio test offers a good basis for FEATURECOLLECTION.

p_list <- lapply(list(c(3.2,4), c(3,4.6), c(3.8,4.4)), st_point)
pt_sfc <- st_sfc(p_list)
pt_sf <- st_sf(x = c("a", "b", "c"), pt_sfc)

st_as_sfc(geojson_list(pt_sfc))
st_as_sfc(geojson_list(pt_sf))

Do the geojsonio folks somewhere specify how they want to represent GeoJSON in R?

pinging @sckott and @ateucher. In the geojson_list tests, I only see coordinates as a list in this line. For reference, jsonlite treats an unnamed list and a vector equivalently as an array in JavaScript if auto_unbox=TRUE.

> jsonlite::toJSON(list(1,2,3))
[[1],[2],[3]] 
> jsonlite::toJSON(list(1,2,3), auto_unbox=TRUE)
[1,2,3] 
> jsonlite::toJSON(c(1,2,3), auto_unbox=TRUE)
[1,2,3] 
> jsonlite::toJSON(c(1,2,3))
[1,2,3] 

@timelyportfolio i think geojson as a list or string both make sense - i'd like to soon depend on our geojson pkg to have a unified geojson representation with simple s3 print methods, see https://github.com/ropensci/geojson - the workflow and print methods in geojson only handle strings right now - but may add methods that handle lists in geojson - then we'd have similar s3 print methods for geojson as lists and strings

@sckott correct me if I'm wrong, but because jsonlite is so clever and flexible (can just as easily handle coordinates as lists or vectors/matrices or combos thereof), we've been pretty lax about enforcing a strict structure on the geo_list. Probably something we should remedy...

@ateucher right, we should specify a structure for that type - i thin kit would make most sense to include list support in geojson - and have checker fxns and docs there for proper format - then that can be used elsewhere

Since your geometries are a constrained subset of simple features, you could as well follow the sf pattern (if you don't already): point is a vector, linestring/multipoint is matrix (with each point a row), multilinestring/polygon are a list of matrices, multipolygon is a list of list of matrices, geometrycollection is a list of any of the previous.

I think that's a great idea @edzer

This is some ugly code, but I'll post as a starting point to converting FEATURECOLLECTION and FEATURE.

# test geojson conversion functions
library(sf)
st_as_sfc.geo_list = function(x, ...) {
  x = switch(x$type,
             Point = st_point(x$coordinates),
             MultiPoint = st_multipoint(x$coordinates),
             LineString = st_linestring(x$coordinates),
             MultiLineString = st_multilinestring(x$coordinates),
             Polygon = st_polygon(x$coordinates),
             MultiPolygon = st_multipolygon(x$coordinates),
             GeometryCollection = st_geometrycollection(
               lapply(x$geometries, function(y) st_as_sfc.geo_list(y)[[1]])),
             stop("unkown class")
  )
  st_sfc(x, crs = st_crs(4326))
}

p_list <- lapply(list(c(3.2,4), c(3,4.6), c(3.8,4.4)), st_point)
pt_sfc <- st_sfc(p_list)
pt_sf <- st_sf(x = c("a", "b", "c"), pt_sfc)

# use geojsonio to convert sf to geojson
gj_pt_sf <- geojsonio::geojson_list(pt_sf)

# start by trying to convert a single geojson feature to s
feat <- gj_pt_sf$features[[1]]

# manually translate one geojson feature
st_sf(
  feat$properties,
  st_as_sfc.geo_list(feat$geometry)
)


# now try to convert all of the geojson features
feats <- lapply(
  gj_pt_sf$features,
  function(ft) {
    st_sf(
      ft$properties,
      st_as_sfc.geo_list(ft$geometry)
    )
  }
)

do.call(sf::rbind.sf, feats)

I'm not very familiar with json, but it seems that multiple R representations could be valid?

Example with OSM data:

response = httr::GET("https://nominatim.openstreetmap.org/search?format=json&limit=1&polygon_geojson=1&q=hawai")

js_h <- httr::content(response, "parsed")[[1]]
js_l <- jsonlite::fromJSON(httr::content(response, as="text"))

str(js_h$geojson, 2)
#> List of 2
#> $ type       : chr "MultiPolygon"
#> $ coordinates:List of 17
#> ..$ :List of 1
#> ..$ :List of 1
#> ...
str(js_l$geojson, 3)
#> 'data.frame':    1 obs. of  2 variables:
#>   $ type       : chr "MultiPolygon"
#> $ coordinates:List of 1
#> ..$ :List of 17
#> .. ..$ : num [1, 1:248, 1:2] -157 -157 -157 -157 -157 ...
#> .. ..$ : num [1, 1:235, 1:2] -157 -157 -157 -157 -157 ...
#> .. ..$ : num [1, 1:192, 1:2] -171 -171 -171 -171 -171 ...
#> .. ..$ : num [1, 1:156, 1:2] -166 -166 -166 -166 -166 ...
#> .. ..$ : num [1, 1:134, 1:2] -156 -156 -156 -156 -156 ...
#> ...

I couldn't get either to work right of the bat with st_as_sfc.geo_list. As I said, I'm not very familiar with json, so is there a preferred way to cast it to an R object? R Nominatim package prefers the second approach.

Looks like you found one I have not discovered yet. Here is my first attempt at converting coordinates from json to the sf-expected https://github.com/timelyportfolio/mapedit/blob/3cb976875799ba1312b04bddc5d0fa959e5b8ab5/R/edit_map_return_sf.R#L43-L68.

Any reason to keep this issue open?

@edzer even though I opened the issue, I am probably the least qualified to answer "open or close?" Whether easy geojson -> sf belongs here or somewhere else like geojsonio, I am not sure, but I do think that this is an incredibly valuable gap missing currently in sf. While geojsonio does sf -> geojson, nothing fills this geojson->sf gap currently. Since it is not technically difficult with a little testing and a lot of work has already been done, I'd love to know that somewhere it gets addressed. @sckott, @ateucher, @tim-salabim any additional thoughts?

One concern I do have is that sf assumes character is WKT, and often geojson is just text with no additional class clues, so this might be difficult to plug into S3.

st_read can read geojson from text files, does that help?

If geojson comes as a text string, why wouldn't it have a class?

for wkt, sf has an st_as_sfc.character method, but there is no st_as_sf.character - we could use that for geojson?

st_read can actually read geojson from a character string, as can readOGR. I don't think it's documented anywhere but it works:

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3

geojson_txt <- "{\"type\":\"MultiPoint\",\"coordinates\":[[3.2,4],[3,4.6],[3.8,4.4],[3.5,3.8],[3.4,3.6],[3.9,4.5]]}"

st_read(geojson_txt)
#> Reading layer `OGRGeoJSON' from data source `{"type":"MultiPoint","coordinates":[[3.2,4],[3,4.6],[3.8,4.4],[3.5,3.8],[3.4,3.6],[3.9,4.5]]}' using driver `GeoJSON'
#> converted into: MULTIPOINT
#> Simple feature collection with 1 feature and 0 fields
#> geometry type:  MULTIPOINT
#> dimension:      XY
#> bbox:           xmin: 3 ymin: 3.6 xmax: 3.9 ymax: 4.6
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs

Super cool, didn't know that!

For this case, it pays off to use read_sf instead, which won't print the data source:

> x = read_sf(geojson_txt)
> x
Simple feature collection with 1 feature and 0 fields
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: 3 ymin: 3.6 xmax: 3.9 ymax: 4.6
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
                        geometry
1 MULTIPOINT(3.2 4, 3 4.6, 3....

at @timelyportfolio : I guess the geojson -> sf is now sufficiently answered, and only in need of better documentation?

Guess so :) Thanks @ateucher for answering from Runconf! Now I need to see if I should do this instead of my hacked conversion in mapedit.

What if the geoJSON character string is a column in a data.frame? I can read a single cell using st_read() or read_sf() as described above, but ultimately I'm trying to turn this data.frame into a sf object by converting the character column of geoJSONs to an sfc, which should then let me convert the data.frame with st_as_sf()?

Can the mechanics of read_sf() that allow it to handle character string geoJSON be used in st_as_sfc()?

Yes; a reprex would be helpful here! I'd try to unlist it, with recursive = FALSE. Maybe twice.

Okay! We know that read_sf() and st_read() can parse a character representation of a geoJSON:

library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3

geo <- "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68152563269095,36.43764870908927]}"

read_sf(geo)
#> Simple feature collection with 1 feature and 0 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: -118.6815 ymin: 36.43765 xmax: -118.6815 ymax: 36.43765
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
#>                         geometry
#> 1 POINT (-118.681525632691 36...

But if we have a data.frame with a column that is a character representation of a geoJSON...

geo <- c("{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68152563269095,36.43764870908927]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67408758213843,36.43366018922779]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67708346361097,36.44208638659282]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67886661944996,36.44110273135671]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68089232041565,36.44173155205561]}")

test <- data.frame(numbers = 1:5, letters = letters[1:5], geometry = geo)

test
#>   numbers letters
#> 1       1       a
#> 2       2       b
#> 3       3       c
#> 4       4       d
#> 5       5       e
#>                                                                                 geometry
#> 1 {"geodesic":true,"type":"Point","coordinates":[-118.68152563269095,36.43764870908927]}
#> 2 {"geodesic":true,"type":"Point","coordinates":[-118.67408758213843,36.43366018922779]}
#> 3 {"geodesic":true,"type":"Point","coordinates":[-118.67708346361097,36.44208638659282]}
#> 4 {"geodesic":true,"type":"Point","coordinates":[-118.67886661944996,36.44110273135671]}
#> 5 {"geodesic":true,"type":"Point","coordinates":[-118.68089232041565,36.44173155205561]}

...we might like to be able to turn that whole data.frame into a simple feature by making the character column the geometry column of the simple feature with st_as_sfc():

library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, proj.4 4.9.3

geo <- c("{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68152563269095,36.43764870908927]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67408758213843,36.43366018922779]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67708346361097,36.44208638659282]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67886661944996,36.44110273135671]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68089232041565,36.44173155205561]}")

test <- data.frame(numbers = 1:5, letters = letters[1:5], geometry = geo)

test$geometry <- st_as_sfc(test$geometry)
#> OGR: Unsupported geometry type
#> Error in CPL_sfc_from_wkt(x): OGR error

test <- st_as_sf(test)
#> Error in st_sf(x, ..., agr = agr): no simple features geometry column present

Yeah, right now you have to do sth like

> do.call(rbind, lapply(as.character(test$geometry), read_sf))
Simple feature collection with 5 features and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -118.6815 ymin: 36.43366 xmax: -118.6741 ymax: 36.44209
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
                        geometry
1 POINT (-118.681525632691 36...
2 POINT (-118.674087582138 36...
3 POINT (-118.677083463611 36...
4 POINT (-118.67886661945 36....
5 POINT (-118.680892320416 36...

and it would be quite simple to modify the st_as_sfc method for character such that it not only takes WKT, but, with an option set, also GeoJSON.

We now have

> geo <- c("{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68152563269095,36.43764870908927]}",
+          "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67408758213843,36.43366018922779]}",
+          "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67708346361097,36.44208638659282]}",
+          "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67886661944996,36.44110273135671]}",
+          "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68089232041565,36.44173155205561]}")
> st_as_sfc(geo, GeoJSON=TRUE)
Geometry set for 5 features 
geometry type:  POINT
dimension:      XY
bbox:           xmin: -118.6815 ymin: 36.43366 xmax: -118.6741 ymax: 36.44209
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
POINT (-118.681525632691 36.4376487090893)
POINT (-118.674087582138 36.4336601892278)
POINT (-118.677083463611 36.4420863865928)
POINT (-118.67886661945 36.4411027313567)
POINT (-118.680892320416 36.4417315520556)

FWIW, you can also use the wellknown package to achieve this.

library(sf)
library(wellknown)

geo <- c("{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68152563269095,36.43764870908927]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67408758213843,36.43366018922779]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67708346361097,36.44208638659282]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.67886661944996,36.44110273135671]}",
         "{\"geodesic\":true,\"type\":\"Point\",\"coordinates\":[-118.68089232041565,36.44173155205561]}")

wkt = sapply(geo, wellknown::geojson2wkt, USE.NAMES = FALSE)

test <- data.frame(numbers = 1:5, letters = letters[1:5], geometry = wkt)
st_as_sf(test, wkt = "geometry")

Simple feature collection with 5 features and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -118.6815 ymin: 36.43366 xmax: -118.6741 ymax: 36.44209
epsg (SRID):    NA
proj4string:    NA
  numbers letters                       geometry
1       1       a POINT (-118.681525632691 36...
2       2       b POINT (-118.674087582138 36...
3       3       c POINT (-118.677083463611 36...
4       4       d POINT (-118.67886661945 36....
5       5       e POINT (-118.680892320416 36...

Though it doesn't infer crs as @edzer 's new implementation.

"infer" is too much honour: the json doesn't contain it, but the GeoJSON rfc says it only can be 4326, GDAL adheres to that I suppose.

Thank you both! These are great and straightforward solutions.

Hi @edzer , @timelyportfolio, and @tim-salabim.

I'm aware this issue has long been closed. I was wondering if I could get some assistance in converting this data to a spatial object.

https://gist.github.com/denironyx/9fa87a93c10ffa92608362ae770bd7bc

@denironyx looks like you are working on a very neat project. Here is one solution, but I am sure there are others.

 result$data$placesByName %>%
  # convert coordinates to a list of points
  mutate(geom.coordinates = lapply(geom.coordinates, function(x) list(list(matrix(x,ncol=2,byrow=FALSE))))) %>%
  mutate(geometry = mapply(
    function(type, points) {
      # don't know if this is a valid assumption
      #  but get st_* function based on type column
      get(paste0("st_",tolower(type)))(points)
    },
    geom.type,
    geom.coordinates,
    SIMPLIFY = FALSE
  )) %>%
  select(-contains("geom")) %>%
  mutate(geometry = geo$geometry) %>%
  st_sf(crs = "WGS84")

Thanks @timelyportfolio for the assistance. It worked magic!

Was this page helpful?
0 / 5 - 0 ratings