geoops

Project Status: Active - The project has reached a stable, usable state and is being actively developed. Build Status codecov cran checks rstudio mirror downloads cran version

geoops does spatial operations on GeoJSON.

geoops is inspired by the JS library turf. It’s tagline is Advanced geospatial analysis for browsers and node. Turf works only with GeoJSON, as does geoops. I don’t know JS that well, but it’s easy enough to understand the language, so I’ve been porting Turf to C++ wrapped up in R. The C++ so we can have fast performance. We’ve also wrapped the Turf JS library itself in the package lawn, but we should be able to get better performance out of C++.

geoops has a ways to go to include all the methods that Turf has, but we’ll get there eventually.

This package is alpha stage, expect bugs and changes.

All data is expected to be in WGS-84.

We use a library from Niels Lohmann for working with JSON in C++.

See also:

Package API:

#>  - geo_bearing
#>  - geo_midpoint
#>  - geo_bbox_polygon
#>  - geo_pointgrid
#>  - geo_area
#>  - geo_get_coords
#>  - version
#>  - geo_nearest
#>  - geo_along
#>  - geo_distance
#>  - geo_destination
#>  - geo_trianglegrid
#>  - geo_planepoint
#>  - geo_line_distance

Installation

Stable version

install.packages("geoops")

Dev version

devtools::install_github("ropensci/geoops")
library("geoops")

get json c++ library version info

geoops::version()
#> [1] "{\"compiler\":{\"c++\":\"201103\",\"family\":\"clang\",\"version\":\"10.0.0 (clang-1000.10.25.5)\"},\"copyright\":\"(C) 2013-2017 Niels Lohmann\",\"name\":\"JSON for Modern C++\",\"platform\":\"apple\",\"url\":\"https://github.com/nlohmann/json\",\"version\":{\"major\":3,\"minor\":1,\"patch\":2,\"string\":\"3.1.2\"}}"

distance

Calculate distance between two GeoJSON points

pt1 <- '{
  "type": "Feature",
  "properties": {
    "marker-color": "#f00"
   },
   "geometry": {
      "type": "Point",
      "coordinates": [-75.343, 39.984]
   }
}'
#'
pt2 <- '{
  "type": "Feature",
  "properties": {
     "marker-color": "#0f0"
   },
   "geometry": {
      "type": "Point",
      "coordinates": [-75.534, 39.123]
    }
}'
geo_distance(pt1, pt2)
#> [1] 97.15958

bearing

Calculate bearing between two points

geo_bearing(pt1, pt2)
#> [1] -170.233

destination

geo_destination(pt1, 50, 90, 'miles')
#> [1] "{\"geometry\":{\"coordinates\":[-74.398884,39.98017],\"type\":\"Point\"},\"properties\":{},\"type\":\"Feature\"}"

line distance

Calculate length of GeoJSON LineString or Polygon

line <- '{
  "type": "Feature",
  "properties": {},
  "geometry": {
    "type": "LineString",
    "coordinates": [
      [-77.031669, 38.878605],
      [-77.029609, 38.881946],
      [-77.020339, 38.884084],
      [-77.025661, 38.885821],
      [-77.021884, 38.889563],
      [-77.019824, 38.892368]
    ]
  }
}'
geo_line_distance(line, units = "kilometers")
#> [1] 2.637684

nearest

Calculate nearest point to a reference point

point1 <- '{
  "type": "Feature",
  "properties": {
     "marker-color": "#0f0"
   },
  "geometry": {
     "type": "Point",
     "coordinates": [28.965797, 41.010086]
  }
}'
#'
points <- '{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Point",
        "coordinates": [28.973865, 41.011122]
      }
    }, {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Point",
        "coordinates": [28.948459, 41.024204]
      }
    }, {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "Point",
        "coordinates": [28.938674, 41.013324]
      }
    }
  ]
}'
geo_nearest(point1, points)
#> [1] "{\"geometry\":{\"coordinates\":[28.973865,41.011122],\"type\":\"Point\"},\"properties\":{},\"type\":\"Feature\"}"

comparison to rgeos

distance

library(rgeos)
rgeospt1 <- rgeos::readWKT("POINT(0.5 0.5)")
rgeospt2 <- rgeos::readWKT("POINT(2 2)")
microbenchmark::microbenchmark(
  rgeos = rgeos::gDistance(rgeospt1, rgeospt2),
  geoops = geoops::geo_distance(pt1, pt2, units = "miles"),
  times = 10000L
)
#> Unit: microseconds
#>    expr    min      lq     mean  median      uq      max neval
#>   rgeos 31.780 34.6845 41.64295 35.9510 37.8825 3582.783 10000
#>  geoops 25.588 26.9945 31.30968 28.9785 30.1410 1816.468 10000

nearest

point1 <- '{"type":["Feature"],"properties":{"marker-color":["#0f0"]},"geometry":{"type":["Point"],"coordinates":[28.9658,41.0101]}}'
point2 <- '{"type":["FeatureCollection"],"features":[{"type":"Feature","properties":{},"geometry":{"type":"Point","coordinates":[28.9739,41.0111]}},{"type":"Feature","properties":{},"geometry":{"type":"Point","coordinates":[28.9485,41.0242]}},{"type":"Feature","properties":{},"geometry":{"type":"Point","coordinates":[28.9387,41.0133]}}]}'
g1 <- readWKT("MULTILINESTRING((34 54, 60 34), (0 10, 50 10, 100 50))")
g2 <- readWKT("POINT(100 30)")
microbenchmark::microbenchmark(
  rgeos = rgeos::gNearestPoints(g1, g2),
  geoops = geoops::geo_nearest(point1, points),
  times = 10000L
)
#> Unit: microseconds
#>    expr     min       lq     mean   median       uq       max neval
#>   rgeos 412.753 433.4635 529.4310 447.8365 511.4355 45261.304 10000
#>  geoops  92.894 101.9460 123.2057 116.4675 125.4725  2189.419 10000

Example use case

expand

Get some GeoJSON data, a FeatureCollection of Polygons

file <- system.file("examples/zillow_or.geojson", package = "geoops")
x <- paste0(readLines(file), collapse = "")

Break each polygon into separate JSON string

library("jqr")
polys <- unclass(jq(x, ".features[]"))

Using geo_area, calculate the area of the polygon

areas <- vapply(polys, geo_area, 1, USE.NAMES = FALSE)

Visualize area of the polygons as a histogram

hist(areas, main = "")
plot of chunk unnamed-chunk-20
plot of chunk unnamed-chunk-20

Visualize some of the polygons, all of them

library(leaflet)
leaflet() %>%
  addProviderTiles(provider = "OpenStreetMap.Mapnik") %>%
  addGeoJSON(geojson = x) %>%
  setView(lng = -123, lat = 45, zoom = 7)
plot of chunk unnamed-chunk-21
plot of chunk unnamed-chunk-21

Just one of them

leaflet() %>%
  addProviderTiles(provider = "OpenStreetMap.Mapnik") %>%
  addGeoJSON(geojson = polys[1]) %>%
  setView(lng = -122.7, lat = 45.48, zoom = 13)
plot of chunk unnamed-chunk-22

Meta

ropensci_footer