Introduction

This exercise compares isochrones created by Open Source Routing Machine (OSRM) with different specified resolutions. The isochrones were built from the starting coordinates of City Hall in Philadelphia, PA, and include catchment areas for a 0 to 5 minute drive, a 5 to 10 minute drive, and a 10 to 15 minute drive from City Hall. The effectiveness of the isochrones created by each service were analyzed with calculated drive times to 10 points of interest that were either within the city of Philadelphia or the neighboring city, Camden, NJ. Drive times between City Hall and each point of interest were calculated with the Bing Maps Routing API during a time that should have light to moderate traffic.

Step 1: Install and load libraries

# Use pacman package to load packages used for exercise
library(pacman)

p_load(dplyr, rjson, lubridate, kableExtra, httr, geojsonR, sf, sp, osrm, gtools, slickR)

Step 2: Load points

# Load Points from working directory
points <- read.csv("./data/Philadelphia-Points.csv")

Step 3: Calculate drive times

The first point in the data set was City Hall in Philadelphia, which was used for creating the isochrones and was the starting point for drive times to each point of interest. All other points in the data set were for the other destinations in Philadelphia or Camden. The travel time in minutes for a car was calculated for each point of interest as the end point, calculated using traffic conditions at 10AM on 11/02/2021.

Bing maps API key available here.

# Calculate drive time in minutes between first point (City Hall) and other points
for (i in 1:nrow(points)){
  url <- paste("http://dev.virtualearth.net/REST/V1/Routes/",
               "Driving",
               "?wp.0=",paste0(points$lat[1],"+",points$lon[1]),
               "&wp.1=",paste0(points$lat[i],"+",points$lon[i]),
               "&dt=",paste0(format(ymd(Sys.Date()), "%m/%d/%Y"),"+10:00:00"),
               "&tt=","Departure",
               "&key=",key,
               sep="")
  
  suppressWarnings(json_bing <- fromJSON(paste(readLines(url), collapse=""), 
                                         unexpected.escape = "skip"))
  
  travelDuration <- json_bing$resourceSets[[1]]$resources[[1]]$travelDuration
  
  points$time[i] <- round(travelDuration/60,1)
}

# Add correct catchment area to points
points$iso_goal <- ifelse(points$time <= 5, 5,
                   ifelse(points$time <= 10, 10, 15))

# Display points data (only relevant columns)
points[-1,c(2:3,7:8)] %>%
  kbl(caption = "Points of Interest with Drive Time in Minutes from 
      Philadelphia City Hall and Catchment Area Goal", align = "r",
      row.names = FALSE) %>%
  kable_minimal("hover")
Points of Interest with Drive Time in Minutes from Philadelphia City Hall and Catchment Area Goal
location address time iso_goal
HipCityVeg 214 S 40th St, Philadelphia, PA 19103 12.4 15
20th Street Pizza 108 S 20th St, Philadelphia, PA 19103 6.3 10
Grindcore House West 4134 Chester Ave, Philadelphia, PA 19104 11.7 15
Dottie’s Donuts 509 S 6th St, Philadelphia, PA 19147 7.8 10
Grindcore House 1515 S 4th St, Philadelphia, PA 19147 9.5 10
Blackbird 614 N 2nd St, Philadelphia, PA 19123 9.7 10
Rodin Museum 2151 Benjamin Franklin Pkwy, Philadelphia, PA 19130 5.1 10
RCA Pier Water St, Camden, NJ 08103 10.1 15
Fresh Mex 1 Market St, Camden, NJ 08102 8.2 10
Hello Donuts + Coffee 2557 Amber St, Philadelphia, PA 19125 13.3 15

Step 4: Create spatial points

Create spatial points using the coordinates of the destinations in the data set. These will be used for evaluating if each isochrone generated by OSRM is accurate or inaccurate with their catchment areas.

# Create spatial points data frame for points of interest (exclude starting location)
coords    <- points[2:11,c(5,6)]
points.sf <- st_as_sf(x = points[2:11,], coords = c("lon", "lat"),
                      crs = "+proj=longlat +datum=WGS84 +no_defs")

points.sp <- as(points.sf, "Spatial")

Step 5: Isochrone building function

# Create table for storing build time
times <- data.frame(resolution = c(seq(from = 3, to = 9, by = 1), seq(from = 10, to = 150, by = 10)),
                    time = NA, polygons = NA)

# OSRM isochrone function
osrm_fun <- function(res){
  
  # Track how long it takes to build the isochrone
  start.time <- Sys.time()
  
  # Download data from OSRM API
  osrm <- osrmIsochrone(loc = c(points$lon[1], points$lat[1]),
                        breaks = seq(from = 0,to = 15, by = 5), res = res)
  
  # Add information for graphics and further analyses
  osrm$name <- osrm$max

  # Stop timer for isochrone creation
  stop.time                         <- Sys.time()
  times[times$resolution == res,2] <<- round(difftime(stop.time, start.time, 
                                                        units = c("mins")), 2)
  
  # Store number of polygons for each isochrone
  times[time$reoslution == res,3] <<- length(osrm$name)

  # Save OSRM isochrone
  return(osrm)
  
}

Step 6: Build isochrones

# Run osrm_fun from Step 5 for resolutions 3 through 10, and then from 20 to 150 by +10
osrm3   <- osrm_fun(3)
osrm4   <- osrm_fun(4)
osrm5   <- osrm_fun(5)
[...]
osrm150 <- osrm_fun(150)

Step 7: Display isochrones

# Times to build each isochrone (and how many polygons were built)
times %>%
  kbl(format = "html", escape = F, align = "r",
      caption = "Time in Minutes to Generate OSRM Isochrone by Resolution") %>%
  kable_minimal("hover", full_width = F, position = "left") %>%
  column_spec(column = 1:3, width = "2in")
Time in Minutes to Generate OSRM Isochrone by Resolution
resolution time polygons
3 0.01 3
4 0.10 1
5 0.10 3
6 0.10 1
7 0.10 3
8 0.10 2
9 0.10 3
10 0.10 2
20 0.12 3
30 0.37 3
40 0.65 3
50 0.85 3
60 1.48 3
70 1.89 3
80 3.51 3
90 3.00 3
100 4.90 3
110 4.27 3
120 5.25 3
130 7.71 3
140 7.63 3
150 22.40 3 
# Low resolution isochrones
par(mfrow = c(2,4))
plot(osrm3,  main = "Resolution = 3")
plot(osrm4,  main = "Resolution = 4")
plot(osrm5,  main = "Resolution = 5")
plot(osrm6,  main = "Resolution = 6")
plot(osrm7,  main = "Resolution = 7")
plot(osrm8,  main = "Resolution = 8")
plot(osrm9,  main = "Resolution = 9")
plot(osrm10, main = "Resolution = 10")

# Medium resolution isochrones
par(mfrow = c(2,4))
plot(osrm20, main = "Resolution = 20")
plot(osrm30, main = "Resolution = 30")
plot(osrm40, main = "Resolution = 40")
plot(osrm50, main = "Resolution = 50")
plot(osrm60, main = "Resolution = 60")
plot(osrm70, main = "Resolution = 70")
plot(osrm80, main = "Resolution = 80")
plot(osrm90, main = "Resolution = 90")

# High resolution isochrones
par(mfrow = c(2,3))
plot(osrm100, main = "Resolution = 100")
plot(osrm100, main = "Resolution = 110")
plot(osrm100, main = "Resolution = 120")
plot(osrm100, main = "Resolution = 130")
plot(osrm100, main = "Resolution = 140")
plot(osrm100, main = "Resolution = 150")

Step 8: Add catchment areas

This step added the name of each catchment area that a destination was located within for each isochrone created by OSRM. The best case scenario for each destination was that it’s calcaulted drive time from City Hall was within the range of the catchment area for the isochrone. It was also okay if the point of interest was placed in the next catchment area out (i.e. it’s okay if a destination with a 9.5 minute drive time ended up in the 10-15 minute drive time catchment area for an isochrone). The concerning outcome occurred whenever an isochrone placed a destination in a catchment area that’s maximum drive time was less than the calculated drive time to that destination.

points.sp$osrm3   <- over(points.sp, osrm3)$name
points.sp$osrm4   <- over(points.sp, osrm4)$name
points.sp$osrm5   <- over(points.sp, osrm5)$name
[...]
points.sp$osrm150 <- over(points.sp, osrm150)$name
Destination Catchment Areas by OSRM Isocrhone Resolution
location time iso_goal osrm3 osrm4 osrm5 osrm6 osrm7 osrm8 osrm9 osrm10
HipCityVeg 12.4 15 5 NA 5 15 5 NA 5 15
20th Street Pizza 6.3 10 5 NA 5 15 5 NA 5 15
Grindcore House West 11.7 15 5 NA 10 15 5 NA 10 15
Dottie’s Donuts 7.8 10 5 NA 5 15 5 NA 5 15
Grindcore House 9.5 10 5 NA 5 15 5 NA 5 15
Blackbird 9.7 10 5 NA 5 15 5 NA 5 15
Rodin Museum 5.1 10 5 NA 5 15 5 NA 5 15
RCA Pier 10.1 15 5 NA 5 15 5 NA 5 15
Fresh Mex 8.2 10 5 NA 5 15 5 NA 5 15
Hello Donuts + Coffee 13.3 15 5 NA 10 15 5 NA 10 15
location time iso_goal osrm20 osrm30 osrm40 osrm50 osrm60 osrm70 osrm80 osrm90
HipCityVeg 12.4 15 10 10 10 10 10 10 10 10
20th Street Pizza 6.3 10 5 5 5 5 5 5 5 5
Grindcore House West 11.7 15 10 10 10 10 10 10 10 10
Dottie’s Donuts 7.8 10 10 5 10 10 5 5 5 10
Grindcore House 9.5 10 10 10 10 10 10 10 10 10
Blackbird 9.7 10 5 5 5 5 5 5 5 10
Rodin Museum 5.1 10 5 5 5 5 5 5 5 5
RCA Pier 10.1 15 10 10 10 10 10 10 10 10
Fresh Mex 8.2 10 10 10 10 10 10 10 10 10
Hello Donuts + Coffee 13.3 15 10 10 10 10 10 10 10 10
location time iso_goal osrm100 osrm110 osrm120 osrm130 osrm140 osrm150
HipCityVeg 12.4 15 10 10 10 10 10 10
20th Street Pizza 6.3 10 5 5 5 5 5 5
Grindcore House West 11.7 15 10 10 10 10 10 10
Dottie’s Donuts 7.8 10 5 10 5 5 5 10
Grindcore House 9.5 10 10 10 10 10 10 10
Blackbird 9.7 10 5 5 5 5 10 5
Rodin Museum 5.1 10 5 5 5 5 5 5
RCA Pier 10.1 15 10 10 10 15 10 15
Fresh Mex 8.2 10 10 10 10 10 10 10
Hello Donuts + Coffee 13.3 15 10 10 10 10 10 10

Note: Cells in red indicate isochrones that are too generous with their catchment areas.

Isochrones created by OSRM’s API appeared to plateau somewhere between 60 and 90 (on the higher end of the mid resolutions above), with a diminishing return to improved polygon shape at the cost of time in the higher resolution isochrones above.

Step 9: Compare OSRM isochrones

# Function to plot lower resolution isochrones over resolution = 150 isochrone
o.plot <- function(osrm, res){
  png(file = paste0("./img/test",res,".png"), height = 600)
  plot(osrm150, col = "red", lwd = 2, 
       main = paste0("Resolution 150 in Red\nResolution ",res," in Blue"),
       sub = "(Points of Interest are Yellow +)",
       cex.main = 2, cex.sub = 2)
  plot(osrm,   col = rgb(red = 0, green = 0, blue = 1, alpha = .5), add = TRUE)
  points(points.sp, pch = 3, cex = 1.25, lwd = 3, col = "yellow")
  dev.off()
}
# Create and store images for isochrone comparisons
o.plot(osrm3, 3)
o.plot(osrm4, 4)
o.plot(osrm5, 5)
[...]
o.plot(osrm150, 150)

# Load list of images with isochrones
imgs <- list.files("./img/", pattern = "png", full.names = TRUE)
imgs <- mixedsort(sort(imgs))

# Create image carousel
slickR(imgs, height = 600)

Conclusion

Isochrones created by OSRM with a resolution below 20 or 30 do not seem useful. Higher resolution isochrones (resolution above 90) and even some of the medium resolution isochrones (resolution between 60 and 90) appeared to create somewhat stable polygons that refined as the resolution increased. However, those refinements came at the cost of rapidly increasing creation times. If a project required the development of many isochrones, then going above a resolution of 60 or 70 would add significantly more computing time with possibly very little difference in the placement of points in catchment areas. Unfortunately, even with a resolution above 100, OSRM’s isochrones created catchment areas that overestimated how far a person could drive from the starting location, which was more problematic than underestimating how far a person could drive, especially in a case like analyzing access to a medical professional or distribution site.