Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Convert features of a 'multifeature' GeoJSON into R spatial objects

Normally you can read geojson files into R with trusty readOGR, as illustrated here.

However, this fails for multifeature geojsons.

Reproducible example:

downloader::download("https://github.com/Robinlovelace/Creating-maps-in-R/raw/master/data/test-multifeature.geojson", "test.geojson")
test <- rgdal::readOGR("test.geojson", "OGRGeoJSON") # fails with:

Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv,  : 
  Multiple incompatible geometries: wkbPoint: 98; wkbLineString: 660; wkbPolygon: 23

The error message is clear-enough and indicates a solution: split the features. Aside from doing this with regex, I don't know how, however.

Any ideas very much welcome.

The amazing thing: GitHub displays the data natively on the browser, whereas R cannot even (seemingly) read it in!

Alternative way to a solution:

test <- geojsonio::geojson_read("test.geojson")
like image 958
RobinLovelace Avatar asked Jun 01 '15 20:06

RobinLovelace


3 Answers

You can use the require_geomType parameter for various GDAL functions to extract the features that you need:

library(rgdal)

ogrListLayers("test.geojson")
## [1] "OGRGeoJSON"
## attr(,"driver")
## [1] "GeoJSON"
## attr(,"nlayers")
## [1] 1

# This fails but you can at least see the geoms it whines about
ogrInfo("test.geojson", "OGRGeoJSON")
## Error in ogrInfo("test.geojson", "OGRGeoJSON") : 
##   Multiple incompatible geometries: wkbPoint: 98; wkbLineString: 660; wkbPolygon: 23


ogrInfo("test.geojson", "OGRGeoJSON", require_geomType="wkbPoint")
## NOTE: keeping only 98 wkbPoint of 781 features
##     note that extent applies to all features
## Source: "test.geojson", layer: "OGRGeoJSON"
## Driver: GeoJSON number of rows 781 
##   selected geometry type: wkbPoint with 98 rows
## Feature type: wkbPoint with 2 dimensions
## Extent: (12.48326 41.88355) - (12.5033 41.89629)
## CRS: +proj=longlat +datum=WGS84 +no_defs  
## Number of fields: 78 
##                        name type length typeName
## 1                      area    4      0   String
## 2                   bicycle    4      0   String
## ...
## LONG LIST - 78 total


ogrInfo("test.geojson", "OGRGeoJSON", require_geomType="wkbLineString")
## NOTE: keeping only 660 wkbLineString of 781 features
##     note that extent applies to all features
## Source: "test.geojson", layer: "OGRGeoJSON"
## Driver: GeoJSON number of rows 781 
##   selected geometry type: wkbLineString with 660 rows
## Feature type: wkbLineString with 2 dimensions
## Extent: (12.48326 41.88355) - (12.5033 41.89629)
## CRS: +proj=longlat +datum=WGS84 +no_defs  
## Number of fields: 78 
##                        name type length typeName
## 1                      area    4      0   String
## 2                   bicycle    4      0   String
## ...
## LONG LIST - 78 total (same as above)


ogrInfo("test.geojson", "OGRGeoJSON", require_geomType="wkbPolygon")
## NOTE: keeping only 23 wkbPolygon of 781 features
##     note that extent applies to all features
## Source: "test.geojson", layer: "OGRGeoJSON"
## Driver: GeoJSON number of rows 781 
##   selected geometry type: wkbPolygon with 23 rows
## Feature type: wkbPolygon with 2 dimensions
## Extent: (12.48326 41.88355) - (12.5033 41.89629)
## CRS: +proj=longlat +datum=WGS84 +no_defs  
## Number of fields: 78 
##                        name type length typeName
## 1                      area    4      0   String
## 2                   bicycle    4      0   String
## ...
## LONG LIST - 78 total (same as above)


points <- readOGR("test.geojson", "OGRGeoJSON", require_geomType="wkbPoint")
## OGR data source with driver: GeoJSON 
## Source: "test.geojson", layer: "OGRGeoJSON"
## with 781 features;
## Selected wkbPoint feature type, with 98 rows
## It has 78 fields
## NOTE: keeping only 98 wkbPoint of 781 features

lines <- readOGR("test.geojson", "OGRGeoJSON", require_geomType="wkbLineString")
## OGR data source with driver: GeoJSON 
## Source: "test.geojson", layer: "OGRGeoJSON"
## with 781 features;
## Selected wkbLineString feature type, with 660 rows
## It has 78 fields
## NOTE: keeping only 660 wkbLineString of 781 features

polygons <- readOGR("test.geojson", "OGRGeoJSON", require_geomType="wkbPolygon")
## OGR data source with driver: GeoJSON 
## Source: "test.geojson", layer: "OGRGeoJSON"
## with 781 features;
## Selected wkbPolygon feature type, with 23 rows
## It has 78 fields
## NOTE: keeping only 23 wkbPolygon of 781 features

# prove they red in things
plot(lines, col="#7f7f7f")
plot(polygons, add=TRUE)
plot(points, add=TRUE, col="red")

enter image description here

like image 140
hrbrmstr Avatar answered Nov 19 '22 19:11

hrbrmstr


You can use ogr2ogr on the command line to split this monstrous chimera into sensible things:

ogr2ogr -where "OGR_GEOMETRY='LINESTRING'" \
     -f "GeoJSON" lines.geojson  test.geojson

and similarly for points and polygons.

There was some chatter a few years ago about implementing OGR_SQL into readOGR, at which point you would be able to do this from R, but Roger didn't want to do it and nobody wanted to help :(

Once you've created the split geojson files you can read them into a single rgeos::SpatialCollections object:

points=readOGR("points.geojson","OGRGeoJSON")
polys=readOGR("polygons.geojson","OGRGeoJSON")
lines=readOGR("lines.geojson","OGRGeoJSON")
require(rgeos)
g = SpatialCollections(points=points, lines=lines, polygons=polys)
plot(g)

If you want to try something with the geojsonio then you can use Filter to select list elements of a given geometry from the Geometry Collection

polygon_features = Filter(
    function(f){f$geometry$type=="Polygon"},
    test$features)

but then you still have to build something you can get into separate R entities....

like image 41
Spacedman Avatar answered Nov 19 '22 19:11

Spacedman


A few years later, a couple of alternatives - library(geojsonsf) and library(sf) will both read the geojson and convert to sf objects

url <- 'https://github.com/Robinlovelace/Creating-maps-in-R/raw/master/data/test-multifeature.geojson'

## these give the same result
sf <- geojsonsf::geojson_sf( url )
sf <- sf::st_read( url )

Let's take a look

library(mapdeck)

set_token( "MAPBOX_TOKEN" )

mapdeck( style = mapdeck_style("light") ) %>%
    add_sf( sf )

enter image description here

like image 35
SymbolixAU Avatar answered Nov 19 '22 20:11

SymbolixAU