Get started

openrouteservice R package provides easy access to the openrouteservice (ORS) API from R. It allows you to painlessly consume the following services:

Disclaimer

By using this package, you agree to the ORS terms and conditions.

Installation

The package is not yet available from CRAN, but you can install the development version directly from GitHub.

# install.packages("remotes")
remotes::install_github("GIScience/openrouteservice-r")

Setting up API key

In order to start using ORS services you first need to set up your personal API key, which you can get for free. Once you are signed up, go to https://openrouteservice.org/dev/#/home -> TOKENS. At the bottom of the page you can request a free token (name can be anything).

This will save the key in the default keyring of your system credential store. Once the key is defined, it persists in the keyring store of the operating system. This means that it survives beyond the termination of the R session, so you don’t need to set it again each time you start a new R session. To retrieve the key just call ors_api_key() without the key argument.

Alternatively, they key can be provided in the environment variable ORS_API_KEY. The value from the environment variable takes precedence over the former approach allowing to bypass the keyring infrastructure.

Directions

ors_directions() interfaces the ORS directions service to compute routes between given coordinates.

library(openrouteservice)

coordinates <- list(c(8.34234, 48.23424), c(8.34423, 48.26424))

x <- ors_directions(coordinates)

Way points can be provided as a list of coordinate pairs c(lon, lat), or a 2-column matrix-like object such as a data frame.

coordinates <- data.frame(lon = c(8.34234, 8.34423), lat = c(48.23424, 48.26424))

The response formatting defaults to geoJSON which allows to easily visualize it with e.g. leaflet.

library(leaflet)

leaflet() %>%
  addTiles() %>%
  addGeoJSON(x, fill=FALSE) %>%
  fitBBox(x$bbox)

Other output formats, such as GPX, can be specified in the argument format. Note that plain JSON response returns the geometry as Google’s encoded polyline,

x <- ors_directions(coordinates, format = "json")

geometry <- x$routes[[1]]$geometry
str(geometry)
##  chr "mtkeHuv|q@~@VLHz@\\PR|@hBt@j@^n@L\\NjALv@Jh@NXW`Ay@jDM^kArAoAfBs@^WFY?{Be@[?WJWRi@t@o@`BQRULoAPWHOL]h@mA`C_@d@["| __truncated__

so an additional postprocessing step might be necessary.

library(googlePolylines)
str(decode(geometry))
## List of 1
##  $ :'data.frame':    182 obs. of  2 variables:
##   ..$ lat: num [1:182] 48.2 48.2 48.2 48.2 48.2 ...
##   ..$ lon: num [1:182] 8.34 8.34 8.34 8.34 8.34 ...

The API offers a wide range of profiles for multiple modes of transport, such as: car, heavy vehicle, different bicycle types, walking, hiking and wheelchair. These can be listed with

##                car                hgv               bike           roadbike       mountainbike 
##      "driving-car"      "driving-hgv"  "cycling-regular"     "cycling-road" "cycling-mountain" 
##             e-bike            walking             hiking         wheelchair 
## "cycling-electric"     "foot-walking"      "foot-hiking"       "wheelchair"

Each of these modes uses a carefully compiled street network to suite the profiles requirements.

x <- ors_directions(coordinates, profile="cycling-mountain")

leaflet() %>%
  addTiles() %>%
  addGeoJSON(x, fill=FALSE) %>%
  fitBBox(x$bbox)

Any optional query parameters can be specified by providing them as additional ... arguments to ors_directions. For example, in order to plot the elevation profile of a route colored by steepness use elevation = TRUE to add height to the coordinates of the points along the route and query for steepness in extra_info.

library("sf")

x <- ors_directions(coordinates, profile = "cycling-mountain", elevation = TRUE,
                    extra_info = "steepness", output = "sf")

height <- st_geometry(x)[[1]][, 3]

Here we use simple features output for the sake of easy postprocessing which includes finding the length of individual route segments and their distance relative to the starting point. These can be computed with st_distance() upon converting the LINESTRING to a list of POINTs,

points <- st_cast(st_geometry(x), "POINT")
n <- length(points)
segments <- cumsum(st_distance(points[-n], points[-1], by_element = TRUE))
## st_as_s2(): dropping Z and/or M coordinate
## st_as_s2(): dropping Z and/or M coordinate

while their steepness can be extracted from the requested metadata.

steepness <- x$extras$steepness$values
steepness <- rep(steepness[,3], steepness[,2]-steepness[,1])
steepness <- factor(steepness, -5:5)

palette = setNames(rev(RColorBrewer::brewer.pal(11, "RdYlBu")), levels(steepness))

For the final plot we use ggplot2 in combinations with ggforce which supports handling of length units associated with the data.

library("ggplot2")
library("ggforce")
library("units")

units(height) <- as_units("m")

df <- data.frame(x1 = c(set_units(0, "m"), segments[-(n-1)]),
                 x2 = segments,
                 y1 = height[-n],
                 y2 = height[-1],
                 steepness)

y_ran = range(height) * c(0.9, 1.1)

n = n-1

df2 = data.frame(x = c(df$x1, df$x2, df$x2, df$x1),
                 y = c(rep(y_ran[1], 2*n), df$y2, df$y1),
                 steepness,
                 id = 1:n)

ggplot() + theme_bw() +
  geom_segment(data = df, aes(x1, y1, xend = x2, yend = y2), size = 1) +
  geom_polygon(data = df2, aes(x, y, group = id), fill = "white") +
  geom_polygon(data = df2, aes(x, y , group = id, fill = steepness)) +
  scale_fill_manual(values = alpha(palette, 0.8), drop = FALSE) +
  scale_x_unit(unit = "km", expand = c(0,0)) +
  scale_y_unit(expand = c(0,0), limits = y_ran) +
  labs(x = "Distance", y = "Height")

Advanced options are natively formatted as JSON objects, but can be passed as their R list representation.

polygon = list(
    type = "Polygon",
    coordinates = list(
      list(
        c(8.330469, 48.261570),
        c(8.339052, 48.261570),
        c(8.339052, 48.258227),
        c(8.330469, 48.258227),
        c(8.330469, 48.261570)
      )
    ),
    properties = ""
  )

options <- list(
  avoid_polygons = polygon
)

x <- ors_directions(coordinates, profile="cycling-mountain", options=options)

leaflet() %>%
  addTiles() %>%
  addGeoJSON(polygon, color="#F00") %>%
  addGeoJSON(x, fill=FALSE) %>%
  fitBBox(x$bbox)

Isochrones

Reachability has become a crucial component for many businesses from all different kinds of domains. ors_isochrones() helps you to determine which areas can be reached from certain location(s) in a given time or travel distance. The reachability areas are returned as contours of polygons. Next to the range provided in seconds or meters you may as well specify the corresponding intervals. The list of optional arguments to ors_isochrones() is similar as to ors_directions().

library(mapview)

# embed data in the output file
mapviewOptions(fgb = FALSE)

coordinates <- data.frame(lon = c(8.34234, 8.34234), lat = c(48.23424, 49.23424))

## 1 hour range split into 20 minute intervals
res <- ors_isochrones(coordinates, range = 3600, interval = 1200, output = "sf")
res
## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 7.704047 ymin: 47.7234 xmax: 9.28589 ymax: 49.94667
## Geodetic CRS:  WGS 84
##   group_index value              center                       geometry
## 1           0  1200 8.344268, 48.233826 POLYGON ((8.191246 48.23723...
## 2           0  2400 8.344268, 48.233826 POLYGON ((8.042676 48.19035...
## 3           0  3600 8.344268, 48.233826 POLYGON ((7.884197 48.34021...
## 4           1  1200 8.343705, 49.234211 POLYGON ((8.139146 49.26354...
## 5           1  2400 8.343705, 49.234211 POLYGON ((7.941564 49.20463...
## 6           1  3600 8.343705, 49.234211 POLYGON ((7.704047 49.44966...
values <- levels(factor(res$value))
ranges <- split(res, values)
ranges <- ranges[rev(values)]

names(ranges) <- sprintf("%s min", as.numeric(names(ranges))/60)

mapview(ranges, alpha.regions = 0.2, homebutton = FALSE, legend = FALSE)

Here we have used sf output for the sake of some further postprocessing and visualization. By grouping the isochrones according to ranges we gain the ability of toggling individual ranges when displayed in mapview. Another option could be to group by locations. The following example illustrates a possible approach to applying a custom color palette to the non-overlapping parts of isochrones.

locations = split(res, res$group_index)

locations <- lapply(locations, function(loc) {
  g <- st_geometry(loc)
  g[-which.min(values)] <- st_sfc(Map(st_difference,
                                      g[match(values[-which.min(values)], loc$value)],
                                      g[match(values[-which.max(values)], loc$value)]))
  st_geometry(loc) <- g
  loc
})

isochrones <- unsplit(locations, res$group_index)

pal <- setNames(heat.colors(length(values)), values)
mapview(isochrones, zcol = "value", col = pal, col.regions = pal,
        alpha.regions = 0.5, homebutton = FALSE)

Matrix

One to many, many to many or many to one: ors_matrix() allows you to obtain aggregated time and distance information between a set of locations (origins and destinations). Unlike ors_directions() it does not return detailed route information. But you may still specify the transportation mode and compute routes which adhere to certain restrictions, such as avoiding specific road types or object characteristics.

coordinates <- list(
  c(9.970093, 48.477473),
  c(9.207916, 49.153868),
  c(37.573242, 55.801281),
  c(115.663757,38.106467)
)

# query for duration and distance in km
res <- ors_matrix(coordinates, metrics = c("duration", "distance"), units = "km")

# duration in hours
(res$durations / 3600) %>% round(1)
##       [,1]  [,2] [,3]  [,4]
## [1,]   0.0   1.6 24.2 109.5
## [2,]   1.6   0.0 24.0 109.3
## [3,]  24.2  24.0  0.0  86.3
## [4,] 109.9 109.7 86.5   0.0
# distance in km
res$distances %>% round
##       [,1]  [,2] [,3]  [,4]
## [1,]     0   154 2411 10076
## [2,]   154     0 2386 10051
## [3,]  2358  2334    0  7689
## [4,] 10046 10021 7708     0

Geocoding

ors_geocode() transforms a description of a location provided in query, such as the place’s name, street address or postal code, into a normalized description of the location with a point geometry. Additionally, it offers reverse geocoding which does exactly the opposite: It returns the next enclosing object which surrounds the coordinates of the given location. To obtain more relevant results you may also set a radius of tolerance around the requested coordinates.

## locations of Heidelberg around the globe
x <- ors_geocode("Heidelberg")

leaflet() %>%
  addTiles() %>%
  addGeoJSON(x) %>%
  fitBBox(x$bbox)
## set the number of results returned
x <- ors_geocode("Heidelberg", size = 1)

## search within a particular country
x <- ors_geocode("Heidelberg", boundary.country = "DE")

## structured geocoding
x <- ors_geocode(list(locality="Heidelberg", county="Heidelberg"))

## reverse geocoding
location <- x$features[[1L]]$geometry$coordinates

y <- ors_geocode(location = location, layers = "locality", size = 1)

POIs

This service allows you to find places of interest around or within given geographic coordinates. You may search for given features around a point, path or even within a polygon specified in geometry. To list all the available POI categories use ors_pois('list').

geometry <- list(
  geojson = list(
    type = "Point",
    coordinates = c(8.8034, 53.0756)
  ),
  buffer = 500
)

ors_pois(
  request = 'pois',
  geometry = geometry,
  limit = 2000,
  sortby = "distance",
  filters = list(
    category_ids = 488,
    wheelchair = "yes"
  ),
  output = "sf"
)
## Simple feature collection with 3 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 8.80294 ymin: 53.07452 xmax: 8.80903 ymax: 53.0757
## Geodetic CRS:  WGS 84
##                                                 category_ids  distance             osm_tags
## 1 {"488":{"category_name":"kiosk","category_group":"shops"}}  32.95127 {"wheelchair":"yes"}
## 2 {"488":{"category_name":"kiosk","category_group":"shops"}} 292.89365 {"wheelchair":"yes"}
## 3 {"488":{"category_name":"kiosk","category_group":"shops"}} 395.87884 {"wheelchair":"yes"}
##   osm_type     osm_id                 geometry
## 1        1  726580652  POINT (8.80294 53.0757)
## 2        1 2525463951 POINT (8.80777 53.07559)
## 3        1 4827632819 POINT (8.80903 53.07452)

You can gather statistics on the amount of certain POIs in an area by using request='stats'.

ors_pois(
  request = 'stats',
  geometry = geometry,
  limit = 2000,
  sortby = "distance",
  filters = list(category_ids = 488)
  )
## <ors_pois>
## List of 1
##  $ :List of 2
##   ..$ places     :List of 2
##   .. ..$ total_count: int 8
##   .. ..$ shops      :List of 3
##   .. .. ..$ group_id   : int 420
##   .. .. ..$ categories :List of 1
##   .. .. .. ..$ kiosk:List of 2
##   .. .. .. .. ..$ count      : int 8
##   .. .. .. .. ..$ category_id: int 488
##   .. .. ..$ total_count: int 8
##   ..$ information:List of 4
##   .. ..$ attribution: chr "openrouteservice.org | OpenStreetMap contributors"
##   .. ..$ version    : chr "0.1"
##   .. ..$ timestamp  : int 1634138573
##   .. ..$ query      :List of 5
##   .. .. ..$ request : chr "stats"
##   .. .. ..$ limit   : int 2000
##   .. .. ..$ sortby  : chr "distance"
##   .. .. ..$ filters :List of 1
##   .. .. .. ..$ category_ids: int 488
##   .. .. ..$ geometry:List of 2
##   .. .. .. ..$ geojson:List of 2
##   .. .. .. .. ..$ type       : chr "Point"
##   .. .. .. .. ..$ coordinates: num [1:2] 8.8 53.1
##   .. .. .. ..$ buffer : int 500

Elevation

Given a point or line geometry you can use ors_elevation to query for its elevation.

x <- ors_geocode("Königstuhl", output = "sf")

ors_elevation("point", st_coordinates(x))
## <ors_elevation>
##  num [1:3] 8.72 49.4 560

Optimization

The optimization endpoint solves the vehicle routing problem (VRP) of finding an optimal set of routes for a fleet of vehicles to traverse in order to deliver to a given set of locations. The service is based on Vroom and can be used to schedule multiple vehicles and jobs respecting time windows, capacities and required skills. VRP generalizes the classic traveling salesman problem of finding the fastest or shortest possible route that visits a given list of locations.

The following example involves a 2-vehicle fleet carrying out deliveries across 6 locations.

home_base <- data.frame(lon = 2.370658, lat = 48.721666)

vehicles = vehicles(
  id = 1:2,
  profile = "driving-car",
  start = home_base,
  end = home_base,
  capacity = 4,
  skills = list(c(1, 14), c(2, 14)),
  time_window = c(28800, 43200)
)

Both vehicles share the start/end points and have the same capacity, but differ in the set of skills assigned. We are interested in using them to serve a number of jobs with certain skills requirements between locations. These skills are mandatory, which means a given job can only be served by a vehicle that has all its required skills.

locations <- list(
  c(1.98806, 48.705),
  c(2.03655, 48.61128),
  c(2.39719, 49.07611),
  c(2.41808, 49.22619),
  c(2.28325, 48.5958),
  c(2.89357, 48.90736)
)

jobs = jobs(
  id = 1:6,
  service = 300,
  amount = 1,
  location = locations,
  skills = list(1, 1, 2, 2, 14, 14)
)

The helper functions vehicles and jobs produce data.frames which have the format appropriate for ors_optimization. Route geometries are enabled by setting the corresponding flag in options.

res <- ors_optimization(jobs, vehicles, options = list(g = TRUE))

The geometries are returned as Google’s encoded polylines, so for visualization in leaflet they need to be decoded. Furthermore, we extract the job locations from the response such that we can label them in the order in which they are visited along the routes.

lapply(res$routes, with, {
  list(
    geometry = googlePolylines::decode(geometry)[[1L]],
    locations = lapply(steps, with, if (type=="job") location) %>%
      do.call(rbind, .) %>% data.frame %>% setNames(c("lon", "lat"))
  )
  }) -> routes

## Helper function to add a list of routes and their ordered waypoints
addRoutes <- function(map, routes, colors) {
  routes <- mapply(c, routes, color = colors, SIMPLIFY = FALSE)
  f <- function (map, route) {
    with(route, {
      labels <- sprintf("<b>%s</b>", 1:nrow(locations))
      markers <- awesomeIcons(markerColor = color, text = labels, fontFamily = "arial")
      map %>%
        addPolylines(data = geometry, lng = ~lon, lat = ~lat, col = ~color) %>%
        addAwesomeMarkers(data = locations, lng = ~lon, lat = ~lat, icon = markers)
    })
  }
  Reduce(f, routes, map)
}

leaflet() %>%
  addTiles() %>%
  addAwesomeMarkers(data = home_base, icon = awesomeIcons("home")) %>%
  addRoutes(routes, c("purple", "green"))