Major-Assignment 1

Task description

Park-and-Ride Simulation

In this assignment, you will simulate a journey that begins at a Census Tract centroid, drives to the nearest MARTA station, and then take the MARTA to the Midtown station. You will then compare how travel time varies by different Census Tracts. The steps and data required for this analysis are outlined below.

Step 1. Download the required GTFS data. Convert it to sf format, extract MARTA rail stations, and clean the stop names to remove duplicates. Also, extract the destination station.

Step 2. Download the required Census data. Convert Census Tract polygons into centroids and create a subset for analysis.

Step 3. Download the required OSM data. Convert it into an sfnetwork object and clean the network.

Step 4. Simulate a park-and-ride trip (from a test origin → closest station → Midtown station).

Step 5. Turn the simulation process from Step 4 into a reusable function.

Step 6. Apply the function from Step 5 iteratively to all home locations.

Step 7.Finally, create maps and plots to analyze whether there are any disparities in transit accessibility when commuting to Midtown.

Step 0. Packages

Importing the necessary packages is part of this assignment. Add any required packages to the code chunk below as you progress through the tasks.

library(tidyverse)
library(gtfsrouter) 
library(sf)
library(osmdata)
library(sfnetworks)
library(tidycensus)
library(ggplot2) 
library(purrr)
library(hms)
library(tidytransit)
library(tmap)
library(here)
library(units)
library(leaflet)
library(plotly)
library(tidygraph)
library(leafsync)
library(leafpop)
library(magrittr)
library(dbscan)
library(tigris)
library(nominatimlite)

Step 1. GTFS data

# TASK ////////////////////////////////////////////////////////////////////////
# Download MARTA GTFS data and assign it to `gtfs` object
gtfs <- read_gtfs("https://www.itsmarta.com/google_transit_feed/google_transit.zip")
# //TASK //////////////////////////////////////////////////////////////////////



# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# Edit stop_name to append serial numbers (1, 2, etc.) to remove duplicate names
stop_dist <- stop_group_distances(gtfs$stops, by='stop_name') %>%
  filter(dist_max > 200)

gtfs$stops <- gtfs$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))

# Create a transfer table
gtfs <- gtfsrouter::gtfs_transfer_table(gtfs, 
                                        d_limit = 200, 
                                        min_transfer_time = 120)

# NOTE: Converting to sf format uses stop_lat and stop_lon columns contained in gtfs$stops.
#       In the conversion process, stop_lat and stop_lon are converted into a geometry column, and
#       the output sf object do not have the lat lon column anymore.
#       But many other functions in tidytransit look for stop_lat and stop_lon.
#       So I re-create them using mutate().
gtfs <- gtfs %>% gtfs_as_sf(crs = 4326)

gtfs$stops <- gtfs$stops %>% 
  ungroup() %>% 
  mutate(stop_lat = st_coordinates(.)[,2],
         stop_lon = st_coordinates(.)[,1]) 

# Get stop_id for rails and buses
rail_stops <- gtfs$routes %>% 
  filter(route_type %in% c(1)) %>% 
  inner_join(gtfs$trips, by = "route_id") %>% 
  inner_join(gtfs$stop_times, by = "trip_id") %>% 
  inner_join(gtfs$stops, by = "stop_id") %>% 
  group_by(stop_id) %>% 
  slice(1) %>% 
  pull(stop_id)

# Extract MARTA rail stations
station <- gtfs$stops %>% filter(stop_id %in% rail_stops)

# Extract Midtown Station
midtown <- gtfs$stops %>% filter(stop_id == "134")

# Create a bounding box to which we limit our analysis
bbox <- st_bbox(c(xmin = -84.45241, ymin = 33.72109, xmax = -84.35009, ymax = 33.80101), 
                 crs = st_crs(4326)) %>% 
  st_as_sfc()

# =========== NO MODIFY ZONE ENDS HERE ========================================

Step 2. Census data

# TASK ////////////////////////////////////////////////////////////////////////
# Specify Census API key whichever you prefer using census_api_key() function
census_api_key(Sys.getenv("CENSUS_API"))
# //TASK //////////////////////////////////////////////////////////////////////

# TASK ////////////////////////////////////////////////////////////////////////
# Using get_acs() function, download Census Tract level data for 2022 for Fulton, DeKalb, and Clayton in GA.
# and assign it to `census` object.
# Make sure you set geometry = TRUE.

# Required data from the Census ACS:
#  1) Median Household Income (name the column `hhinc`)
#  2) Minority Population (%) (name the column `pct_minority`)
# Note: You may need to download two or more Census ACS variables to calculate minority population (%). "Minority" here can refer to either racial minorities or racial+ethnic minorities -- it's your choice.
acs_vars <- c(
  hhinc = "B19013_001",
  total_pop = "B03002_001E",
  white_non_hisp = "B03002_003E"
)

census <- get_acs(
  geography = "tract",
  variables = acs_vars,
  state = "GA",
  county = c("Fulton", "DeKalb", "Clayton"),
  year = 2022,
  geometry = TRUE,
  output = "wide"
) %>%
  mutate(
    pct_minority = (total_pop - white_non_hisp) / total_pop * 100
  ) %>%
  select(
    GEOID, NAME,
    hhinc = hhincE,
    pct_minority,
    geometry
  )

# //TASK //////////////////////////////////////////////////////////////////////


# TASK ////////////////////////////////////////////////////////////////////////
# Convert the CRS of `census` to GCS

census <- census %>% 
  st_transform(crs = 4326)
# //TASK //////////////////////////////////////////////////////////////////////


# TASK ////////////////////////////////////////////////////////////////////////
# get centroids of `census` polygons, 
# extract centroids that fall inside the `bbox`, 
# and assign it to `home` object.

home <- census %>%
  st_centroid() %>%
  st_filter(bbox)

# //TASK //////////////////////////////////////////////////////////////////////

Step 3. OSM data

# TASK ////////////////////////////////////////////////////////////////////////
# 1. Get OSM road data for all road types for the `bbox` area.
# 3. Convert the OSM data into an sf object
# 4. Convert osmdata polygons into lines

osm_road <- opq(bbox=bbox) %>%
  add_osm_feature(key = 'highway', 
                  value = c("motorway", "motorway_link",
                            "trunk", "trunk_link", 
                            "primary", "primary_link",
                            "secondary", "secondary_link",
                            "tertiary", "residential")) %>%
  osmdata_sf() %>% 
  osm_poly2line()

names(osm_road)

## [1] "bbox"              "overpass_call"     "meta"             
## [4] "osm_points"        "osm_lines"         "osm_polygons"     
## [7] "osm_multilines"    "osm_multipolygons"

# //TASK //////////////////////////////////////////////////////////////////////


# TASK ////////////////////////////////////////////////////////////////////////
# 1. Convert osm_road$osm_line into sfnetwork
# 2. Activate edges
# 3. Clean the network using the methods we learned
# 4. Assign the cleaned network to an object named 'osm'

osm <- osm_road$osm_line[bbox,] %>% select(osm_id, highway)
# converting the line component of OSM data into a network
net <- sfnetworks::as_sfnetwork(osm, directed = FALSE)
print(net)

## # A sfnetwork with 6215 nodes and 5589 edges
## #
## # CRS:  EPSG:4326 
## #
## # An undirected multigraph with 1196 components with spatially explicit edges
## #
## # Node data: 6,215 × 1 (active)
##               geometry
##            <POINT [°]>
## 1 (-84.34923 33.74581)
## 2  (-84.35038 33.7454)
## 3 (-84.35015 33.74516)
## 4   (-84.34924 33.745)
## 5 (-84.37242 33.74325)
## 6 (-84.36738 33.74324)
## # ℹ 6,209 more rows
## #
## # Edge data: 5,589 × 5
##    from    to osm_id  highway                                           geometry
##   <int> <int> <chr>   <chr>                                     <LINESTRING [°]>
## 1     1     2 9164335 motorway_link (-84.34923 33.74581, -84.3491 33.74624, -84…
## 2     3     4 9165104 motorway_link (-84.35015 33.74516, -84.34948 33.74517, -8…
## 3     5     6 9186247 motorway      (-84.37242 33.74325, -84.37048 33.74326, -8…
## # ℹ 5,586 more rows

print(paste0('Which one is active?: ', sfnetworks::active(net)))

## [1] "Which one is active?: nodes"

# simplify our network
simple_net <- net %>%
  activate("edges") %>%
  filter(!edge_is_multiple()) %>% # remove multiple edges
  filter(!edge_is_loop()) # remove loops

# print out the differences in the edge count.
message(str_c("Before simplification, there were ", net %>% st_as_sf("edges") %>% nrow(), " edges. \n",
            "After simplification, there are ", simple_net %>% st_as_sf("edges") %>% nrow(), " edges."))

## Before simplification, there were 5589 edges. 
## After simplification, there are 5545 edges.

# subdivide edges
subdiv_net <- convert(simple_net, sfnetworks::to_spatial_subdivision)

## Warning: to_spatial_subdivision assumes attributes are constant over geometries

subdiv_net <- subdiv_net %>% 
  activate("nodes") %>% 
  mutate(custom_id = seq(1, subdiv_net %>% st_as_sf("nodes") %>% nrow()),
         is.new = case_when(is.na(.tidygraph_node_index) ~ "new nodes",
                            TRUE ~ "existing nodes"),
         is.new = factor(is.new)) 

# using spatial morpher
smoothed_net <- convert(subdiv_net, sfnetworks::to_spatial_smooth)

# remove
removed <- setdiff(subdiv_net %>% st_as_sf('nodes') %>% pull(custom_id), 
                   smoothed_net %>% st_as_sf('nodes') %>% pull(custom_id))

osm <- smoothed_net

# //TASK //////////////////////////////////////////////////////////////////////


# TASK ////////////////////////////////////////////////////////////////////////
# Add a new column named 'length' to the edges part of the object `osm`.
osm %<>%
  st_transform(4326) %>% 
  activate("edges") %>% 
  mutate(length = edge_length())
# //TASK //////////////////////////////////////////////////////////////////////

Step 4. Simulate a park-and-ride trip

# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# Extract the first row from `home` object and store it `home_1`
home_1 <- home[1,]
# =========== NO MODIFY ZONE ENDS HERE ========================================


# TASK ////////////////////////////////////////////////////////////////////////
# Find the shortest paths from `home_1` to all other stations
paths <- st_network_paths(osm, 
                  from = home_1, 
                  to = station, 
                  type = "shortest")
# //TASK //////////////////////////////////////////////////////////////////////

  
# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# Get shortest network distances from `home_1` to all other stations.
dist_all <- map_dbl(1:nrow(paths), function(x){
  osm %>%
    activate("nodes") %>% 
    slice(paths$node_paths[[x]]) %>% 
    st_as_sf("edges") %>% 
    pull(length) %>% 
    sum()
}) %>% unlist()

# Replace zeros with a large value.
if (any(dist_all == 0)){
  dist_all[dist_all == 0] <- max(dist_all)
}

# TASK ////////////////////////////////////////////////////////////////////////
# Find the closest station and assign the value to `closest_station`

closest_station_index <- which.min(dist_all)
closest_station <- station[closest_station_index, ]
# //TASK //////////////////////////////////////////////////////////////////////


# TASK ////////////////////////////////////////////////////////////////////////
# Get the distance to the closest station.

closest_dist <- min(dist_all)
# //TASK //////////////////////////////////////////////////////////////////////


# TASK ////////////////////////////////////////////////////////////////////////
# Calculate travel time based on the `closest_dist` assuming we drive at 20 miles/hour speed.
# Assign the value to `trvt_osm_m` object.

car_speed <- set_units(20, mile/h)
trvt_osm_m <- closest_dist/set_units(car_speed, m/min) %>%
  as.vector(.)
# //TASK //////////////////////////////////////////////////////////////////////


# TASK ////////////////////////////////////////////////////////////////////////
# Create a subset of stop_times data table for date = 2025-11-10, 
# minimum departure time of 7AM, maximum departure time of 10AM.
# Assign the output to `am_stop_time` object
am_stop_time <- filter_stop_times(gtfs_obj = gtfs, 
                                  extract_date = "2025-11-10", # can't go back in time
                                  min_departure_time = 3600*7,
                                  max_arrival_time = 3600*10)
# //TASK //////////////////////////////////////////////////////////////////////


# TASK ////////////////////////////////////////////////////////////////////////
# 1. Calculate travel times from the `closest_station` to all other stations 
#    during time specified in `am_stop_time`. Allow ONE transfer.
# 2. Filter the row where the destination is Midtown station. 
# 3. Assign it to `trvt` object.

trvt <- travel_times(filtered_stop_times = am_stop_time,
                     stop_name = closest_station$stop_name,
                     time_range = 3600,
                     arrival = FALSE,
                     max_transfers = 1,
                     return_coords = TRUE) %>%
  filter(to_stop_name == midtown$stop_name[1])
# //TASK //////////////////////////////////////////////////////////////////////


# =========== NO MODIFICATION ZONE STARTS HERE ===============================
# Divide the calculated travel time by 60 to convert the unit from seconds to minutes.
trvt_gtfs_m <- trvt$travel_time/60

# Add the travel time from home to the nearest station and
# the travel time from the nearest station to Midtown station
total_trvt <- trvt_osm_m + trvt_gtfs_m
# =========== NO MODIFY ZONE ENDS HERE ========================================

Step 5. Convert Step 4 into a function

# Function definition (do not modify other parts of the code in this code chunk except for those inside the TASK section)

get_trvt <- function(home, osm, station, midtown){
  
  # TASK ////////////////////////////////////////
  # If the code in Step 4 runs fine,
  # Replace where it says **YOUR CODE HERE..** below with 
  # the ENTIRETY of the previous code chunk (i.e., Step 4)

  paths <- st_network_paths(osm, 
                  from = home, 
                  to = station, 
                  type = "shortest")
  dist_all <- map_dbl(1:nrow(paths), function(x){
    osm %>%
      activate("nodes") %>% 
      slice(paths$node_paths[[x]]) %>% 
      st_as_sf("edges") %>% 
      pull(length) %>% 
      sum()
  }) %>% unlist()

  if (any(dist_all == 0)){
    dist_all[dist_all == 0] <- max(dist_all)
  }

  closest_station_index <- which.min(dist_all)
  closest_station <- station[closest_station_index, ]
  
  closest_dist <- min(dist_all)
  
  car_speed <- set_units(20, mile/h)
  trvt_osm_m <- closest_dist/set_units(car_speed, m/min) %>%
    as.vector(.)
  
  am_stop_time <- filter_stop_times(gtfs_obj = gtfs, 
                                    extract_date = "2025-11-10", # can't go back in time
                                    min_departure_time = 3600*7,
                                    max_arrival_time = 3600*10)
  
  trvt <- travel_times(filtered_stop_times = am_stop_time,
                       stop_name = closest_station$stop_name,
                       time_range = 3600,
                       arrival = FALSE,
                       max_transfers = 1,
                       return_coords = TRUE) %>%
    filter(to_stop_name == midtown$stop_name[1])
  
  trvt_gtfs_m <- trvt$travel_time/60
  total_trvt <- trvt_osm_m + trvt_gtfs_m
  
  # //TASK //////////////////////////////////////

  # =========== NO MODIFICATION ZONE STARTS HERE ===============================
  if (length(total_trvt) == 0) {total_trvt = 0}

  return(total_trvt)
  # =========== NO MODIFY ZONE ENDS HERE ========================================
}

Step 6. Apply the function to all home locations.

# Prepare an empty vector
total_trvt <- vector("numeric", nrow(home))

# Apply the function to all Census Tracts
# Fill `total_trvt` object with the calculated time
for (i in 1:nrow(home)){
  total_trvt[i] <- get_trvt(home[i,], osm, station, midtown)
}

# cbind the calculated travel time to `home`
home <- home %>% 
  cbind(trvt = total_trvt)

Step 7. Create maps and plots

Create two maps and two plots by following the instructions in the code chunk below. Write a brief (maximum 200 words) description summarizing your observations from the maps and plots

• Response: The spatial analysis reveals distinct socioeconomic patterns in park-and-ride accessibility. Higher-income tracts, concentrated in northern and suburban areas, generally exhibit longer total travel times to Midtown. However, these longer times do not necessarily indicate poorer accessibility. Rather, they reflect residential and lifestyle choices: higher-income households often reside in suburban environments characterized by larger housing, lower density, and perceived advantages in safety and environmental quality. The extended travel time arises primarily from longer driving distances to reach rail stations rather than inefficiency in the transit system itself. Conversely, lower-income tracts, typically situated within the urban core and near MARTA rail corridors, display shorter combined travel times. Their proximity reflects structural dependence on public transit and constrained residential choice rather than deliberate preference for centrality. The weak association between minority population share and travel time further suggests that spatial variation is shaped more by socioeconomic and locational preferences than by direct demographic disadvantage. From my perspective, the analysis underscores that longer travel times among higher-income residents are not necessarily indicative of reduced accessibility, but rather a reflection of voluntary suburbanization and lifestyle differentiation—emphasizing that spatial inequality does not always equate to transit inequity.

# TASK ////////////////////////////////////////
# Create an interactive map displaying `census` (polygons) and `home` (points), effectively visualizing household income and travel time to Midtown Station, respectively.
tmap_mode('view')

## ℹ tmap mode set to "view".

census_map <- census %>% st_filter(bbox)
home_map <- home %>% st_filter(bbox)

pal_income   <- colorNumeric(palette = "YlGnBu", domain = census$hhinc)
pal_minority <- colorNumeric(palette = "Blues",   domain = census$pct_minority)
pal_travel   <- colorNumeric(palette = "Reds",   domain = home$trvt)

m_income <- leaflet() %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  # census polygons
  addPolygons(
    data = census_map,
    fillColor = ~pal_income(hhinc),
    fillOpacity = 0.7,
    weight = 0.2,
    color = "#333333",
    popup = ~paste0("<b>", NAME, "</b><br>",
                    "Median HH Income: $", scales::comma(hhinc))
  ) %>%
  # home points, travel time
  addCircles(
    data = home_map,
    fillColor = ~pal_travel(trvt),
    fillOpacity = 0.9,
    stroke = FALSE,
    radius = 70,
    popup = ~paste0("<b>", NAME, "</b><br>",
                    "Travel time to Midtown: ", round(trvt, 1), " min")
  ) %>%
  addCircles(data = midtown, color = "yellow", radius = 200, weight = 3) %>%
  addLegend(
    position = "bottomright",
    pal = pal_income,
    values = census$hhinc,
    title = "Median HH Income ($)"
  ) %>%
  addLegend(
    position = "bottomleft",
    pal = pal_travel,
    values = home$trvt,
    title = "Travel Time (min)"
  )

m_income

# Create an interactive map displaying `census` (polygons) and `home` (points) effectively visualizing the percentage of minority population and travel time to Midtown Station, respectively.

m_minority <- leaflet() %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(
    data = census_map,
    fillColor = ~pal_minority(pct_minority),
    fillOpacity = 0.7,
    weight = 0.2,
    color = "#333333",
    popup = ~paste0("<b>", NAME, "</b><br>",
                    "Minority Population: ", round(pct_minority, 1), "%")
  ) %>%
  addCircles(
    data = home_map,
    fillColor = ~pal_travel(trvt),
    fillOpacity = 0.9,
    stroke = FALSE,
    radius = 70,
    popup = ~paste0("<b>", NAME, "</b><br>",
                    "Travel time to Midtown: ", round(trvt, 2), " min")
  ) %>%
  addCircles(data = midtown, color = "yellow", radius = 200, weight = 3) %>%
  addLegend(
    position = "bottomright",
    pal = pal_minority,
    values = census$pct_minority,
    title = "Minority Population (%)"
  ) %>%
  addLegend(
    position = "bottomleft",
    pal = pal_travel,
    values = home$trvt,
    title = "Travel Time (min)"
  )

m_minority

# Create a scatter plot with a trend line showing the relationship between household income (x-axis) and travel time to Midtown Station (y-axis).

plot_income <- census %>%
  st_drop_geometry() %>%
  left_join(home %>% st_drop_geometry() %>% select(GEOID, trvt), by = "GEOID") %>%
  ggplot(aes(x = hhinc, y = trvt)) +
  geom_point(color = "black", alpha = 0.6, size = 1.8) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 0.8) +
  theme_minimal(base_family = "Helvetica") +
  theme(
    plot.background  = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA),
    axis.text        = element_text(color = "black"),
    axis.title       = element_text(color = "black"),
    plot.title       = element_text(color = "black", face = "bold", size = 13),
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "Income vs. Travel Time to Midtown",
    x = "Median Household Income ($)",
    y = "Travel Time (min)"
  )

plot_income

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

# Create a scatter plot with a trend line showing the relationship between the percentage of minority population (x-axis) and travel time to Midtown Station (y-axis).

plot_minority <- census %>%
  st_drop_geometry() %>%
  left_join(home %>% st_drop_geometry() %>% select(GEOID, trvt), by = "GEOID") %>%
  ggplot(aes(x = pct_minority, y = trvt)) +
  geom_point(color = "black", alpha = 0.6, size = 1.8) +
  geom_smooth(method = "lm", se = FALSE, color = "blue", linewidth = 0.8) +
  theme_minimal(base_family = "Helvetica") +
  theme(
    plot.background  = element_rect(fill = "white", color = NA),
    panel.background = element_rect(fill = "white", color = NA),
    axis.text        = element_text(color = "black"),
    axis.title       = element_text(color = "black"),
    plot.title       = element_text(color = "black", face = "bold", size = 13),
    panel.grid.minor = element_blank()
  ) +
  labs(
    title = "Minority Share vs. Travel Time to Midtown",
    x = "Minority Population (%)",
    y = "Travel Time (min)"
  )

plot_minority

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

# //TASK //////////////////////////////////////