knitr::opts_chunk$set(echo = TRUE)

For this tutorial, we need two packages. leaflet allows you to create an interactive map. rgdal allows you to import shapefile layers, such as transit lines or transit stops.

NOTE FOR MACOS USERS: you might need to first manually install GDAL. You can use this link for guidance https://tlocoh.r-forge.r-project.org/mac_rgeos_rgdal.html

#Load packages
library(leaflet)
library(rgdal)

Draw a map

First let’s just figure out how to get leaflet to make a map. To start, you need to determine a centerpoint latitude and longitude. To do that, use googlemaps, single click at a specific spot on the map and the lat,long values will appear on the bottom. Cut and paste those values below. I have used a point in Kendall Square. A zoom value of 12 is a good place to start.

You can choose the map background from a variety of sources. You can explore the options here: http://leaflet-extras.github.io/leaflet-providers/preview/. Move the map so you are over Boston. The default one you see is “OpenStreetMap.Mapnik”. Then on the right, you can Scroll to explore options. I chose “CartoDB.Positron” because it does not overly highlight roadways, making it easier to overlay transit lines.

# Store values into variables
myLAT <- 42.366740
myLNG <- -71.084760
myZOOM <- 12
myTILES <- "CartoDB.Positron"

# Create a map using leaflet and store that map in the variable mymap
mymap <- leaflet() %>%
  addProviderTiles(myTILES) %>%
  setView(myLNG, myLAT, zoom = myZOOM)

# Display the map by just using the variable name itself
mymap

Add subway lines

Now lets add MBTA subway lines. In the same directory as your .Rmd file, make a subdirectory called “transitlayers”. Go to the MassGIS site and download the .zip file, and extract the contents into a subfolder “mbtarapidtransit” in the “transitlayers” folder. https://docs.digital.mass.gov/dataset/massgis-data-mbta-rapid-transit

Use the readOGR command to load the shapefile into the variable called subway. You will then need to make sure that the shapefile is using the correct “projection” (one of the more confusing aspects of GIS).

Then use the leaflet function to make the map, and store the map into mymap variable. Then put the mymap on a line by itself to display the map.

# Import subway lines shape file
subway <- readOGR(dsn="transitlayers/mbtarapidtransit", layer="MBTA_ARC")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\rosen\Google Drive\CIVE4567 summer 2020\transitlayers\mbtarapidtransit", layer: "MBTA_ARC"
## with 138 features
## It has 4 fields
subway <- spTransform(subway, CRS("+proj=longlat +ellps=GRS80"))

# Create map
mymap <- leaflet() %>%
  addProviderTiles(myTILES) %>%
  setView(myLNG, myLAT, zoom = myZOOM) %>%
  addPolylines(data = subway)

# Draw map
mymap

Color the subway lines

We want to color the lines. The shapefile has some additional information about the layers, including the name (color) of each line. This information is also on the MassGIS webpage referenced above. We see that the there is a variable LINE that has the color of the subway line. We can see all the possibilities by using the levels and factors commands.

# Show all the possibilities for the variable LINE
levels(factor(subway$LINE))
## [1] "BLUE"   "GREEN"  "ORANGE" "RED"    "SILVER"

R deals with colors either by using a hex number for a color, or a number of built-in color names. This cheatsheet us useful. https://www.nceas.ucsb.edu/sites/default/files/2020-04/colorPaletteCheatsheet.pdf. We need to make a mapping of the line names to color. In this case, it happens to be that they align. We can use the colorFactor to map R colors (the palette items) with the domain (the options in the shapefile dataset). We store that mapping into the variable pal.

When we call leaflet to make the map, we set the color of the lines based on the variable LINE using the mapping function pal.

# Map colors to the different subway lines
pal <- colorFactor(palette = c('blue','green','orange','red','grey'),
                   domain = c('BLUE','GREEN','ORANGE','RED','SILVER'))

# Create the map
mymap <- leaflet() %>%
  addProviderTiles(myTILES) %>%
  setView(myLNG, myLAT, zoom = 12) %>%
  addPolylines(data = subway, color=~pal(LINE))

# Draw the map
mymap

Add stations to the map

Now lets add the stations. The shapefile for stations is not lines but points. We could choose to draw those points as circles, but we could also just draw those points as markers. We can define a new MBTA icon that can be used on the map.

# Define new MBTA icon
mbtaIcon <- makeIcon(
  iconUrl = "https://upload.wikimedia.org/wikipedia/commons/thumb/6/64/MBTA.svg/240px-MBTA.svg.png",
  iconWidth = 15, iconHeight = 15
)

# Load the MBTA stops shapefile
stops <- readOGR(dsn="transitlayers/mbtarapidtransit", layer="MBTA_NODE")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\rosen\Google Drive\CIVE4567 summer 2020\transitlayers\mbtarapidtransit", layer: "MBTA_NODE"
## with 166 features
## It has 4 fields
stops <- spTransform(stops, CRS("+proj=longlat +ellps=GRS80"))

# Create the map
mymap <- leaflet() %>%
  addProviderTiles(myTILES) %>%
  setView(myLNG, myLAT, zoom = 12) %>%
  addPolylines(data = subway, color=~pal(LINE)) %>%
  addMarkers(data = stops, icon=mbtaIcon)

# Draw the map
mymap

Add commuter rail lines to the map

Go to MassGIS and download the MBTA trains shapefiles. Put them into a subfolder called “trains” within your “transitlayers” folder. https://docs.digital.mass.gov/dataset/massgis-data-trains

This layer contains all trains, so We need to limit it to only the commuter rail. There is a variable COMMRAIL for the line layer and C_RAILSTAT for the station layer. When those variables are set to “Y” it indicates they are part of the commuter rail system. So we subset the shapefile to only include commuter rail related items.

# Import rail lines
trains <- readOGR(dsn="transitlayers/trains", layer="TRAINS_ARC")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\rosen\Google Drive\CIVE4567 summer 2020\transitlayers\trains", layer: "TRAINS_ARC"
## with 6595 features
## It has 28 fields
## Integer64 fields read as strings:  TRACK LINE_CODE ASSET_ID
trains <- spTransform(trains, CRS("+proj=longlat +ellps=GRS80"))

# Subset to only include commuter rail
trains <- subset(trains, COMMRAIL == "Y", na.rm=TRUE)

# Import rail stations
trainstops <- readOGR(dsn="transitlayers/trains", layer="TRAINS_NODE")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\rosen\Google Drive\CIVE4567 summer 2020\transitlayers\trains", layer: "TRAINS_NODE"
## with 387 features
## It has 6 fields
trainstops <- spTransform(trainstops, CRS("+proj=longlat +ellps=GRS80"))
# Subset to only include commuter rail stations
trainstops <- subset(trainstops, trainstops@data$C_RAILSTAT == "Y", na.rm=TRUE)

# Make a purple circle icon for the commuter rail stations
crIcon <- makeIcon(
  iconUrl = "https://upload.wikimedia.org/wikipedia/commons/thumb/8/80/Eo_circle_purple_circle.svg/240px-Eo_circle_purple_circle.svg.png", iconWidth = 15, iconHeight = 15)

# Create the map
mymap <- leaflet() %>%
  addProviderTiles(myTILES) %>%
  setView(myLNG, myLAT, zoom = 12) %>%
  addPolylines(data = trains, color="purple", weight=2) %>%
  addPolylines(data = subway, color=~pal(LINE)) %>%
  addMarkers(data = stops, icon=mbtaIcon) %>%
  addMarkers(data = trainstops, icon=crIcon)

# Display the map
mymap

Add bus lines to the map

Go to MassGIS and download the MBTA bus shapefiles. Put them into a subfolder called “mbtabus” within your “transitlayers” folder. https://docs.digital.mass.gov/dataset/massgis-data-mbta-bus-routes-and-stops

We want to be able to create a map that highlights the “key bus routes”. Those are high ridership and frequency routes, additional information here. https://en.wikipedia.org/wiki/MBTA_key_bus_routes. You can see them on this map. https://cdn.mbta.com/sites/default/files/2020-05/subway-map-june2020-v34a-GLX-shuttle.pdf

Here, we make the key routes yellow and slightly thicker lines with a line weight of 2, while the other bus routes are in grey with a thinner line weight of 1.

# Import bus routes shapefile
bus <- readOGR(dsn="transitlayers/mbtabus", layer="MBTABUSROUTES_ARC")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\rosen\Google Drive\CIVE4567 summer 2020\transitlayers\mbtabus", layer: "MBTABUSROUTES_ARC"
## with 923 features
## It has 11 fields
## Integer64 fields read as strings:  CTPS_ROU_2
bus <- spTransform(bus, CRS("+proj=longlat +ellps=GRS80"))

# Subset the bus routes into "key routes" and "non-key-routes"
keybuslist <- as.character(c(1,15,22,23,28,32,39,57,66,71,73,77,111,116,117))
bus.key <- subset(bus, MBTA_ROUTE %in% keybuslist, na.rm=TRUE)
bus.nokey <- subset(bus, !(MBTA_ROUTE %in% keybuslist), na.rm=TRUE)

# Create the map
mymap <- leaflet() %>%
  addProviderTiles(myTILES) %>%
  setView(myLNG, myLAT, zoom = 12) %>%
  addPolylines(data = bus.nokey, color="grey", weight=1) %>%
  addPolylines(data = bus.key, color="yellow", weight=2) %>%
  addPolylines(data = trains, color="purple", weight=2) %>%
  addPolylines(data = subway, color=~pal(LINE)) %>%
  addMarkers(data = stops, icon=mbtaIcon) %>%
  addMarkers(data = trainstops, icon=crIcon)

# Display the map
mymap

Note that in the chunk above, several variables were defined in previous chunks that are needed. When you do this on your own, combining all the pieces together into one chunk, make sure you include all the necessary building blocks and variable definitions.