This content was presented to Nelson\Nygaard Staff at a Lunch and Learn webinar on Thursday, August 11th, 2022, and is available as a recording here and embedded below.
Today’s Agenda
- Introduction (10 min)
- Routing & Origin-Destination Analysis (20 min)
- Map Matching (15 min)
- Isochrones (15 min)
Introduction
Motivation
Welcome back to the R & Python webinar series! I hope to keep these going from here on out – the first quarter of the year was extremely busy and then it is always difficult to get things started again. But hopefully this extended break has given us all some time to think about topics we would like to pursue in future webinars – please post in #r-python-users about topics you would like us to discuss in future sessions!
Today we are going to be talking about Valhalla – no, not the one from Norse mythology – an open source routing engine for the OpenStreetMap street network. Recently, I have set up NN’s very own Valhalla instance (hosted in the cloud on DigitalOcean) to be utilized via an HTTP API - e.g., with R or Python (and hopefully eventually ArcGIS through either R or Python). I originally was motivated to set this up for map matching (which we will talk about), but Valhalla can also be used for a variety of other street network routing related tasks. Today we will talk about a) routing and origin-destination analysis, b) map matching, and c) isochrones. There are a few other tools available through Valhalla we can discuss in a future session if there is interest – the full documentation is here.
One major limitation of Valhalla I would like to solve for in the future by creating our own instance of OpenTripPlanner (OTP)(a multi-modal trip planner) is that Valhalla only works for single modes of travel using the street network, so it does not work for transit routing, which is by its nature always a multi-modal trip. So it cannot be used to generate transit-mode isochrones, for example, because these would require the incorporation of both transit and walk elements. An additional note – as of this writing, Valhalla’s documentation also makes reference to multimodal trips being newly available in its beta version, which I can experiment with as well. So stay tuned for an NN implementation of OpenTrip Planner, or the deployed beta version of Valhalla!
Valhalla can be accessed through the HTTP API directly (though this requires a bit more know-how), or you can use one of the functions included in the {valhallr} package as an API wrapper. As of this writing, the CRAN version of {valhallr} does not have a wrapper function for the map matching feature, so Bryan wrote that new function into [a forked version of the package}(https://github.com/bpb824/valhallr). For now, to access the map matching function, you will have to install the version of {valhallr} from Bryan’s repository – Bryan has filed a pull request to get his map matching wrapper integrated into the overall package, but that will depend on the original package’s author to accept.
Today’s Example Data
Today we are going to use TriMet’s most recent GTFS feed to get some sample data to work with – download the data yourself if you want to test locally!
We are going to filter out TriMet’s transit centers from the stop data, and use those locations both for origin-destination routing and travel time matrix analysis and for pedestrian isochrone generation. Additionally, we’re going to go through how to match a bus route to the OSM street network using TriMet’s route 12, which runs along Barbur and Sandy Boulevards (primarily) between Tigard Transit Center and Parkrose Transit Center. I chose this route because it is a cross-town route that travels through a variety of different regional contexts, and the map matching holds up pretty well.
#To install Bryan's version of `{valhallr}`
#devtools::install_github("bpb824/valhallr")
library(valhallr)
library(tidyverse)
library(sf)
library(tidytransit)
library(leaflet)
library(osmdata)
library(scales)
nn_valhalla_hostname <- "128.199.8.29"
coord_global = 4326
coord_local = 2269 #see here: https://nelsonnygaard.shinyapps.io/coord-system-reference/
trimet_gtfs <- read_gtfs("sample-data/trimet-2022-07-23.zip")
trimet_stops <- trimet_gtfs$stops
#Filtering down to just transit center stops
tc_stops <- trimet_stops %>%
filter(str_detect(tts_stop_name,"transit.center"))
trimet_tc_names = c("Barbur", "Rose Quarter", "Gresham", "Hillsboro",
"Lake Oswego", "N Lombard", "Tigard", "Beaverton",
"Gateway", "Hollywood", "Oregon City", "Willow Creek",
"Sunset", "Washington Square", "Clackamas Town Center",
"Parkrose")
grouped_tc_locs <- tc_stops %>%
mutate(transit_center_name = map_chr(stop_name, function(stop_name){
matched_name <- trimet_tc_names[str_detect(stop_name,trimet_tc_names)]
if(length(matched_name)>0){
return(matched_name)
}else{
return(NA_character_)
}
})) %>%
group_by(transit_center_name) %>%
summarise(stop_lon = mean(stop_lon),
stop_lat = mean(stop_lat)) %>%
st_as_sf(coords = c("stop_lon","stop_lat"), crs=coord_global)
#Filtering down to a route shape to match
top_rt_12_shp_id <- trimet_gtfs$trips %>%
filter(route_id == 12) %>%
group_by(shape_id) %>%
summarise(n=n()) %>%
arrange(desc(n)) %>%
pull(shape_id) %>%
head(1)
rt_shp <- trimet_gtfs$shapes %>%
filter(shape_id == top_rt_12_shp_id) %>%
select(shape_pt_lon,shape_pt_lat) %>%
as.matrix() %>%
st_linestring() %>%
st_sfc(crs = coord_global)
leaflet() %>%
addTiles() %>%
addAwesomeMarkers(data = grouped_tc_locs, label=~transit_center_name,
icon = makeAwesomeIcon(icon = "bus", library = "fa")) %>%
addPolylines(data = rt_shp, label="Route 12 Sample Shape")Routing & Origin-Destination Analysis
To demonstrate the routing functionality of Valhalla, we will use the termini of Route 12 – Parkrose Transit Center and Tigard Transit Center. To keep it simple, I will demonstrate how to use the "auto" routing (input using the costing parameter), but you can also use the "pedestrian" or "bicycle" costing models as well. Note that Valhalla does not have any traffic model associated with the roadways, so travel times are not reflective of traffic conditions at a particular time of day – the speed used is the speed assigned to each link in the underlying OpenStreetMap data. There are also additional options to tweak the route selection, the full API documentation has descriptions of all possible parameters, and then the R documentation has some high level descriptions.
Some key aspects of the below code to note:
- The IP address of NN’s Valhalla routing service is
128.199.8.29– this is stored in a variable at the top of this tutorial and used throughout. It must be used in each Valhalla API wrapper function with thehostnameparameter. Additional security may be added to this service in the future, but for now it is going to be left open for ease of use, unless any issues are encountered. - Valhalla, like Google, compresses polylines into an encoded string for data compression purposes, and so to access any polylines coming out of Valhalla, they muse be decoded using the
decode()function. - Valhalla’s routing function offers a number of parameters for altering how routes are selected. As an example below, I query an auto route avoiding highways, but there are many other parameters than can be set, the full list is in the documentation.
origin_loc <- grouped_tc_locs %>%
filter(transit_center_name == "Parkrose") %>%
st_coordinates() %>%
as_tibble() %>%
#Valhalla requires location data in tibble format, with the columns `lon` and `lat`
rename(lon = X, lat = Y)
dest_loc <- grouped_tc_locs %>%
filter(transit_center_name == "Tigard") %>%
st_coordinates() %>%
as_tibble() %>%
#Valhalla requires location data in tibble format, with the columns `lon` and `lat`
rename(lon = X, lat = Y)
#Regular auto route
auto_od_route_res_obj <- route(from = origin_loc, to = dest_loc,
costing = "auto", hostname = nn_valhalla_hostname)
#Extracting shape geometry
auto_od_route_shp <- auto_od_route_res_obj$legs$shape %>%
decode() %>% #Decode Valhalla polyline string
select(lng,lat) %>%
as.matrix() %>%
st_linestring() %>%
st_sfc(crs = coord_global)
#Auto route avoiding highways
auto_od_route_res_obj_avoid_highways <- route(from = origin_loc, to = dest_loc,
costing = "auto", hostname = nn_valhalla_hostname,
costing_options = list(
use_highways = 0
))
#Extracting shape geometry
auto_od_route_shp_avoid_highways <- auto_od_route_res_obj_avoid_highways$legs$shape %>%
decode() %>% #Decode Valhalla polyline string
select(lng,lat) %>%
as.matrix() %>%
st_linestring() %>%
st_sfc(crs = coord_global)
#Generate map of route shapes and origin and destination
leaflet() %>%
addTiles() %>%
addPolylines(data = rt_shp, color="green", label="Route 12 (WB) Alignment") %>%
addPolylines(data = auto_od_route_shp, color="red", label="Valhalla Auto Preferred Route") %>%
addPolylines(data = auto_od_route_shp_avoid_highways, color="purple", label="Valhalla Auto Route Avoiding Highways") %>%
addAwesomeMarkers(data = origin_loc, icon = makeAwesomeIcon(icon = "play", library = "fa", markerColor = "green")) %>%
addAwesomeMarkers(data = dest_loc, icon = makeAwesomeIcon(icon = "stop", library = "fa", markerColor = "red"))#check out summary object
auto_od_route_res_obj$summary$has_time_restrictions
[1] FALSE
$min_lat
[1] 45.43042
$min_lon
[1] -122.7702
$max_lat
[1] 45.5615
$max_lon
[1] -122.5605
$time
[1] 1343.255
$length
[1] 28.512
$cost
[1] 1851.151#Check out maneuvers object
auto_od_route_res_obj$legs$maneuvers[[1]] type
1 2
2 20
3 23
4 9
5 19
6 20
7 24
8 24
9 20
10 15
11 15
12 5
instruction
1 Drive northeast on I 205 North.
2 Take exit 23B on the right onto US 30 Bypass West toward Killingsworth Street.
3 Keep right to take US 30 Bypass West toward Killingsworth Street.
4 Bear right onto Northeast Lombard Street/US 30 Bypass.
5 Turn left to take the I 205 South ramp toward I 84/Portland/Salem/Oregon City.
6 Take exit 21B on the right onto I 84 West toward Portland.
7 Keep left to take I 5 South toward Salem.
8 Keep left to take I 5 South toward Salem.
9 Take exit 294 on the right onto OR 99W South toward Tigard/Newberg/McMinnville/Lincoln City.
10 Turn left onto Southwest Main Street.
11 Turn left onto Southwest Commercial Street.
12 Your destination is on the right.
verbal_succinct_transition_instruction
1 Drive northeast. Then Take exit 23B on the right.
2 <NA>
3 <NA>
4 Bear right. Then Turn left to take the I 205 South ramp.
5 <NA>
6 <NA>
7 <NA>
8 <NA>
9 <NA>
10 Turn left.
11 Turn left.
12 <NA>
verbal_pre_transition_instruction
1 Drive northeast on I 205 North. Then Take exit 23B on the right.
2 Take exit 23B on the right onto US 30 Bypass West toward Killingsworth Street.
3 Keep right to take US 30 Bypass West toward Killingsworth Street.
4 Bear right onto Northeast Lombard Street, US 30 Bypass. Then Turn left to take the I 205 South ramp.
5 Turn left to take the I 205 South ramp toward I 84, Portland.
6 Take exit 21B on the right onto I 84 West toward Portland.
7 Keep left to take I 5 South toward Salem.
8 Keep left to take I 5 South toward Salem.
9 Take exit 294 on the right onto OR 99W South toward Tigard, Newberg.
10 Turn left onto Southwest Main Street.
11 Turn left onto Southwest Commercial Street.
12 Your destination is on the right.
verbal_post_transition_instruction
1 Continue for 30 meters.
2 <NA>
3 <NA>
4 Continue for 200 meters.
5 Continue for 3 kilometers.
6 Continue for 9 kilometers.
7 <NA>
8 Continue for 12 kilometers.
9 Continue for 3 kilometers.
10 Continue for 400 meters.
11 Continue for 100 meters.
12 <NA>
street_names
1 I 205 North
2 NULL
3 NULL
4 Northeast Lombard Street, US 30 Bypass
5 I 205 South
6 I 84 West, US 30, Banfield Freeway, T. H. Banfield Expressway
7 NULL
8 I 5 South
9 Southwest Pacific Highway, OR 99W
10 Southwest Main Street
11 Southwest Commercial Street
12 NULL
time length cost begin_shape_index end_shape_index
1 1.339 0.033 1.540 0 2
2 40.106 0.444 46.057 2 21
3 18.675 0.236 20.976 21 34
4 8.298 0.165 40.737 34 40
5 135.313 3.199 207.577 40 137
6 359.951 8.692 476.436 137 375
7 23.833 0.529 29.161 375 407
8 511.300 11.978 609.205 407 681
9 178.151 2.758 276.853 681 786
10 49.015 0.350 97.493 786 819
11 17.270 0.124 45.112 819 828
12 0.000 0.000 0.000 828 828
verbal_multi_cue travel_mode travel_type
1 TRUE drive car
2 NA drive car
3 NA drive car
4 TRUE drive car
5 NA drive car
6 NA drive car
7 NA drive car
8 NA drive car
9 NA drive car
10 NA drive car
11 NA drive car
12 NA drive car
verbal_transition_alert_instruction
1 <NA>
2 Take exit 23B on the right.
3 Keep right to take US 30 Bypass West.
4 Bear right onto Northeast Lombard Street.
5 Turn left to take the I 205 South ramp.
6 Take exit 21B on the right.
7 Keep left to take I 5 South.
8 Keep left to take I 5 South.
9 Take exit 294 on the right.
10 Turn left onto Southwest Main Street.
11 Turn left onto Southwest Commercial Street.
12 Your destination will be on the right.
sign.exit_number_elements sign.exit_branch_elements
1 NULL NULL
2 23B US 30 Bypass West, 1
3 NULL US 30 Bypass West, 1
4 NULL NULL
5 NULL I 205 South, 1
6 21B I 84 West, US 30 West, 1, NA
7 NULL I 5 South, 2
8 NULL I 5 South, 2
9 294 OR 99W South
10 NULL NULL
11 NULL NULL
12 NULL NULL
sign.exit_toward_elements
1 NULL
2 Killingsworth Street, 1
3 Killingsworth Street, 1
4 NULL
5 I 84, Portland, Salem, Oregon City
6 Portland
7 Salem, OMSI, CIty Center, 1, NA, NA
8 Salem, Beaverton, 1, NA
9 Tigard, Newberg, McMinnville, Lincoln City
10 NULL
11 NULL
12 NULLSometimes what you may really be interested in routing is actually understanding travel distances or times between a variety of origins and destinations – Valhalla offers an origin-destination matrix function build expressly for this that is quicker than querying each route individually – this is accessible using the sources_to_targets() function in R.
#Creating a tibble with latitude and longitude as columns instead of using an sf object (per Valhalla needs)
grouped_tc_coords <- bind_cols(
grouped_tc_locs %>% st_drop_geometry(),
grouped_tc_locs %>% st_coordinates() %>%
as_tibble() %>%
rename(lon = X, lat = Y)
) %>%
#Valhalla spits out resulting OD travel times and distances by an index starting at 0
mutate(index = seq_along(transit_center_name) -1)
#Calling OD matrix API wrapper
od_res <- sources_to_targets(froms = grouped_tc_coords,
tos = grouped_tc_coords,
costing = "auto",
hostname = nn_valhalla_hostname)
#Joining in names of transit centers based on index, converting travel time and distance to more convenient units
od_res_joined <- od_res %>%
left_join(grouped_tc_coords %>%
select(index,transit_center_name) %>%
rename(from_index = index,
from_transit_center = transit_center_name)) %>%
left_join(grouped_tc_coords %>%
select(index,transit_center_name) %>%
rename(to_index = index,
to_transit_center = transit_center_name)) %>%
mutate(tt_minutes = time/60, #minutes
dist_miles = distance*0.621371) # miles
#Creating a graphical matrix of the travel times, with each cell colored by travel time in minutes
ggplot(od_res_joined,aes(x = from_transit_center, y = to_transit_center, fill = tt_minutes))+
geom_tile()+
geom_text(aes(label = comma(tt_minutes,accuracy=0.1)),
size=2.5,color="white")+
scale_fill_viridis_c(option = "A",direction = -1)+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
labs(x = "Origin Transit Center", y = "Destination Transit Center",
fill="Auto\nTravel Time\n(minutes)",
title = "Auto Travel Times between TriMet Transit Centers")
Map Matching
An issue that comes up across many NN projects, but especially in dealing with transit route shapes, is matching a linear geometry to a street network, so that further calculations can be done for specific street network links. In this case, we will be matching a transit route shape from GTFS (extracted above) to the OSM street network using Valhalla.
Map matching requires an input geometry (a tibble with lon and lat values) in the shape parameter and a costing model for matching – in this case we will use costing = "bus". The output will be a list with a variety of valuable information, including statistics on the performance of the overall match and matching of particular points in the trajectory, and details about every link. However, we will extract just two pieces – a) the matched “trace” (a polyline) representing the path the matching algorithm thinks was taken along the street network and b) the series of OSM way IDs that are used by the matched trace.
Note that in OSM parlance, street network links are known as “ways” and in network mathematics parlance these are known as “edges” – I will use the terms interchangeably below.
#Preparing input geometry
rt_shp_coords <- rt_shp %>%
st_coordinates() %>%
as_tibble() %>%
rename(lon = X, lat = Y)
#Map matching API call
match_res_obj <- trace_attributes(shape = rt_shp_coords, costing = "bus", hostname = nn_valhalla_hostname)
#Transforming matched trace into an sf linestring
matched_rt_shape <- match_res_obj$shape %>%
decode() %>% #Decode Valhalla polyline string
select(lng,lat) %>%
as.matrix() %>%
st_linestring() %>%
st_sfc(crs = coord_global)
#This will return a record of the current edge each time the trace passes through an OSM node
#There will be some repeat edges as OSM ways do not neccessarily break at each OSM node.
edge_tibble <- tibble(
way_id = as.character(match_res_obj$edges$way_id)
)
#Below I am creating an edge sequence --
#I do it using a loop just in case an edge happens to be used twice, which is the case on this route at Barbur transit center
extracted_edges_pre <- edge_tibble %>%
select(way_id) %>%
mutate(edge_sequence =NA)
extracted_edges_pre$edge_sequence[1] =1
for(j in 2:nrow(extracted_edges_pre)){
if(extracted_edges_pre$way_id[j] == extracted_edges_pre$way_id[j-1]){
extracted_edges_pre$edge_sequence[j] = extracted_edges_pre$edge_sequence[j-1]
}else{
extracted_edges_pre$edge_sequence[j] = extracted_edges_pre$edge_sequence[j-1] + 1
}
}
#Finally I can extract the ordered sequence of edges used
extracted_edges <- extracted_edges_pre %>%
distinct(way_id,edge_sequence)
#Now I can fetch the relevant OSM geometry for each edge using the way ID
osm_query_res = osmdata::opq_osm_id(id = str_flatten(extracted_edges$way_id, collapse = ", "),
type = "way") %>%
osmdata_sf()
#Create an OSM link geometry for mapping
osm_link_geoms <-osm_query_res$osm_lines %>%
select(osm_id, name, highway, geometry) %>%
rename(osm_link_id = osm_id) %>%
left_join(extracted_edges %>%
rename(osm_link_id = way_id)) %>%
arrange(edge_sequence) %>%
st_as_sf()
#Map OSM link geometry, Valhalla trace, and original route shape to compare.
#Match looks excellent!
leaflet() %>%
addTiles() %>%
addPolylines(data = rt_shp, color="green", label="Route 12 (WB) Alignment (from GTFS)", group = "GTFS Route Shape") %>%
addPolylines(data = matched_rt_shape, color="black", label="Valhalla Matched Bus Route", group = "Valhalla Matched Trace") %>%
addPolylines(data = osm_link_geoms, label=~ paste0(edge_sequence,": ",osm_link_id), group = "Matched OSM Links",
color="red") %>%
addLayersControl(overlayGroups = c("GTFS Route Shape","Valhalla Matched Trace", "Matched OSM Links"))Isochrones
Finally, Valhalla can be used to generate isochrones for pedestrian, bicycle and auto costing models – stay tuned for use with a multimodal costing model in the future. Isochrones can be generated given an origin location, a costing model, a metric to use (typically time in minutes, but you can also use distance), and then breakpoints to use for each contour. We will generate some walking isochrones from one transit center, and then from all of them at once.
single_iso <- grouped_tc_locs %>%
filter(transit_center_name == "Parkrose") %>%
st_coordinates() %>%
as_tibble() %>%
rename(lon = X, lat = Y) %>%
#Piping origin location directly into isochrone function
isochrone(costing = "pedestrian",
contours = c(5,10,15,20),
hostname = nn_valhalla_hostname)
#Generating map
leaflet() %>%
addTiles() %>%
addPolygons(data =single_iso, fillColor =~ fillColor,
color="white", opacity = 0.9, weight = 2,
fillOpacity = 0.5, label =~ paste0(contour, " min"))all_tc_isos <- grouped_tc_coords %>%
#Using purrr to enable iteration over each transit center
mutate(isochrone_res = pmap(.l = list(lon,lat),
.f = function(lon,lat){
isochrone(tibble(lon=lon,lat=lat), contours = c(5,10,15,20),
costing = "pedestrian",hostname = nn_valhalla_hostname)
})) %>%
unnest(isochrone_res) %>%
st_as_sf()
#Generating map
leaflet() %>%
addTiles() %>%
addPolygons(data =all_tc_isos, fillColor =~ fillColor,
color="white", opacity = 0.9, weight = 2,
fillOpacity = 0.5, label =~ paste0(contour, " min")) %>%
setView(lat = 45.518875680886126, lng=-122.67939233629437, zoom = 11)