What’s inside atl object

Now, let’s take a look the atl object in which we read the GTFA feed from MARTA. This object is a list. In it, names(atl) shows that there are 9 data.frames. Notice that there are many tables displayed in Table 1 above, but only some of them are in atl. Notice in Table 1 that only some are required; others are optional.

typeof(atl)

## [1] "list"

names(atl)

## [1] "agency"         "calendar"       "calendar_dates" "routes"        
## [5] "shapes"         "stop_times"     "stops"          "trips"         
## [9] "."

print(head(atl))

## $agency
## # A tibble: 1 × 7
##   agency_id agency_name                  agenc…¹ agenc…² agenc…³ agenc…⁴ agenc…⁵
##   <chr>     <chr>                        <chr>   <chr>   <chr>   <chr>   <chr>  
## 1 MARTA     Metropolitan Atlanta Rapid … http:/… Americ… en      (404)8… custse…
## # … with abbreviated variable names ¹​agency_url, ²​agency_timezone,
## #   ³​agency_lang, ⁴​agency_phone, ⁵​agency_email
## 
## $calendar
## # A tibble: 4 × 10
##   service_id monday tuesday wednesday thursday friday saturday sunday start_date
##   <chr>       <int>   <int>     <int>    <int>  <int>    <int>  <int> <date>    
## 1 2               0       0         0        0      0        0      0 2021-08-14
## 2 3               0       0         0        0      0        1      0 2021-08-14
## 3 4               0       0         0        0      0        0      1 2021-08-14
## 4 5               1       1         1        1      1        0      0 2021-08-14
## # … with 1 more variable: end_date <date>
## # ℹ Use `colnames()` to see all variable names
## 
## $calendar_dates
## # A tibble: 3 × 3
##   service_id date       exception_type
##   <chr>      <date>              <int>
## 1 5          2021-09-06              2
## 2 5          2021-11-25              2
## 3 5          2021-11-26              2
## 
## $routes
## # A tibble: 118 × 8
##    route_id route_short_name route_lon…¹ route…² route…³ route…⁴ route…⁵ route…⁶
##    <chr>    <chr>            <chr>       <chr>     <int> <chr>   <chr>   <chr>  
##  1 15664    1                Marietta B… ""            3 https:… FF00FF  ""     
##  2 15665    2                Ponce de L… ""            3 https:… 008000  ""     
##  3 15666    3                Martin Lut… ""            3 https:… FF8000  ""     
##  4 15667    4                Moreland A… ""            3 https:… FF00FF  ""     
##  5 15668    5                Piedmont R… ""            3 https:… 00FFFF  ""     
##  6 15669    6                Clifton Ro… ""            3 https:… 008080  ""     
##  7 15670    8                North Drui… ""            3 https:… 008000  ""     
##  8 15671    9                Boulevard … ""            3 https:… 00FFFF  ""     
##  9 15672    12               Howell Mil… ""            3 https:… FF00FF  ""     
## 10 15673    14               14th Stree… ""            3 https:… 00FF00  ""     
## # … with 108 more rows, and abbreviated variable names ¹​route_long_name,
## #   ²​route_desc, ³​route_type, ⁴​route_url, ⁵​route_color, ⁶​route_text_color
## # ℹ Use `print(n = ...)` to see more rows
## 
## $shapes
## # A tibble: 370,776 × 4
##    shape_id shape_pt_lat shape_pt_lon shape_pt_sequence
##    <chr>           <dbl>        <dbl>             <int>
##  1 119402           33.8        -84.6                 1
##  2 119402           33.8        -84.6                 2
##  3 119402           33.8        -84.6                 3
##  4 119402           33.8        -84.6                 4
##  5 119402           33.8        -84.6                 5
##  6 119402           33.8        -84.6                 6
##  7 119402           33.8        -84.6                 7
##  8 119402           33.8        -84.6                 8
##  9 119402           33.8        -84.6                 9
## 10 119402           33.8        -84.6                10
## # … with 370,766 more rows
## # ℹ Use `print(n = ...)` to see more rows
## 
## $stop_times
## # A tibble: 1,033,990 × 5
##    trip_id arrival_time departure_time stop_id stop_sequence
##    <chr>   <time>       <time>         <chr>           <int>
##  1 8044645 14:35:00     14:35:00       903320              1
##  2 8044645 14:36:02     14:36:02       903448              2
##  3 8044645 14:36:51     14:36:51       901144              3
##  4 8044645 14:37:43     14:37:43       904219              4
##  5 8044645 14:38:19     14:38:19       903850              5
##  6 8044645 14:38:39     14:38:39       903664              6
##  7 8044645 14:39:39     14:39:39       213268              7
##  8 8044645 14:40:12     14:40:12       213101              8
##  9 8044645 14:41:35     14:41:35       212926              9
## 10 8044645 14:41:53     14:41:53       211886             10
## # … with 1,033,980 more rows
## # ℹ Use `print(n = ...)` to see more rows

Step 1: Extracting lines

GTFS from MARTA contains information about various transit types, including rails and buses. The routes table in atl has a column named route_type, which contains integer values denoting types of transit service (e.g., bus vs. rail transit). To see what each integer means, you will need to look at the link to the Reference provided above (or see below for a screen capture from the Reference).

Reference from Google

routes table contains route_types, which is needed for separating rail transit from bus transit.shapes table contains the geographic information needed for geo-operations such as buffer, intersections, etc. Because routes and shapes table do not contain a common key, they need to be joined through trips table as intermediate table.

Notice that there can be multiple trips that runs on the same routes. If we join trips and routes, one rows in routes table will be matched with multiple rows in trips table and create duplicated.

# Join routes table with trips table with shapes table
trip_routes <- atl$trips %>% 
  full_join(atl$routes, by = "route_id")

trip_shape <- atlsf$shapes %>% 
  full_join(atl$trips %>% 
              select(shape_id, trip_id), 
            by = "shape_id")

# Merging the two into one and then taking only one row for each
# unique combination of route_id and shape_id.
route_trip_shape <- trip_shape %>% 
  select(-shape_id) %>% 
  full_join(trip_routes, by = c("trip_id")) %>% 
  group_by(shape_id, route_id) %>% 
  slice(1)

# Route type is not really intuitive - let's fix that
route.shape <- route_trip_shape %>% 
  mutate(route_type = case_when(
    route_type == "0" ~ 'Tram, Streetcar',
    route_type == "1" ~ 'Subway, Metro',
    route_type == "2" ~ 'Rail',
    route_type == "3" ~ 'Bus'
  ))

pal <- leaflet::colorFactor(c("red", "orange", "pink"), domain = route.shape$route_type)

route.shape %>% 
  leaflet::leaflet(data = .) %>% 
    leaflet::addProviderTiles(providers$CartoDB.DarkMatter) %>% 
    leaflet::addPolylines(color = ~pal(route_type), 
                 weight = 3,
                 opacity = 0.9,
                 popup = paste0("Route type: ", route.shape$route_type))

Step 2: Drawing buffers

We can use st_buffer() function from sf package to draw buffers from the transit routes. To be more precise, we would need to draw buffers around each stop, as people cannot just hop on bus or train between stations. For this class, however, we will simplify the process by drawing buffers around the line geometry given the fact that bus stations are so closely located to one another and their 400 meter buffer (i.e., commonly used distance for walking) would be almost the same as buffer around the line.

Drawing buffers is a simple process, but there is a caveat to it which is rooted in the recent change in sf package. Until Version 1.0, the sf package used ‘equirectangular projection’ using GEOS library for many operations that involve, for example, st_intersects, st_intersection, st_union, st_nearest_point, st_join, etc. What this means is that sf package carried out these operations assuming that the Earth is a flat surface. It did so even when the given data is in geographical coordinates.

From Version 1.0, sf package started using S2 spherical geometry for spatial operations, based on S2 geometry library written by Google. For our purposes, the difference introduced in Version 1.0. is unlikely to have huge impacts. However, we may see some glitches here and there. For example, in the code chunks below I transform route.shape from crs = 4326 to crs = 26967. If we use crs = 4326, the buffer polygons will have rough outlines with zig-zag patterns. You can try it yourself. For more information on what all of these means, you can read this post.

# If we do not convert CRS to a projected one, the buffer may generate a ragged boundaries. 
# https://r-spatial.github.io/sf/articles/sf7.html#buffers-1
# https://r-spatial.org/r/2020/06/17/s2.html#sf-10-goodbye-flat-earth-welcome-s2-spherical-geometry

# Buffering LINESTRING ---------------------------------------
sf::sf_use_s2(TRUE) # This is default to TRUE when we load SF package. So this code is not really needed. I added it here just to make it more explicit.
MARTA_buffer <- route.shape %>% 
  sf::st_transform(crs = 26967) %>% 
  sf::st_buffer(dist = units::set_units(400, "m")) 

# To union the buffer polygons by route_type
MARTA_buffer_group <- MARTA_buffer %>% 
  group_by(route_type) %>% 
  summarise()

# Buffering POINTS ---------------------------------------
# Just to show that lines & points overlap well
MARTA_stop_buffer <- atlsf$stops %>% 
  st_transform(crs = 26967) %>% 
  st_buffer(dist = units::set_units(400, "m")) %>% 
  st_union()

pal <- colorFactor(palette = c("red", "yellow", "blue"), domain = MARTA_buffer_group$route_type)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  
  addPolygons(data = MARTA_buffer_group %>% st_transform(crs = 4326), col = ~pal(route_type),
              popup = MARTA_buffer_group$route_type,
              group = "Buffer from Line") %>%
  
  addPolygons(data = MARTA_stop_buffer %>% st_transform(crs = 4326),
              weight = 0.5, 
              color = "white",
              group = "Buffer from Point") %>% 
  addLayersControl(
    overlayGroups = c("Buffer from Line", "Buffer from Point"),
    options = layersControlOptions(collapsed = FALSE)
  )

Step 3: Intersect the buffer with Census Tracts,

We will use tidycensus package again to download Census data at Tract-level.

And then, of course, we do some data cleaning. In the code, I use st_area() to calculate the area of the polygon. Note that I unclass() the output from st_area() because st_area() does not return a plain numeric value; the return is of ‘units’ class. This is to provide not only the area value in numeric but also the unit in which the area is measured in. While this is highly useful, having something other than plain numeric values may cause some complications later. So I drop the class ‘unit’ and convert it to plain numeric values.

# Cleaning Census data
acs2020c <- acs2020 %>% 
  select(GEOID,
         hhinc = hhincE,
         r_tot = r_totE,
         r_wh = r_whE,
         r_bl = r_blE,
         tot_hh = tot_hhE,
         own_novhc = own_novhcE,
         rent_novhc = rent_novhcE) %>% 
  mutate(pct_wh = r_wh / r_tot,
         pct_bl = r_bl / r_tot,
         pct_novhc = (own_novhc + rent_novhc)/tot_hh) %>% 
  # Calculate area of the Census Tract polygons
  mutate(area = st_area(.) %>% unclass) %>% 
  mutate(ln_pop_den = log((r_tot / (area/1000^2)) + 1)) %>% 
  filter(!is.na(hhinc), !is.na(r_tot), !is.na(own_novhc))

pal1 <- colorNumeric(palette = "YlOrRd", domain = acs2020c$hhinc)
pal2 <- colorFactor("Spectral", domain = MARTA_buffer_group$route_type)

# Mapping
leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(group = "ACS",
              data = acs2020c, 
              color = "grey", 
              fillColor = ~pal1(hhinc), 
              fillOpacity = 0.5,
              weight = 1, 
              popup = leafpop::popupTable(round(st_drop_geometry(acs2020c[,c("hhinc", "pct_wh", "pct_bl")]),2))) %>% 
  addPolygons(group = "MARTA",
              data = MARTA_buffer_group %>% 
                st_transform(crs = st_crs(acs2020c)), 
              color = ~pal2(route_type),
              weight = 1,
              opacity = 0.9) %>% 
  addLayersControl(
    overlayGroups = c("ACS", "MARTA"),
    options = layersControlOptions(collapsed = FALSE)
  )

## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'

## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'

Then we merge the Census data and the buffers from the GTFS. This requires st_intersection(), not st_join(), because we want to ‘cut’ the Census polygons by the shape of GTFS buffers so that we will be able to calculate the area of the intersection.

Notice that the output from st_intersection() can be different from the input. While st_join returns the same geometry as the input data (although the same geometry can be duplicated), st_intersection() directly edits the geometry (see the image below).

# Intersect buffer with tract 
buffer_intersect_tract <- acs2020c %>% 
  # Unify the CRS
  st_transform(crs = st_crs(MARTA_buffer_group)) %>% 
  # Intersection
  st_intersection(MARTA_buffer_group) %>% 
  # Extract bus routes only
  filter(route_type == "Bus")

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

# Visualize
pal_buf <- colorFactor(palette = "Spectral", domain = buffer_intersect_tract$route_type)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(data = acs2020c %>% st_transform(crs = 4326), fillOpacity = 0.2, color = "white", weight = 0.5, group = 'Census') %>% 
  addPolygons(data = MARTA_buffer_group %>% st_transform(crs = 4326), fillOpacity = 0.2, color = "yellow", weight  = 0.5, group = 'GTFS') %>% 
  addPolygons(data = buffer_intersect_tract %>% st_transform(crs = 4326), group = "Intersection",
              fillColor = ~pal_buf(route_type), fillOpacity = 0.9, weight = 1, opacity = 0.3, color = "grey") %>%     # leaflet takes 4326
  addLayersControl(overlayGroups = c("Census", "GTFS", "Intersection"), options = layersControlOptions(collapsed = FALSE))

Step 4. Calculate the proportion of each Census Tract that are within the buffers.

Finally, we need to calculate the proportion of each Census Tract that overlaps with the GTFS buffers. We will need to know the (1) total area of Census Tracts, which we calculated above and (2) the area of the intersected polygons. We will then divide the intersected polygons’ area by the total area and put that into a column named pct_served.

# Area of intersected area
buffer_intersect_tract2 <- buffer_intersect_tract %>% 
  mutate(subarea = unclass(st_area(.)),
         pct_served = subarea/area) %>% 
  # Back to WGS84
  st_transform(crs = 4326)

# Joining the intersected information back to the original Census Data.
# This is needed because we may have lost a few polygons during intersection.
buffer_by_tract <- acs2020c %>% 
  # In left_join(A, B), B must be non-sf object.
  left_join(buffer_intersect_tract2 %>% 
              select(GEOID, pct_served) %>% 
              st_set_geometry(NULL),
            by = "GEOID") %>% 
  # There are many NAs in the pct_served column because there are many Census Tracts
  # that do not overlap with the GTFS buffer. They would have NA in pct_served column.
  # This code converts the NAs into 0.
  mutate(pct_served = case_when(is.na(pct_served) ~ 0,
                                 TRUE ~ pct_served))

# Map it!
pal_serv <- colorNumeric(palette = "Spectral", domain = buffer_by_tract$pct_served, reverse = TRUE)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(data = buffer_by_tract, 
              fillColor = ~pal_serv(pct_served), fillOpacity = 0.9, 
              color = "white", opacity = 0.2, 
              weight = 1) %>% 
  addLegend("bottomright",
            pal = pal_serv,
            values = buffer_by_tract$pct_served,
            title = "% Area within 400m from transit line") 

## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'

Merge frequency by stop information with ACS data

We use spatial join to merge GTFS service frequency data to Census.

stop_tract <- acs2020c %>% 
  st_transform(crs = 4326) %>% 
  st_join(stop_freq_sf) %>% 
  group_by(GEOID) %>% 
  summarise(avg_n_departures = mean(n_departures),
            sum_n_departures = sum(n_departures),
            hhinc = mean(hhinc),
            r_tot = mean(r_tot),
            pct_wh = mean(pct_wh),
            pct_bl = mean(pct_bl),
            pct_novhc = mean(pct_novhc),
            r_tot = mean(r_tot),
            ln_pop_den = mean(ln_pop_den))

knitr::kable(head(stop_tract))

sum_n_depart <- colorQuantile(palette = "Reds", domain = stop_tract$sum_n_departures)

sum_depart_map <- leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(data = stop_tract, 
              fillColor = sum_n_depart(stop_tract$sum_n_departures),
              fillOpacity = 0.8,
              color = 'grey',
              weight = 1, 
              popup = paste0("ID: ", stop_tract$GEOID, ", Value: ", stop_tract$sum_n_departures)) %>% 
  addLegend(position = "bottomright",
            values = stop_tract$sum_n_departures,
            pal = sum_n_depart)

sum_depart_map

message(str_c("there are ", sum(is.na(stop_tract$sum_n_departures)), " NAs in n_departures"))

## there are 84 NAs in n_departures

# Is there any correlation between SES and service area?
stop_tract_plot <- stop_tract %>% 
  mutate(hhinc = log(hhinc),
         pct_novhc = log(pct_novhc + 0.05),
         n_departure_log = log(sum_n_departures + 1)) %>% 
  pivot_longer(cols = c('hhinc', 'pct_wh', 'pct_bl', 'pct_novhc', "r_tot", "ln_pop_den"), 
               names_to = "variable", values_to = "value") %>%
  mutate(variable = factor(variable, labels = c('Household Income', "Population Density", "% Black", 
                                      "% No car", "% White", "Total Population")))

departure_count <- stop_tract_plot %>% 
  ggplot(aes(x = n_departure_log, y = value)) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "lm") +
  facet_wrap(~variable, scales = "free_y") +
  labs(x = "Count of departures between 7AM and 10AM", title = "Departure Count VS. Socio-demographics") +
  theme_bw() 

plotly::ggplotly(departure_count)

## `geom_smooth()` using formula 'y ~ x'

## Warning: Removed 504 rows containing non-finite values (stat_smooth).

Calculate travel times from one station to all other stations (Optional)

The tidytransit package offers a very convenient function that calculates the shortest travel times from a stop to all other stops. Try ?tidytransit::travel_times to see what this function does and what arguments it takes. Notice that only two out of nine arguments are required arguments, and the rest are optional. The required arguments are ‘filtered_stop_times’ and ‘stop_name.’ If you do not provide the optional arguments, the function will use the pre-populated values for the calculation.

But if you run the function right away, there will be an error. The error is caused by the fact that there are stations that share names but are far from each other. There may be more than one stations in different parts of the city that share their names, which sometimes actually happen in real world. If duplicate names in our dataset is caused by these legitimate cases, the duplicate names should not be considered as errors. However, if in reality there are no such stations, the duplicate names would be actual flaws in your data.

The travel_times() function calculates travel time based on station names. This function is designed to throw an error if there are stops with the same name regardless of whether the duplicate names are legitimate reflection of the real world or flaws in the data. One exceptional case when duplicated names don’t cause error is when the stations with duplicate names are within 300 meters from each other – in this case, the function considers the two stations may actually be the same station and merges them.

For all other cases (duplicate names that are 200-meters apart), we need to modify the duplicate names so that the names are actually different. We do that through the following steps:

1. For each unique station names in stops table, calculate how many duplicates there are and, if there are duplicates, calculate distances to the duplicate stops. stop_group_distances() function can do this.
2. Get names of the stops that have duplicates far from each other. If stops with duplicate names are close (e.g., within 200 meters), they may actually be same stations.
3. For all stops with duplicates > 200-meter distances, change their station names by appending (1), (2), … (n) where n is the number of duplicates. For example, if there are two “MIDTOWN STATION”s, they will each be “MIDTOWN STATION (1)” and “MIDTOWN STATION (2)”.
4. One of the requirements for travel_times() function is that atl object contains a transfer table, which is absent when freshly imported from the hard drive. This is NOT always the case – the transfer table is an optional table in GTFS definition. Some entities include transfer table by default, some don’t. MARTA is among those who don’t include transfer table. We can use gtfs_transfer_table() function from gtfsrouter package to create this table and include it to the atl object.
5. Finally, one of the two required arguments, ‘filtered_stop_times’, takes an output from a function called tidytransit::filter_stop_times(). This function takes two required arguments: (1) our atl object and (2) a date for which we’d like calculate travel time. You can also choose to specify (3) minimum departure time and (4) maximum arrival time.

Finally, you are ready to calculate travel times!

# Step 1:
stop_dist <- stop_group_distances(atlsf$stops, by='stop_name') %>% 

# Step 2:  
  filter(dist_max > 200)

# Step 3
atl$stops <- atl$stops %>% 
  group_by(stop_name) %>% 
  mutate(stop_name = case_when(stop_name %in% stop_dist$stop_name ~ paste0(stop_name, " (", seq(1,n()), ")"),
                               TRUE ~ stop_name)) 

# Step 4
atl$transfers <- gtfsrouter::gtfs_transfer_table(atl, 
                                     d_limit = 200, 
                                     min_transfer_time = 120)
# Step 5
am_stop_time <- filter_stop_times(gtfs_obj = atl, 
                                  extract_date = "2021-08-14",
                                  min_departure_time = 3600*7, # input unit is in second. So 3500*7 is 7AM
                                  max_arrival_time = 3600*10)  # similarly, this is 10AM 

# travel_times
trvt <- travel_times(filtered_stop_times = am_stop_time,
                     stop_name = "MIDTOWN STATION",
                     return_coords = TRUE)

# ..and visualize the output
trvt_pal <- colorQuantile(palette = "Reds", domain = trvt$travel_time)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addCircles(data = trvt %>% st_as_sf(coords = c("to_stop_lon", "to_stop_lat"), crs = 4326), # converting trvt into sf object on the fly.
             fillColor = ~trvt_pal(travel_time), # Define color
             stroke = F, # Turn off the outer border lines of circles
             radius = 300, # Size of the circle
             fillOpacity = 0.7, # Transparency
             popup = paste0("Travel Time from <br> <strong> MIDTOWN: ", round(trvt$travel_time/60, 2), " minutes <strong>") %>% # Defines what's displayed on popup
               lapply(htmltools::HTML))

On average, you can see that it takes 50.1 minutes from the Midtown Station to other stations. On the map, you can see that it takes longer time to get to stops that are further away from the Midtown Station.

# Average time to get to other stations
print(str_c("On average, it takes ", 
            round(mean(trvt$travel_time)/60,1), 
            " minutes to travel from the Midtown Station to all other stations."))

## [1] "On average, it takes 50.1 minutes to travel from the Midtown Station to all other stations."