geoops

Build Status codecov 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\":\"9.0.0 (clang-900.0.39.2)\"},\"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\":1,\"string\":\"3.1.1\"}}"

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 25.166 28.2415 35.69903 29.4915 31.7615 18693.60 10000
#>  geoops 23.240 24.9910 33.07648 26.5655 28.7375 41458.95 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 370.119 391.0525 463.1019 406.9535 453.951 14023.597 10000
#>  geoops  85.386  94.9875 113.6604 108.6060 117.964  1083.361 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