3.3 Routing

We are all generally familiar with routing algorithms through the ubiquity of Google Maps. In fact, Google Maps has an API that could be used from within R to generate routes between two points, but this is not a free service. At the moment, the most accessible free service (up to a certain limit) appears to be Mapbox, with which you should make an account. A free account (without a credit card entered) will allow you to do 100,000 directions calls and 100,000 isochrone calls per month, which is more than enough for the purposes of this chapter.

Once you’ve created an account, on the Accounts page, create a new token and select all the “secret scopes”. Once created, you will have one chance to note down the access key, which you should save on your computer.

Install mapboxapi, created by the same developer who created tigris and tidycensus. The mb_access_token() line should be run once, preferably in your console, similar to census_api_key(); in this case you should even more particularly NOT include your actual key in any published code because somebody could use up your credits. Recall that this kind of function puts the key in your .Renviron file on your computer, which you can check again after you’ve run this to confirm it worked (and then restart R or run readRenviron("~/.Renviron") for RStudio to read the key).

devtools::install_github("walkerke/mapboxapi")

library(mapboxapi)

mb_access_token("YOUR_TOKEN_HERE", install = T)

Let’s use epa_od from the last section as an input for routing, where we’ll use the centroid of each tract as the origin and destination points. There are 458 pairs in total. Below, we create two separate dataframes holding the lat/long coordinates of the origin and destination tract centroids, respectively, in the format that mapboxapi will need:

epa_od_origin <-
  epa_od %>% 
  st_drop_geometry() %>% 
  select(h_tract) %>% 
  left_join(ca_tracts %>% select(h_tract = GEOID)) %>% 
  st_as_sf() %>% 
  st_centroid() %>% 
  st_coordinates()

epa_od_destination <-
  epa_od %>% 
  st_centroid() %>% 
  st_coordinates()

Note that we had originally joined tract geometries to epa_od based on the w_tract or destination field. So when creating epa_od_origin we need to drop those destination geometries and add origin geometries, based on the h_tract field. st_centroid() and st_coordinates() get us the lat/long of the tract in a matrix format, which is what we’ll need for routing.

Now, let’s demonstrate a single route using mb_directions(). Multiple additional arguments can be provided, as described in the online documentation, but the simplest version requires just the origin/destination inputs and a profile= which can be “driving”, “driving-traffic”, “cycling”, or “walking”.

route_test <-
  mb_directions(
    origin = epa_od_origin[1, ],
    destination = epa_od_destination[1, ],
    profile = "driving-traffic"
  )
route_test %>% 
  st_as_sf() %>% 
  leaflet() %>%
  addMapboxTiles(
    style_id = "streets-v11",
    username = "mapbox"
  ) %>%
  addPolylines()

Notes:

  • addMapboxTiles() is provided through mapboxapi and can be used in the leaflet pipeline to reference publicly available basemaps from MapBox, or, as is a core feature of MapBox, basemaps you create yourself using your MapBox account. You might find these nicer looking than the other default options we’ve used in leaflet thus far, especially the satellite maps.
  • Given the distance traversed here, you might notice that the granularity of the route is not that accurate, seeming to be connecting the dots between disparate points. However, if we switch to a shorter route (EPA to Stanford), we’ll find that the route geometry is better (this is a limitation to keep in mind):
route_test <-
  mb_directions(
    origin = epa_od_origin[min(which(epa_od$w_tract == "06085513000")),],
    destination = epa_od_destination[min(which(epa_od$w_tract == "06085513000")),],
    profile = "driving-traffic"
  )

Here, we use which() to figure out the row(s) in the original epa_od dataframe where the destination is a specific tract we identified to be on Stanford campus. That index can then be fed to epa_od_origin or epa_od_destination, which only holds coordinates, to extract the correct coordinates. But since mb_directions() can only do one O-D pair at a time, we arbitrarily just keep the min() or first index we found.

route_test %>% 
  st_as_sf() %>% 
  leaflet() %>%
  addMapboxTiles(
    style_id = "streets-v11",
    username = "mapbox"
  ) %>%
  addPolylines()

Now, let’s wrap mb_directions() in map_dfr() to loop through our full epa_od dataframe:

epa_od_route <- 
  1:nrow(epa_od_origin) %>%
  map_dfr(function(x){
    mb_directions(
      origin = epa_od_origin[x, ],
      destination = epa_od_destination[x, ],
      profile = "driving-traffic"
    )
  })
epa_od_route %>% 
  st_as_sf() %>% 
  leaflet() %>%
  addMapboxTiles(
    style_id = "light-v9",
    username = "mapbox"
  ) %>%
  addPolylines()

From here, the output geometries, or the fields duration and distance, can be used for any number of analyses, including estimates of vehicle miles traveled (incorporating ACS data to estimate the % of workers in each O-D pair who are driving vs. using other modes).