Introduction

This exercise compares isochrones created by Mapbox and Open Source Routing Machine (OSRM). The isochrones were built from the 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.

Step 1: Install and load libraries

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

p_load(rjson, lubridate, kableExtra, httr, geojsonR, sp, osrm, leaflet, sf)

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.

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: Build Mapbox isochrone

Mapbox token available here.

# Track how long it takes to build the isochrone
start.time <- Sys.time()

# Compile coordinates for the starting point (City Hall)
coord <- paste0(points$lon[1],",",points$lat[1])

# Create catchment area breaks (drive time minutes)
sub1 <- "5"
sub2 <- "10"
sub3 <- "15"
time <- paste0(sub1,",",sub2,",",sub3)

# Download data from Mapbox API
iso.url <- paste("https://api.mapbox.com/isochrone/v1/mapbox/driving/",coord,
                 "?contours_minutes=",time,"&polygons=true&access_token=",token,sep = "")

# Compile each individual catchment area polygon
r         <- GET(url = iso.url)
r.header  <- headers(r)
rr        <- content(r, "text")
rrr       <- FROM_GeoJson(rr)
sr1       <- Polygon(rrr[["features"]][[1]][["geometry"]][["coordinates"]], hole = TRUE)
srs1      <- Polygons(list(sr1), sub1)
sr2       <- Polygon(rrr[["features"]][[2]][["geometry"]][["coordinates"]], hole = TRUE)
srs2      <- Polygons(list(sr2), sub2)
sr3       <- Polygon(rrr[["features"]][[3]][["geometry"]][["coordinates"]], hole = TRUE)
srs3      <- Polygons(list(sr3), sub3)

# Combine catchment area polygons into one shapefile
poly    <- SpatialPolygons(list(srs3, srs2, srs1), 1:3)
mapbox  <- SpatialPolygonsDataFrame(poly,data = data.frame(name = c(3,2,1),
           row.names = row.names(poly)))

# Add projection
proj4string(mapbox) <- "+proj=longlat +datum=WGS84 +no_defs"

# Add information for graphics and further analyses
mapbox@data$name        <- c(sub1, sub2, sub3)
mapbox@data$id          <- 1:nrow(mapbox@data)
mapbox@data$drive_times <- factor(paste0(mapbox@data$id,
                           ". ", as.numeric(mapbox@data$name) - 5,
                           " to ", mapbox@data$name, " min"))

# Stop timer for isochrone creation
stop.time     <- Sys.time()
mb.time.diff  <- round(difftime(stop.time, start.time, units = c("mins")), 2)

# Create a color palette for the Mapbox catchment areas
mb.pal <- colorFactor(rev(topo.colors(3)), mapbox@data$drive_times)

# Remove data that is no longer needed
rm(coord, sub1, sub2, sub3, time, iso.url, r, r.header, rr, rrr, 
   sr1, srs1, sr2, srs2, sr3, srs3, poly, start.time, stop.time)

# Check isochrone polygons
plot(mapbox)

Step 5: Build OSRM isochrone

There is no API Key required to build isochrones with OSRM.

# Track how long it takes to build the isochrone
start.time <- Sys.time()

# Set resolution to a number between 30 and 100 (30 is default)
res <- 70

# 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,
                      returnclass = "sp")

# Add information for graphics and further analyses
osrm$name             <- osrm$max
osrm@data$drive_times <- factor(paste0(osrm@data$id,
                         ". ", osrm@data$min,
                         " to ", osrm@data$max, " min"))

# Stop timer for isochrone creation
stop.time       <- Sys.time()
osrm.time.diff  <- round(difftime(stop.time, start.time, units = c("mins")), 2)

# Create a color palette for the OSRM catchment areas
osrm.pal <- colorFactor(rev(heat.colors(3)), osrm@data$drive_times)

# Remove data that is no longer needed
rm(start.time, stop.time)

#Check OSRM polygons
plot(osrm)

Step 6: Mapbox isochrone map

In the interactive map below, created with the Leaflet package for R, select a point of interest on the map to see additional details for that destination. These details included the calculated drive time from City Hall from Step 2 above and the street address for the point of interest.

#Plot mapbox isochrone with points of interest
leaflet() %>% 
  setView(points$lon[1], points$lat[1], zoom = 11) %>%
  addProviderTiles("CartoDB.Positron", group="Greyscale") %>% 
  addMarkers(points$lon[-1], points$lat[-1],
             popup = paste0("<b><a href = '",points$website[-1],
             "'>",points$location[-1],"</a></b>",
             "<br>", points$address[-1],
             "<br>", "Drive time: ",
             points$time[-1], " minutes")) %>% 
  addPolygons(fill=TRUE, stroke=TRUE, color = "black",
              fillColor = ~mb.pal(mapbox@data$drive_times),
              weight=0.5, fillOpacity=0.2,
              data = mapbox, popup = mapbox@data$drive_times,
              group = "Drive Time") %>% 
  
# Add a legend
addLegend("bottomleft", pal = mb.pal, values = mapbox@data$drive_time, opacity=0.2,
          title = "Mapbox Isochrones")

Step 7: OSRM isochrone map

The interactive map below behaves like the graphic in Step 6 above, but displays the isochrone created by OSRM.

# Plot OSRM isochrone with points of interest
leaflet() %>% 
  setView(points$lon[1], points$lat[1], zoom = 11) %>%
  addProviderTiles("CartoDB.Positron", group="Greyscale") %>% 
  addMarkers(points$lon[-1], points$lat[-1],
             popup = paste0("<b><a href = '",points$website[-1],
             "'>",points$location[-1],"</a></b>",
             "<br>", points$address[-1],
             "<br>", "Drive time: ",
             points$time[-1], " minutes")) %>% 
  addPolygons(fill=TRUE, stroke=TRUE, color = "black",
              fillColor = ~osrm.pal(osrm@data$drive_times),
              weight=0.5, fillOpacity=0.2,
              data = osrm, popup = osrm@data$drive_times,
              group = "Drive Time") %>% 

# Add a legend
addLegend("bottomleft", pal = osrm.pal, values = osrm@data$drive_time, opacity=0.2,
          title = paste0("OSRM Isochrones (Resolution = ",res,")"))

Step 8: Identify catchment areas

Each of the 10 destinations were spatially joined with the Mapbox and OSRM isochrones to identify which catchment area they were located within for each isochrone. Then, using the calculated drive times, each point of interest was noted as either being correctly located by the Mapbox isochrone, the OSRM isochrone, both isochrones, or neither isochrone in the variable “accurate” below. If the drive time was less than the maximum time of the catchment area range assigned to that destination for an isochrone, then that destination was considered correctly located by that isochrone. However, if the drive time for that destination was greater than the maximum time of the catchment area range assigned to that destination for an isochrone, then that destination was considered incorrectly located by that isochrone.

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

points.sp                       <- as(points.sf, "Spatial")
points.sp@proj4string@projargs  <- mapbox@proj4string@projargs

# Add catchment areas from Mapbox and OSRM isochrones
points.sp@data$mapbox <- as.numeric(over(points.sp, mapbox)$name)
points.sp$osrm        <- as.numeric(over(points.sp, osrm)$name)

# Identify which isochrones were accurate (or both or neither)
points.sp$accurate <- ifelse(points.sp$mapbox == points.sp$iso_goal & 
                             points.sp$osrm == points.sp$iso_goal,"both",
                      ifelse(points.sp$mapbox == points.sp$iso_goal &
                             points.sp$osrm != points.sp$iso_goal,"mapbox",
                      ifelse(points.sp$mapbox != points.sp$iso_goal &
                             points.sp$osrm == points.sp$iso_goal,"osrm", "neither")))

# Display results
points.sp@data[-1,-c(1,3,4)] %>%
  kbl(caption = "Catchment Areas for Each Location by Isochrone API Service",
      align = "r", row.names = FALSE) %>%
  kable_minimal("hover")
Catchment Areas for Each Location by Isochrone API Service
location time iso_goal mapbox osrm accurate
HipCityVeg 12.4 15 15 10 mapbox
20th Street Pizza 6.3 10 10 5 mapbox
Grindcore House West 11.7 15 15 10 mapbox
Dottie’s Donuts 7.8 10 10 5 mapbox
Grindcore House 9.5 10 15 10 osrm
Blackbird 9.7 10 10 5 mapbox
Rodin Museum 5.1 10 10 5 mapbox
RCA Pier 10.1 15 15 10 mapbox
Fresh Mex 8.2 10 10 10 both
Hello Donuts + Coffee 13.3 15 15 10 mapbox

Note: Mapbox included the 6th location “Grindcore House” within its 10-15 minute catchment area, but Bing Maps calculated the drive time as 9.5 minutes. This is technically incorrect, but a disagreement by about 30 seconds seems trivial for most applications of catchment areas, especially if it is an overestimate for travel time.

Step 9: Results

# Remove City Hall from points data for analysis
for_counts <- points.sp@data[-1,]

# Count of points that were correctly identified by Mapbox or OSRM or neither
counts              <- data.frame(dplyr::count(for_counts, for_counts$accurate))
colnames(counts)[1] <- "iso"

mapbox.c <- as.numeric(sum(counts[counts$iso %in% c("both", "mapbox"),"n"]))
if(length(mapbox.c) == 0){mapbox.c <- 0}

osrm.c <- as.numeric(sum(counts[counts$iso %in% c("both", "osrm"), "n"]))
if(length(osrm.c) == 0){osrm.c <- 0}

neither.c <- nrow(for_counts) - max(mapbox.c, osrm.c)

# Plot isochrones over each other
plot(osrm, col = rgb(red = 1, green = 0, blue = 0), lwd = 3,
     main = "Isochrone Comparison Plot:\nOSRM in Red, Mapbox in Blue",
     sub = "(Points of Interest are Yellow +)")

plot(mapbox, col = rgb(red = 0, green = 0, blue = 1, alpha = 0.75),
     add = TRUE, lwd = 3)

points(points.sp, pch = 3, cex = 1, lwd = 3, col = "yellow")


Catchment area accuracy

Mapbox correctly identified the drive time catchment area for9 of 10 destinations.

OSRM correctly identified the drive time catchment area for 2 of 10 destinations.

There was 1 destination that was not correctly identified by either method.


Time to calulate each isochrone:

Mapbox: 0.05 minutes

OSRM: 6.69 minutes

Conclusion

Mapbox appeared to generate more accurate catchment areas than OSRM when compared to Bing Maps drive times. Additionally, the Mapbox isochrone was created in a fraction of the time needed to create the OSRM isochrone. Both methods included bodies of water within their catchment area, but OSRM had more problems with that aspect. Mapbox did not include Petty Island–an island East of Philadelphia in the Delaware River that is only accessible from the Camden, NJ side–in its catchment areas, which was correct. However, OSRM included chunks of that island in its catchment areas, which does not seem feasible even with empty roads in less than 15 minutes. Mapbox created more accurate catchment areas in this analysis, and may be the more reliable service to use if the analysis fits Mapbox’s paramaters: four or less catchment areas with a maximum drive time of 60 minutes. OSRM can create more catchment areas at once, and can go a further maximum drive time than Mapbox, which makes it useful for different purposes, though with the caveat that it appears to create generous catchment areas.