1 Intro

In times of disaster, such as the devastating flood event that occurred in the Ahr valley in Germany in July 2021, effective and efficient routing becomes crucial for emergency responders, humanitarian organizations, and disaster management teams. Openrouteservice, an open-source routing engine built on OpenStreetMap data, offers a powerful and flexible solution for navigation and logistics planning in such scenarios.

This training aims to familiarize participants with the capabilities of openrouteservice within R and demonstrate how it can be leveraged to address routing challenges during a disaster. Through a series of code examples we will explore the key features and functionalities of openrouteservice. You will learn how easy it is to interact with openrouteservice from R and be able to integrate it directly into your own workflow.

2 Setup

Before we delve into the details of utilizing openrouteservice, let’s ensure that we have all the necessary prerequisites in place. The following code chunk provides instructions on how to install the required packages. The openrouteservice R library is currently not published on CRAN, but can be conveniently installed through the devtools packages install_github function:

install.packages(c("devtools", "dplyr","sf","tidyr","geojsonsf","jsonlite",
                   "rjson","units","classInt","tmap","ggplot2",,"mapview","mapedit"))
devtools::install_github("GIScience/openrouteservice-r")

Additional libraries we will need for this training:

# main libraries
library(dplyr) # manipulating tabular data
library(sf) # geospatial classes and datastructures - the swiss army knife on GIS in R
library(openrouteservice) # interface to openrouteservice API

# utils
library(tidyr) # data wrangling
library(geojsonsf) # sf classes to/from geojsons
library(jsonlite) # handling json in R
library(rjson) # another json package for R
library(units) # add units to numerics
library(classInt) # to classify ranges of numerics

# visualization
library(mapview) # simple, quick interactive map visualization
library(mapedit) # create geometries from interactive maps
library(tmap) # visualize spatial data
library(ggplot2) # visualize (non-) spatial data

The main libraries at one glance:

dplyr for the ease of manipulating data.frame classes Wickham H, François R, Henry L, Müller K, Vaughan D (2023). dplyr: A Grammar of Data Manipulation. R package version 1.1.2, https://CRAN.R-project.org/package=dplyr.
sf the dplyr equivalent to work with spatial data.frames. Pebesma, E., & Bivand, R. (2023). Spatial Data Science: With Applications in R (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016
openrouteservice is the package that provides a user firendly R interface to the openrouteservice API. The package is not published on CRAN, therefore we install it via the devtools package from github. https://github.com/GIScience/openrouteservice-r

In order to use the openrouteservice public API you need to register yourself an account and create an API key.

Sign up here: openrouteservice.org/dev/#/login
Or consult our step by step guide as part of the resources we’ve provided you in the workshop book here: giscience.github.io/openrouteservice-workshop

When you sucessfully created your personal API key, copy it over in the chunk below.

your_api_key <-
  "YOUR_API_KEY"

3 Context

During the Ahr valley flood event, many roads and bridges were damaged or rendered impassable, posing significant challenges for navigation and transportation. In such situations, openrouteservice can assist in finding optimal paths, calculating travel distances and times, and identifying alternative routes to reach affected areas or deliver essential supplies.

4 Simple routing and avoid areas

Let’s start to use the openrouteservice library by performing a simple A-to-B routing task. We start by defining an origin and destination and then call the API to create a route out of it. The coordinates provided in the code chunk point to the Nürburgring, a motorsports complex in rural Rhineland-Palatinate, which served as staging area during the emergency response. In this area, rescue helpers from all over Germany were stationed to sleep, eat and repair their equipment inbetween deployments.

4.1 Regular route from Nürburgring to Bad Neuanahr-Ahrweiler

origin <- c(6.943241, 50.334265)
origin_sf <-
  st_as_sf(
    data.frame(longitude = origin[1], latitude = origin[2]),
    coords = c("longitude", "latitude"),
    crs = 4326
  )
destination <- c(7.119166, 50.548979)
destination_sf <- st_as_sf(
  data.frame(longitude = destination[1],
             latitude = destination[2]),
  coords = c("longitude", "latitude"),
  crs = 4326
)

route <- ors_directions(
  list(origin, destination),
  api_key = your_api_key,
  output = "sf",
  instructions = F
)

The ors_directions function consumes a list of numeric vectors which contain the coordinates of origin and destination. You can also add up to 50 waypoints in between. Via the output parameter you can control the structure of the response. The library supports common spatial classes like sf and sp but text and geojson as well. We set instructions=FALSE to prevent the response to contain text based navigation instructions (these can get quite large). Important parameters we did not set, but used the defaults are:

profile: The default is driving-car, however pedestrian, different cycle, heavy vehicles, wheelchair are possible as well
preference: The default is fastest. shortest and recommended are other possible values.
options: Not set by default. It represents an additional object with further parameters to individualize the request like:
- avoid_borders: Do not route through controlled borders, or none at all
- avoid_countries: Do not route through specified countries.
- avoid_features: Do not route via ferries or tolways
- avoid_polygons: Do not route through a customized area. Must be formatted in GeoJSON as either a Polygon or Multipolygon object. Featurecollections are not supported

Check out all available parameters with the interactive API playground here:

openrouteservice.org/dev/#/api-docs/v2/directions/{profile}/post

For the purpose of disaster aware routing we will now use the lastly mentioned parameter and define an area to be avoided by the route calculation.

4.2 Disaster aware route from Nürburgring to Bad Neuanahr-Ahrweiler

First we need to define an area to be avoided from the route calculation. For this purpose we use the mapedit package, which allows us to draw a polygon in an interactive mapview.

Run the following chunk. In the mapview, first zoom out a bit to see the city of Bad Neunahr-Ahrweiler. Next select the polygon or bbox symbol from the left navigation and draw an area around the Ahr river south of Bad Neunahr-Ahrweiler. When you are happy with your drawing, click on the Done button in the lower right corner.

lf <- mapview(destination_sf)

# draw some polygons that we will select later
drawing <- lf |>
  editMap()

## Warning: `select_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `select()` instead.
## ℹ The deprecated feature was likely used in the
##   mapedit package.
##   Please report the issue at
##   <https://github.com/r-spatial/mapedit/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see
## where this warning was generated.

Great! Now we use the created sf representation of your drawn shape and process it into the geojson format, ready to be consumed by the openrouteservice API in the next chunk

options <- list(
  avoid_polygons =
    drawing$finished |>
    st_union() |>
    st_as_sf() |>
    sf_geojson() |>
    rjson::fromJSON()
)

Now we can execute the ors_directions function again, but with the set avoid areas in the options parameter.

disaster_route <- ors_directions(
  list(origin, destination),
  api_key = your_api_key,
  output = "sf",
  instructions = F,
  options = options # here we pass the set options object with the avoid areas to the options parameter
  ,
  profile = "driving-car"
)

Perfect. Let’s visualize the results.

# add a attribute to both routes to be able to differentiate.
disaster_route$type <- "avoid"
route$type <- "regular"

# merge them into one single sf dataframe
routes <- rbind(disaster_route, route)

# we use tmap for visualization in the interactive "view" mode
tmap_mode("view")

## tmap mode set to interactive viewing

tm_basemap("OpenStreetMap") +
  tm_shape(rbind(destination_sf, origin_sf)) +
  tm_symbols() +
  tm_shape(routes) +
  tm_lines("type",
           palette = c("firebrick3", "darkgrey"),
           lwd = 4)

We see the regular route passes much further along to the left. Whereas the disaster aware route swings far to the right, crosses the river east of Bad Neuenahr-Ahrweiler and enters the city from north.

Feel free to play around with the parameters in the chunks provided above. You can change the coordinates and take a closer look at the accompanying data.

For example, what if the staging area had been set up in the area of the city of Koblenz (east) or Bonn (north) instead of at the Nürburgring? What would have been the difference in time and distance?

5 Multiple Routes

Building upon the concepts and skills acquired, we tackle a more complex task. We utilize data from the Copernicus Emergency Mapping service, which provides information on destroyed or damaged roads and bridges in the Ahr valley, to calculate routes from the disaster relief staging area at the Nürburgring. Our objective is to reach 36 locations representing the affected communities in the valley.

For this purpose we need the following additional data:

affected_roads_bridges All damaged and destroyed roads and bridges provided by the Copernicus Emergency Mapping Service. This dataset as well as others are available here: https://emergency.copernicus.eu/mapping/list-of-components/EMSR517
affected_places All place=* objects from OpenStreetMap that are located within the boundaries of municipalities that where affected.

You find the files also in the repo this R-Markdown file is located in. https://github.com/GIScience/openrouteservice-workshop-r/tree/master/data

affected_roads_bridges <- st_read("data/affected_roads.gpkg")

## Reading layer `affected_roads' from data source
##   `/home/mreinmuth/giscience/FOSS4G/openrouteservice-workshop-r/data/affected_roads.gpkg'
##   using driver `GPKG'
## Simple feature collection with 404 features and 0 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 6.7779 ymin: 50.37014 xmax: 7.249583 ymax: 50.55296
## Geodetic CRS:  WGS 84

affected_places <- st_read("data/affected_places.gpkg")

## Reading layer `affected_places' from data source
##   `/home/mreinmuth/giscience/FOSS4G/openrouteservice-workshop-r/data/affected_places.gpkg'
##   using driver `GPKG'
## Simple feature collection with 35 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 6.777893 ymin: 50.37296 xmax: 7.194586 ymax: 50.5515
## Geodetic CRS:  WGS 84

As before we need to prepare the avoid areas for the openrouteservice query. The roads and bridges are of type LINESTRING. The APIs avoid area parameter however can only consume single polygons or multipolygons. Therefore we need convert it via buffer and union into one continuous multipolygon.

affected_roads_bridges <- affected_roads_bridges |>
  st_transform(25832) |>
  st_buffer(2) |>
  st_transform(4326) |>
  st_union() |>
  st_as_sf() |>
  st_make_valid()


options <- list(avoid_polygons =
                  affected_roads_bridges |> sf_geojson() |> rjson::fromJSON())

Next we loop over all communities in affected_places. Some communities won’t be reachable at all, they were cut off by the disaster. Openrouteservice will throw errors for these destinations, as it won’t be able to create a route. We account for this circumstance with a trycatch block.

affected_places$id <- NA
for (a in 1:nrow(affected_places)) {
  tryCatch({
    # request a regular route from staging to community
    directions <- ors_directions(
      list(origin,
           affected_places[a, ] |> st_coordinates() |> c()),
      api_key = your_api_key,
      output = "sf",
      instructions = F
      ,
      profile = "driving-car"
    )
    # request a disaster aware route from staging to community
    directions_avoid <- ors_directions(
      list(origin,
           affected_places[a, ] |> st_coordinates() |> c()),
      api_key = your_api_key,
      output = "sf",
      instructions = F,
      options = options
      ,
      profile = "driving-car"
    )

    # add attribute to differentiate the results later
    directions$type <- "regular"
    directions_avoid$type <- "avoid"
    directions <- rbind(directions, directions_avoid)

    # add id to relate place and route later on
    directions$id <- a
    affected_places[a, ]$id <- a
    # also add the name of the respective destination of the route
    directions$destination <- affected_places[a, ]$name

    # check for the first iteration to init the result sf data frame
    if (a == 1) {
      routes <- directions
      # if not just append to the result sf data frame
    } else {
      routes <- rbind(routes, directions)
    }
  },
  error = function(e)
  {
    # print info on affected communities that openrouteservice could not route to
    print(
      glue::glue(
        "Error on place: {affected_places[a,]$name}; iteration: {a}/{nrow(affected_places)}"
      )
    )
  })
}

## Error on place: Brück 2; iteration: 24/35
## Error on place: Altenburg; iteration: 25/35
## Error on place: Pützfeld; iteration: 32/35
## Error on place: Rech; iteration: 33/35
## Error on place: Schuld 4; iteration: 35/35

We see, the communities of Altenburg, Brück 2, Schuld4, Pützfeld and Rech are unfortunately were no longer reachable via regular motorized transport. These communities were the first to be set up with temporary bridges.

For the communities, however where could generate disaster aware routes, we are now able to determine access paths. Further we can determine the increase in travel time to the communities. For this purpose we need to do some postprocessing of our results. The chunk below will unpack served attributes by openrouteservice on time and distance information for every route. Then a time difference is calculated as new column.

# unpack distance and duration form a list column
routes_tbl <- routes |> unnest_wider(summary)

# add and convert unit information to the column
routes$distance <-
  routes_tbl$distance |> set_units("m") |> set_units("km")
routes$duration <-
  routes_tbl$duration |> set_units("s") |>  set_units("min")

# keep relevant columns only
routes <-
  routes |> select(c(id, destination, type, distance, duration, geometry))

# calculate a new column dur_diff_abs with the absolute difference in time between regular and disaster aware route
routes <- routes |>  mutate(dur_diff_abs =
                              case_when(
                                type == "regular" ~ abs(duration - lead(duration)),
                                type == "avoid" ~ abs(duration - lag(duration)),
                                TRUE ~ NA
                              ))

6 Analysis & Evaluation

In this code chunk, we will take the first step towards analyzing and evaluating the travel time differences between regular routes and routes that bypass flood-affected infrastructure. We will create classes that allow us to extract more meaning and insights from these differences. For this purpose we use the classInt package.

# filter for disaster routes only
routes_avoid <- routes |> filter(type == "avoid")

# create a factor with fixed hierarchy
routes_avoid$destination <- factor(routes_avoid$destination,
                                   levels = routes_avoid$destination)

# create meaningful classes from the numeric differences
intervals <- classIntervals(
  routes_avoid$dur_diff_abs
  |> drop_units(),
  #convert back to numerics
  n = 5,
  # amount of classes
  style = "fixed",
  fixedBreaks = c(0, 1, 5, 10, 20, 50),
  intervalClosure = "right"
)

# based on the classification breaks add a new attribute
routes_avoid$class <- cut(
  routes_avoid$dur_diff_abs |> drop_units(),
  intervals$brks,
  labels = c("A", "B", "C", "D", "E"),
  include.lowest = T
)

Now we will add the classes for route comparison and incorporate them into the affected communities. Additionally, we will create a regular dataframe from the spatial dataframe to create a more lightweight object for non-spatial evaluation in the next step.

# convert from sf dataframe to non spatial dataframe
st_geometry(routes_avoid) <- NULL

# join class attribute to the routes sf dataframe
routes <- routes |>
  left_join(routes_avoid |> select(c("id", "dur_diff_abs", "class")),
            by = "id")

# join class attribute to the affected places sf dataframe
affected_places <- affected_places |>
  left_join(routes_avoid |> select(c("id", "dur_diff_abs", "class")),
            by = "id")

In the following code chunk, we will create a bar plot for each affected community, showcasing the travel time differences between regular routes and disaster-aware routes. Before plotting, we will sort the data based on a leveled factor x variable to ensure a visually appealing and informative representation.

routes_avoid <- routes_avoid |> arrange(desc(dur_diff_abs))
routes_avoid$destination <-
  factor(routes_avoid$destination, levels = routes_avoid$destination)

routes_avoid |>
  ggplot(aes(
    y = dur_diff_abs,
    x = destination,
    fill = class
  ))  +
  geom_bar(stat = "Identity") + # chart type
  labs(x = "Affected Places",
       y = "duration difference after flood") +
  scale_color_grey() +
  theme_classic() +
  theme(axis.text.x = element_text(
    angle = 45,
    vjust = 1,
    hjust = 1
  ))

The bar plot will visually represent these travel time differences, allowing for easy comparison among the affected communities. By sorting the data based on the leveled factor in the x variable, the plot will display the communities in descending traveltime increase.

Class A consists of 6 communities where there is no travel time difference between the regular and disaster-aware routes
Class shows the highest increase in travel time and bears 9 communities. Although the range is more stretched with ~20 min compared to the oter classes
Overall the mean traveltime difference is at 13.64 min and the median at 8.64 min. Not including the 5 not reachable communities.

The last chunk visualizes the routes and affected communities by class

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(affected_places) +
  tm_symbols(size = 1, col = "class") +
  tm_shape(origin_sf) +
  tm_dots("yellow3", size = 1.25, ) +
  tm_shape(routes) +
  tm_lines(col = "lightblue", lwd = 2) +
  tm_facets("type", sync = T)

Do you see any spatial pattern?

Where are class A communities located, where class E?
The communities represented as missing are the non reachable communities. Do these have something in common?

Disaster aware routing with openrouteservice

Marcel Reinmuth

2023-06-26