library(tigris)
library(acs)
library(leaflet)
library(mapview)
library(ggplot2)
library(maptools)
library(stringr)
library(rgdal)

Load Shapefile of Census Tracts

countylist <- c('17','21','25', '09')
shapefile <- tracts(state='MA', county=countylist, year=2018, class='sp')
## Warning in proj4string(obj): CRS object has comment, which is lost in output

Import MBTA Transit Lines

# Import subway lines shape file
subway <- readOGR(dsn="mbtarapidtransit", layer="MBTA_ARC")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/rachaelburstein/Desktop/R/mbtarapidtransit", layer: "MBTA_ARC"
## with 138 features
## It has 4 fields
subway <- spTransform(subway, CRS("+proj=longlat +ellps=GRS80"))

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

# 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="mbtarapidtransit", layer="MBTA_NODE")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/rachaelburstein/Desktop/R/mbtarapidtransit", layer: "MBTA_NODE"
## with 166 features
## It has 4 fields
stops <- spTransform(stops, CRS("+proj=longlat +ellps=GRS80"))


# Import rail lines
trains <- readOGR(dsn="trains", layer="TRAINS_ARC")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/rachaelburstein/Desktop/R/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="trains", layer="TRAINS_NODE")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/rachaelburstein/Desktop/R/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)


# Import bus routes shapefile
bus <- readOGR(dsn="mbtabus", layer="MBTABUSROUTES_ARC")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/rachaelburstein/Desktop/R/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)

Loading Census Data

# Span
myspan <- 5

# End year
myendyear <- 2018

# Which table
mytable <- acs.lookup(endyear=2018, table.number="B08301")
## Warning in acs.lookup(endyear = 2018, table.number = "B08301"): acs.lookup for endyear>2015: using 2015 variable codes to access 2018 data.
##   (See ?acs.lookup for details)
## Warning in acs.lookup(endyear = 2018, table.number = "B08301"): temporarily downloading and using archived XML variable lookup files;
##   since this is *much* slower, recommend running
##   acs.tables.install()
# Limit the geography
# Convert the list of counties that is in text form to numeric form
countylist2 <- as.numeric(countylist)
# Define the geography to limit the data import
mygeo <- geo.make(state=25, county=countylist2, tract="*")



myvars <- mytable[1]+mytable[10] 
mydata <- acs.fetch(endyear=myendyear, span=myspan, geography=mygeo, variable=myvars)
mystate <- as.character(mydata@geography$state)

mycounty <- as.character(mydata@geography$county)

# add leading zeros to make the total number of characters 3
mycountypadded <- str_pad(mycounty, 3, pad="0")

mytract <- as.character(mydata@geography$tract)

acsgeoid <- paste0(mystate,mycountypadded,mytract)
# Put the dataset into a dataframe format
mydatadf <- data.frame(acsgeoid, mydata@estimate)
# Rename column names
colnames(mydatadf)=c("GEOID", "total", "transit")
# Create new column and calculate the percent commute trips by transit
mydatadf$transitpercent <- mydatadf$transit/mydatadf$total*100
# Merge the ACS data with the shapefile
mydatamerged <- geo_join(shapefile, mydatadf, "GEOID", "GEOID")

# There is some weirdness with some tracts showing up as 100%, so eliminate those by subsetting the results to excluse them.
mydatamerged <- subset(mydatamerged, mydatamerged@data$transitpercent<100)
mypopup <- paste0("GEOID: ", mydatamerged$GEOID, "<br>", "Commute by Public transit: ", round(mydatamerged$transitpercent,1), "%")
mypal <- colorNumeric(
  palette = "YlGnBu",
  domain = NULL
)
# Store values into variables
myLAT <- 42.366740
myLNG <- -71.084760
myZOOM <- 12
myTILES <- "CartoDB.Positron"
mymap <- leaflet() %>%
  addProviderTiles(myTILES) %>%
  setView(myLNG, myLAT, zoom = myZOOM) %>%
  addPolygons(data = mydatamerged, 
              fillColor = ~mypal(transitpercent), # this makes the fill color on the spectrum of the palette
              color = "#b2aeae", # this is the color of the outline of the shapes
              fillOpacity = 0.7, # how see through the shapes are
              weight = 1, 
              smoothFactor = 0.2,
              popup = mypopup) %>%
  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) %>%
  addLegend(pal = mypal, 
            values = mydatamerged$transitpercent, 
            position = "bottomright", 
            title = "Commute by Public Transit",
            labFormat = labelFormat(suffix = "%"),
            bins=5)
# Display the map
mymap