In Lab 5a, we went through various methods for measuring spatial accessibility. In this lab, we complete the journey by learning how to create isochrone maps using the Mapbox platform in R. Mapbox is a mapping and location cloud platform for developers. In this lab, you will

  1. Learn about making maps using Mapbox in R
  2. Learn how to visualize spatial accessibility using isochrones
  3. Learn how to identify vulnerable neighborhoods with low access using a spatial overlay

Our case study for this lab is spatial access to hospitals in the City of Sacramento. This lab guide follows material presented in class Handout 4.

Installing and loading packages


You’ll need to install the following packages in R. You only need to do this once, so if you’ve already installed these packages, skip the code. Also, don’t put these install.packages() in your R Markdown document. Copy and paste the code in the R Console.

install.packages(c("mapboxapi",
                   "fasterize",
                   "leafsync"))

You’ll need to load the following packages. Unlike installing, you will always need to load packages whenever you start a new R session. You’ll also always need to use library() in your R Markdown file.

library(tidyverse)
library(tidycensus)
library(tigris)
library(rmapshaper)
library(sf)
library(leaflet)
library(mapboxapi)
library(fasterize)
library(leafsync)

Setting up your Mapbox account


To access Mapbox services in R, you’ll need a valid Mapbox account with an access token. To get the token, you will need to set up an account. Visit https://account.mapbox.com/auth/signup/ to establish an account - all you need to provide is an email address to sign up! Fill out the form and verify your account through the email Mapbox sends you; you’ll be taken directly to your Mapbox account dashboard page.


Nothing you’ll do today will be intensive enough to incur charges - but to access premium content, including mapping a very large file (e.g. all tracts in the United States) you will need to pay. Scroll down to the Access Tokens section. Copy the default public token that appears on your screen to your clipboard, then return to R.

Now that you have your Mapbox token in hand, you can set it up in R! Save the token into an R object

my_token <- "YOUR TOKEN GOES HERE"

And install it using the function mb_access_token(), which is a part of the mapboxapi package (Walker, 2022).

mb_access_token(my_token, install = TRUE)


The optional argument install = TRUE saves the token to your .Renviron, allowing you to use mapboxapi functions in the future without having to worry about setting your token (similar to the Census API token). To use this feature, you will have to restart your R session now. If this is giving you trouble, take this argument out of mb_access_token().

Read in data


We will need to bring in census tract polygon features and socioeconomic data using the Census API and keep tracts within Sacramento city boundaries. The code for accomplishing these tasks is below. We won’t go through each line of code in detail because we’ve covered all of these operations and functions in prior labs. I’ve embedded comments within the code that briefly explain what each chunk is doing, but go back to prior guides (or RDS/GWR) if you need further help.

# Bring in census tract data. 
ca.tracts <- get_acs(geography = "tract", 
              year = 2019,
              variables = c("B17020_002", "B17020_001","DP04_0058P"),
              state = "CA",
              survey = "acs5",
              output = "wide",
              geometry = TRUE)

# Calculate percent race/ethnicity, and keep essential vars.
ca.tracts <- ca.tracts %>% 
  mutate(ppoor = 100*(B17020_002E/B17020_001E)) %>%
  rename(pnocar = "DP04_0058PE") %>%
  dplyr::select(c(GEOID,ppoor, pnocar))  

# Bring in city boundary data
pl <- places(state = "CA", year = 2019, cb = TRUE)

# Keep Sacramento city
sac.city <- filter(pl, NAME == "Sacramento")

#Clip tracts using Sacramento boundary
sac.tracts <- ms_clip(target = ca.tracts, clip = sac.city, remove_slivers = TRUE)

Let’s bring in Hospital locations in Sacramento city, downloaded from the Health Resources and Services Administration Geospatial Data Warehouse. I zipped up the file and uploaded it onto Github. Use the following code to download and unzip the file. I also uploaded the file in Canvas (Week 5 Lab folder).

download.file(url = "https://raw.githubusercontent.com/crd230/data/master/sacramento_hospitals.zip", destfile = "sacramento_hospitals.zip")
unzip(zipfile = "sacramento_hospitals.zip")


Bring in the data using st_read().

hospitals <-st_read("sacramento_hospitals.shp")

#take a look at the data
glimpse(hospitals)

Mapbox maps in R


The mapboxapi package is an R package that interfaces with Mapbox web services APIs. Its purpose is to help R users incorporate the suite of Mapbox tools into their spatial data science projects. The most well-known feature of Mapbox services is its ability to create stunning web maps which are used on applications all around the world. Mapbox maps are accessed through styles, which are custom design configurations applied to OpenStreetMap or even user-generated vector map tilesets. Mapbox provides a number of their styles to all users with a Mapbox access token. The most recent versions of these styles are as follows:

Mapping using Mapbox in R relies on leaflet functions, which we covered in Lab 3. Let’s get a browseable Leaflet map using Mapbox tiles as a basemap. To add a Mapbox tile, use the function addMapboxTiles().

mapbox_map <- leaflet() %>%
  addMapboxTiles(style_id = "streets-v11",
                 username = "mapbox") 

mapbox_map

Visualizing spatial access


Mapbox has a lot of tools, but most pertinent to this week’s labs is its ability to visualize spatial access. Mapbox uses the Mapbox Isochrone API to draw isochrones around specified locations, which represent the reachable area from those locations within a given travel time by a given travel mode. An isochrone, from the Greek root words iso (equal) and chrone (time), is a line that connects points of equal travel time around a given location. The Mapbox Isochrone API computes areas that are reachable within a specified amount of time from a location, and returns the reachable regions as contours of polygons or lines that you can display on a map.

An isochrone is simply a type of buffer or catchment area around a point. In Lab 5a, we covered catchment area methods for measuring accessibility. There, we relied on distance buffers around centroids. The code below creates a 5km buffer around the University of California, Davis Medical Center by using the argument dist = 5000 in st_buffer().

ucd_med <- hospitals %>%
          filter(name == "UNIVERSITY OF CALIFORNIA DAVIS MEDICAL CENTER")

buf5km <- st_buffer(ucd_med, dist = 5000)

An alternative option is to create network-based isochrones, which are polygons that represent the accessible area around a given location within a given travel time for a given travel mode. Creating and visualizing isochrones is straightforward with the mb_isochrone() function. Supported travel profiles include driving (with and without traffic), cycling, and walking. mb_isochrone() by default returns a simple features polygon object that can be used for visualization. The example below draws a 10-minute driving isochrone around UCD Med for a Wednesday during evening rush hour. mb_isochrone() accepts as an input a coordinate pair, a location description as a character string, or an sf object. Here, we plug in the med center sf point object. We also could have plugged in the address as a string, “2315 Stockton Blvd, Sacramento, CA 95817”. We specify profile = "driving-traffic" to specify driving with traffic as the form of travel.

iso10min <- mb_isochrone(
  ucd_med, 
  time = 10, 
  profile = "driving-traffic",
  depart_at = "2023-01-11T17:00"
  )

We can visualize the comparative extents of these two methods in Sacramento. Run the following code to get a synced interactive map showing the two methods. The makeAwesomeIcon() function in leaflet creates a custom icon appropriate for a medical facility; many other icons are available for common points of interest.

hospital_icon <- makeAwesomeIcon(icon = "ios-medical", 
                                 markerColor = "red",
                                 library = "ion")

# The Leaflet package requires data be in CRS 4326
#buffer access map
map1 <- leaflet() %>% 
  addTiles() %>%
  addPolygons(data = st_transform(buf5km, 4326)) %>% 
  addAwesomeMarkers(data = st_transform(ucd_med, 4326),
                    icon = hospital_icon)

#10 minute travel access map
map2 <- leaflet() %>% 
  addTiles() %>%
  addPolygons(data = iso10min) %>% 
  addAwesomeMarkers(data = st_transform(ucd_med, 4326),
                    icon = hospital_icon)

sync(map1, map2)


The comparative maps illustrate the differences between the two methods quite clearly. Many areas of equal distance to the hospital do not have the same level of access. Conversely, due to the location of highways, there are some areas outside the 5km buffer area that can reach the hospital within 10 minutes.

Let’s draw multiple isochrones around the UC Davis medical center. The argument time = is a vector of isochrone contours, specified in minutes, which defaults to c(5, 10, 15). The maximum time supported is 60 minutes. Here, we create isochrones for each minute from 1 to 45 minutes.

isos <- mb_isochrone(
  location = ucd_med,
  profile = "driving-traffic",
  time = 1:45
)

isos
## Simple feature collection with 45 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -122.0552 ymin: 38.02425 xmax: -120.7421 ymax: 39.13975
## Geodetic CRS:  WGS 84
## # A tibble: 45 × 3
##     time                                                          geometry    id
##  * <int>                                                     <POLYGON [°]> <int>
##  1    45 ((-121.5898 39.13975, -121.5916 39.13769, -121.5907 39.13082, -1…     1
##  2    44 ((-121.5868 39.1317, -121.5869 39.12858, -121.5824 39.12712, -12…     1
##  3    43 ((-121.5718 39.1228, -121.5732 39.12269, -121.572 39.12253, -121…     1
##  4    42 ((-121.5568 39.11375, -121.5586 39.11369, -121.557 39.11348, -12…     1
##  5    41 ((-121.5478 39.10184, -121.5488 39.10169, -121.548 39.09869, -12…     1
##  6    40 ((-121.5418 39.08907, -121.544 39.07769, -121.5435 39.06269, -12…     1
##  7    39 ((-121.5418 39.07421, -121.5425 39.06269, -121.5453 39.06023, -1…     1
##  8    38 ((-121.5448 39.0597, -121.5449 39.05658, -121.5422 39.05631, -12…     1
##  9    37 ((-121.5418 39.0446, -121.5425 39.03569, -121.5542 39.02031, -12…     1
## 10    36 ((-121.5448 39.02989, -121.5448 39.02674, -121.5483 39.02669, -1…     1
## # … with 35 more rows

An sf object of 45 polygons is returned with a time column representing the travel-time around the location. Time is organized in descending order to ensure that overlapping isochrones are plotted correctly, with the shortest time visualized last (on top).

Using leaflet’s addPolygons() function, we can add the isochrones to our Mapbox basemap with a mostly-transparent fill opacity. We use the colorNumeric() function to cut the time variable into bins, which we color using the viridis scheme.

pal <- colorNumeric("viridis", isos$time, na.color = "transparent")

mapbox_map %>%
  addPolygons(data = isos,
              fillColor = ~pal(time),
              stroke = FALSE,
              fillOpacity = 0.1) %>%
  addLegend(values = isos$time,
            pal = pal,
            title = "Drive-time to UCD Med")

The result illustrates some of the wide differences in accessibility between various parts of the region. One notable issue with this visualization approach, however, is that the layering of isochrones in the interior of Sacramento makes it difficult to view the basemap beneath them. This can be resolved by converting to a raster dataset and generating an “accessibility surface” for improved visualization. We’ll do this in the next section.

Practice Exercise: Now that you’ve learned how to use isochrone services in mapboxapi, try it out for yourself! Create an isochrone map around a location of your choice. Times can be specified at 1-minute intervals all the way up to 60 minutes using a vector.

Making an accessibility surface


Accessibility surfaces are commonly used in geographic information systems applications to identify the distance from any particular location to a geographic feature of interest. We can apply this concept to network-based accessibility by using mapboxapi tools. To create the accessibility surface, we will convert our isochrones isos, which is an sf polygon object, to a raster object using the fasterize package. We discussed the raster view of spatial data in Handout 2, and briefly covered how to process raster data in Lab 3. Raster datasets represent geographic information as grid cells defined by a cell size. Higher-resolution raster datasets are represented with smaller cell sizes.

To generate the accessibility surface raster, we will need to reproject isos into a coordinate reference system (CRS) with meters as distance units. Meters is not the distance unit for the current CRS of isos.

st_crs(isos)$units
## NULL

We discussed coordinate reference systems in Lab 4a. Here, we reproject isos into UTM. This will allow us to specify the raster’s resolution in meters.

isos_proj <- st_transform(isos, crs = 32611)

We then generate a 100m resolution raster using the raster() function. This raster represents a “template raster” object defining the extent, resolution and CRS of the output for fasterize().

template <- raster(isos_proj, resolution = 100)

Then use the fasterize() function to allocate the minimum overlapping value from our isochrones to each grid cell. The argument fun = specifies the name of a function by which to combine overlapping polygons. The argument field = specifies the ID in isos_proj providing a value for each of the polygons rasterized.

iso_surface <- fasterize(isos_proj, template, field = "time", fun = "min")

The result can then be mapped with leaflet’s addRasterImage() function.

mapbox_map %>%
  addRasterImage(iso_surface, colors = pal, opacity = 0.5) %>%
  addLegend(values = isos$time, pal = pal,
            title = "Drive-time to UCD Med")

Comparing to the accessibility map from the prior section, accessibility is now represented in a similar way, but with a clearer view of the basemap around the UCD Medical center.

Identifying low-access neighborhoods


The previous example illustrated how to model and visualize accessibility in Sacramento; however, it does not speak directly to who may have difficulties accessing Hospitals. Earlier, we brought in tract poverty rates from the 2015-2019 American Community Survey. Our task in this section is to find neighborhoods with limited access to hospitals (not just UC Med) in Sacramento, and cross-reference this with the city’s neighborhood poverty rates. We revert back to using polygon objects in this section.

We’ll measure accessibility using isochrones as above, and consider a 10 minute drive-time around each hospital location. mb_isochrone() can accept sf objects as input, and will retain an ID from the input sf object if the column name is specified (in our case, the variable for the ID is name). Here, we will create an isochrone representing a 10-minute driving time time = 10.

driving_isos <- mb_isochrone(
  hospitals,
  profile = "driving-traffic",
  time = 10,
  id = "name"
)

These results can be visualized on our Mapbox map:

mapbox_map %>%
  addPolygons(data = driving_isos,
              popup = ~id)

The map represents the reachable area to a hospital within a 10-minute drive, modeled at average driving speed assuming no traffic.


Practice Question: Which spatial data pitfall discussed in Handout 2 might we run into with this analysis?


We want to find the neighborhoods that have low access to hospitals and contains a vulnerable population, which we measure using the poverty rate. Let’s visualize neighborhood poverty in Sacramento on our Mapbox map:

driving_pal <- colorNumeric("viridis", sac.tracts$ppoor)

mapbox_map %>%
  addPolygons(data = sac.tracts,
              fillColor = ~driving_pal(ppoor),
              fillOpacity = 0.5,
              stroke = FALSE,
              smoothFactor = 0.1,
              label = ~round(ppoor, 1)) %>%
  addLegend(values = sac.tracts$ppoor,
            pal = driving_pal,
            title = "% poor")

Next, let’s conduct a spatial overlay of our hospital accessibility layer on Sacramento neighborhood poverty. Spatial overlay is a very common operation when working with spatial data. It can be used to determine which features in one spatial layer overlap with another spatial layer, or extract data from a layer based on geographic information.

In our example, we want to determine the areas in Sacramento city with the greatest proportion of poor households that also are beyond a 10 minute drive from a hospital. To do this, we use the following steps:

  1. We transform the coordinate reference system of our sac.tracts dataset to 4326, the same CRS used by the isochrones;
st_crs(driving_isos)$epsg
## [1] 4326
  1. We extract only those Census tracts with a percentage of poor households that is 30 percent or above;
  2. We use the st_difference() and st_union() functions from the sf package to “cut out” or clip areas from those Census tracts that overlap the 10-minute driving isochrones.

Let’s complete these steps

low_access_dareas <- sac.tracts %>%
    st_transform(crs = 4326) %>%
  filter(ppoor >= 30) %>%
  st_difference(
    st_union(driving_isos)
  )

We can visualize the result on our Mapbox map.

mapbox_map %>%
  addPolygons(data = low_access_dareas)

As the map illustrates, there are several areas that are located beyond a 10-minute drive from a hospital and are high poverty. Notable clusters of neighborhoods that meet this criteria are located in North Sacramento and near Del Paso. Granted, this analysis is not definitive, but gives us some insights into potential issues with hospital accessibility and how we might resolve them.

Guess what? You earned a badge! [Congratulations!](https://www.youtube.com/watch?v=04854XqcfCY

Creative Commons License
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

Website created and maintained by Noli Brazil