Major Assignment_1_

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

knitr::opts_chunk$set(message = FALSE, warning = FALSE)

# Load the packages used in your analysis.
library(tidyverse)   # data wrangling (dplyr, tidyr) and ggplot2
library(sf)          # spatial data
library(tidytransit) # read GTFS, filter stop times, travel_times()
library(gtfsrouter)  # gtfs_transfer_table(), as in the GTFS lab
library(tidycensus)  # ACS estimates
library(osmdata)     # download OSM roads
library(sfnetworks)  # road network as a spatial graph
library(tidygraph)   # graph helpers (edge_is_multiple, group_components)
library(leaflet)     # interactive maps
library(leafsync)    # synchronized side-by-side maps
library(scales)      # axis labels

# The Census API key is set once in the console with
# tidycensus::census_api_key("YOUR KEY", install = TRUE); it is not stored here.

How to Use This Starter File

This file provides the structure for your analysis. Complete the blank code chunks, provide explanations where requested, and add code comments where necessary. Keep the existing code chunks and section headings.

Your implementation may differ from the lab examples, but the checkpoint objects named below must exist because later sections and validation checks use them. Explain important choices with concise code comments and short narrative text.

Work through the document in order and knit as you complete sections. Unfinished code and validation checks may prevent the document from knitting. While working, you may temporarily comment out those lines (Ctrl/Cmd + Shift + C). Uncomment them using the same shortcut and confirm that the entire document knits before submitting. Do not include a Census API key or any other private credential in this file.

Use the methods and coding patterns covered in class and the course labs. Substantial unexplained departures may result in a grade penalty, and you may be asked to explain any part of your code.

Assignment Overview

MARTA wants to understand how long it takes residents in different parts of Atlanta to reach North Avenue Station using a park-and-ride trip. Each simulated journey has two stages:

1. Travel from a Census tract centroid to the reachable MARTA rail station offering free daily parking with the shortest estimated driving time along the selected OSM road network.
2. Travel from that station to North Avenue Station using MARTA service during the assigned morning period.

You will combine OSM, GTFS, and Census data; build a reusable travel-time function; apply it to multiple origins; and evaluate whether estimated travel times differ across neighborhoods.

Simulation Assumptions

Use the following assumptions throughout your analysis:

• Use the supplied road-class speeds to calculate edge travel times. These are assignment assumptions, not observed average speeds or posted speed limits.
• Treat roads as undirected. Ignore one-way restrictions, congestion, traffic signals, and turn delays.
• Select the eligible parking station with the shortest estimated driving time. Do not select by road distance or by combined driving and transit time.
• Routing uses the nearest network nodes to the tract centroids and stations. Ignore travel between these points and their network nodes.
• Assume a parking space is available. Ignore parking-search time and walking time from the parking lot to the station platform.
• Do not add an initial waiting-time allowance before boarding the first train.
• Use the scheduled journey duration returned by travel_times(), allowing at most one transfer. Retain any transfer duration and waiting between trains already included in that result; do not subtract them or add them again.
• For the transfers table, use a maximum connection distance of 200 meters and a minimum transfer duration of 120 seconds, as in the GTFS lab.
• You do not need to restrict GTFS routing to rail-only itineraries. Both endpoints are MARTA rail stations, so rail is the expected connection. Use the journey returned by travel_times(), even if it includes a bus segment.
• Add driving time and transit journey time without synchronizing the train departure with arrival by car.

The result is a simplified driving-plus-transit travel-time estimate, rather than a fully timed door-to-door journey.

Analysis Settings

The bounding box below defines the study area used throughout the assignment.

study_bbox <- sf::st_bbox(
  c(
    xmin = -84.4753,
    ymin = 33.7151,
    xmax = -84.3091,
    ymax = 33.8283
  ),
  crs = sf::st_crs(4326)
)

Provided Parking Data

For this assignment, a rail station is eligible as a park-and-ride station only if MARTA lists it as offering free daily parking. The vector below uses the station names found in the current GTFS feed. It represents whether a station offers free daily parking, not whether a parking space is available in real time. Apply the parking-availability assumption above.

Sources: MARTA Parking Availability and Ashby Station

free_parking_station_names <- c(
  "ASHBY STATION",
  "AVONDALE STATION",
  "BROOKHAVEN-OGLETHORPE STATION",
  "CHAMBLEE STATION",
  "COLLEGE PARK STATION",
  "DORAVILLE STATION",
  "DUNWOODY STATION",
  "EAST LAKE STATION",
  "EAST POINT STATION",
  "EDGEWOOD-CANDLER PARK STATION",
  "HAMILTON E HOLMES STATION",
  "INDIAN CREEK STATION",
  "INMAN PARK-REYNOLDSTOWN STATION",
  "KENSINGTON STATION",
  "LAKEWOOD-FT MCPHERSON STATION",
  "LENOX STATION",
  "LINDBERGH CENTER STATION",
  "MEDICAL CENTER STATION",
  "NORTH SPRINGS STATION",
  "OAKLAND CITY STATION",
  "SANDY SPRINGS STATION",
  "WEST END STATION",
  "WEST LAKE STATION"
)

1. Prepare GTFS Data

Download the MARTA GTFS feed and prepare the information needed for transit routing between rail stations.

Your workflow should:

• identify rail service using the GTFS route type;
• connect routes, trips, stop times, and stops using the appropriate keys;
• represent each physical rail station once using parent_station when available;
• retain stations within the study area using study_bbox;
• restrict potential park-and-ride stations to those listed in free_parking_station_names;
• identify North Avenue Station as the transit destination; and
• create any transfer or filtered stop-time tables needed for routing during the assigned morning period.

Prepare transit routing for November 18, 2026, between 7 AM and 10 AM. For the selected parking station, use the minimum scheduled journey duration to North Avenue among journeys departing during this window and arriving by 10 AM. Set time_range to cover the full three-hour departure window. The relevant GTFS routing functions express times as seconds after midnight.

The study-area restriction applies to station points used as candidate park-and-ride locations. Keep the full GTFS service network for transit routing, including trips and transfers that may extend outside the study area.

Create these checkpoint objects:

• gtfs: the prepared GTFS feed;
• rail_stations: one spatial point per physical rail station in the study area;
• park_ride_stations: the subset of rail_stations that offers free daily parking;
• destination_stop: the North Avenue Station record; and
• morning_stop_times: stop times prepared for the assigned date and time window.

Parent stations and their platforms can have different names. For example, North Avenue’s parent record uses "NORTH AVE STATION", while its platforms use "NORTH AVENUE STATION". Use the supplied code below after downloading gtfs and before creating spatial stop objects or routing tables. It gives platforms their parent station’s name without changing stop IDs or parent-station links.

# Download and prepare the MARTA GTFS feed.
gtfs <- read_gtfs("https://www.itsmarta.com/google_transit_feed/google_transit.zip")

analysis_date <- as.Date("2026-11-18")
am_start <- 7 * 3600   # 7 AM in seconds after midnight
am_end   <- 10 * 3600  # 10 AM in seconds after midnight

# Provided: make platform names consistent with their parent station names.
parent_names <- gtfs$stops %>%
  select(parent_station = stop_id, parent_name = stop_name)

gtfs$stops <- gtfs$stops %>%
  left_join(parent_names, by = "parent_station") %>%
  mutate(stop_name = coalesce(parent_name, stop_name)) %>%
  select(-parent_name)

# Identify physical rail stations, eligible parking stations, and North Avenue Station.
# Rail service: route_type 1 (subway/metro) is MARTA's heavy rail.
rail_routes <- gtfs$routes %>%
  filter(route_type == 1)

# routes -> trips (route_id) -> stop_times (trip_id) -> stops (stop_id)
rail_trips <- gtfs$trips %>%
  semi_join(rail_routes, by = "route_id")

rail_stop_ids <- gtfs$stop_times %>%
  semi_join(rail_trips, by = "trip_id") %>%
  distinct(stop_id)

# Rail trips stop at platforms. Each platform points to its station through
# parent_station, so one station ID per physical station is the parent when it
# exists and the platform itself otherwise.
rail_station_ids <- gtfs$stops %>%
  semi_join(rail_stop_ids, by = "stop_id") %>%
  mutate(station_id = if_else(is.na(parent_station) | parent_station == "",
                              stop_id, parent_station)) %>%
  distinct(station_id)

# Spatial points for the station records, kept only inside the study area.
study_area <- st_as_sfc(study_bbox)

rail_stations <- stops_as_sf(gtfs$stops, crs = 4326) %>%
  filter(stop_id %in% rail_station_ids$station_id) %>%
  distinct(stop_id, .keep_all = TRUE) %>%
  st_filter(study_area) %>%
  select(stop_id, stop_name)

# Eligible park-and-ride stations: rail stations with free daily parking.
park_ride_stations <- rail_stations %>%
  filter(stop_name %in% free_parking_station_names)

# Destination: North Avenue's parent record is named "NORTH AVE STATION".
destination_stop <- rail_stations %>%
  filter(stop_name == "NORTH AVE STATION")

# Transfers and stop times for the assigned date and morning period preparation.
# Travel_times() needs stops that share a name to be within
# 300 m of each other. Some non-rail stops (mostly bus stops) share names but
# are farther apart, so they get numbered names. Rail platforms and stations
# are protected, so their names still match the parking list and destination.
rail_ids_keep <- union(rail_stop_ids$stop_id, rail_station_ids$station_id)

stop_dist <- stop_group_distances(gtfs$stops, by = "stop_name") %>%
  filter(dist_max > 300)

gtfs$stops <- gtfs$stops %>%
  group_by(stop_name) %>%
  mutate(stop_name = if_else(stop_name %in% stop_dist$stop_name &
                               !(stop_id %in% rail_ids_keep),
                             paste0(stop_name, " (", row_number(), ")"),
                             stop_name)) %>%
  ungroup()

# Walking transfers between stops within 200 m, each taking at least 120 s.
gtfs <- gtfsrouter::gtfs_transfer_table(gtfs,
                                        d_limit = 200,
                                        min_transfer_time = 120)

# The full network is kept for routing (no rail or study-area filter here).
# Departures from 7 AM, arrivals by 10 AM, on November 18, 2026.
morning_stop_times <- filter_stop_times(gtfs_obj = gtfs,
                                        extract_date = analysis_date,
                                        min_departure_time = am_start,
                                        max_arrival_time = am_end)

# Checkpoint: these checks should run without errors.
stopifnot(inherits(rail_stations, "sf"))
stopifnot(nrow(rail_stations) > 0)
stopifnot("stop_name" %in% names(rail_stations))
stopifnot(inherits(park_ride_stations, "sf"))
stopifnot(nrow(park_ride_stations) > 0)
stopifnot(all(park_ride_stations$stop_name %in% free_parking_station_names))
stopifnot(nrow(destination_stop) >= 1)

2. Prepare Census Data

Download 2024 ACS 5-year estimates for Census tracts in Fulton, DeKalb, and Clayton counties. Prepare:

• median household income;
• the population variables needed to calculate the proportion of minority population; and
• tract geometry.

Briefly state how you define minority population.

Explanation: Minority population could be defined as everyone who is not non-Hispanic White alone (ACS table B03002). This counts Hispanic or Latino residents of any race and non-Hispanic residents of every other race or combination of races. pct_minority is the total population minus the non-Hispanic White alone population, divided by the total population.

Stored pct_minority as a proportion from 0 to 1. Kept the tract polygons for mapping and created one centroid per tract for routing. Calculated the centroids in a suitable projected CRS, then retained only origins whose centroids fell within the study area using study_bbox.

Create these checkpoint objects:

• census_tracts: tract polygons containing GEOID, hhinc, and pct_minority; and
• origins: tract centroids containing the same identifiers and attributes.

# Download the ACS estimates and calculate the demographic measures.
acs <- get_acs(geography = "tract",
               variables = c(hhinc = "B19013_001",     # median household income
                             pop_total = "B03002_001", # total population
                             nh_white = "B03002_003"), # non-Hispanic White alone
               year = 2024,
               survey = "acs5",
               state = "GA",
               county = c("Fulton", "DeKalb", "Clayton"),
               output = "wide",
               geometry = TRUE,
               progress_bar = FALSE)

acs <- acs %>%
  transmute(GEOID,
            hhinc = hhincE,
            pop_total = pop_totalE,
            nh_white = nh_whiteE,
            # Proportion from 0 to 1; NA where a tract has no residents.
            pct_minority = if_else(pop_total > 0,
                                   (pop_total - nh_white) / pop_total,
                                   NA_real_))

# Create tract centroids and retain origins within the study area.
# Centroids are computed in a projected CRS (UTM 16N, meters), then returned
# to WGS 84 to match the bounding box, stations and road network.
origins <- acs %>%
  st_transform(32616) %>%
  st_centroid() %>%
  st_transform(4326) %>%
  st_filter(study_area) %>%
  select(GEOID, hhinc, pct_minority)

# Keep the corresponding tract polygons and attributes for mapping.
census_tracts <- acs %>%
  st_transform(4326) %>%
  filter(GEOID %in% origins$GEOID) %>%
  select(GEOID, hhinc, pct_minority)

# Checkpoint: these checks should run without errors.
stopifnot(inherits(census_tracts, "sf"))
stopifnot(inherits(origins, "sf"))
stopifnot(all(c("GEOID", "hhinc", "pct_minority") %in% names(census_tracts)))
stopifnot(all(c("GEOID", "hhinc", "pct_minority") %in% names(origins)))
stopifnot(nrow(origins) > 0)

3. Build the OSM Road Network

Download a road network containing motorways, trunk roads, primary roads, secondary roads, and tertiary roads. Include the corresponding connector roads for all five classes so that ramps and connecting segments are retained. Exclude residential roads and other highway classes. Refer to the OSM lab or OSM documentation for the appropriate highway values.

Assign each road segment an assumed speed in miles per hour:

• Motorway: 55
• Trunk: 45
• Primary: 35
• Secondary: 30
• Tertiary: 25
• Connector roads (all five classes): 20

Convert the roads to an undirected sfnetwork and clean it using the methods introduced in the OSM lab. When smoothing the network, 1) preserve road-class changes with require_equal = "highway" so that segments with different assumed speeds are not merged into a single edge and 2) use summarise_attributes = "first" to retain the road class on merged edges.

After cleaning, calculate edge lengths, assign the assumed speeds, and calculate each edge’s travel time. Use consistent units: lengths in meters, speeds in miles per hour, and travel times in minutes. Every edge must have a road class listed above and a valid speed.

Create road_network, an sfnetwork with edge columns named length, speed_mph, and travel_time_min.

# Download the roads within the study area.
road_classes <- c("motorway", "trunk", "primary", "secondary", "tertiary")
link_classes <- paste0(road_classes, "_link")  # connector roads (ramps)

# The public Overpass servers are often busy (HTTP 502/504), so the download
# tries several servers, retries, and saves the result. Later knits reuse the
# saved file instead of downloading again.
osm_cache <- "osm_roads_study_area.rds"

if (file.exists(osm_cache)) {
  osm_roads <- readRDS(osm_cache)
} else {
  overpass_servers <- c("https://overpass-api.de/api/interpreter",
                        "https://overpass.kumi.systems/api/interpreter",
                        "https://overpass.private.coffee/api/interpreter")
  osm_roads <- NULL
  for (attempt in 1:3) {
    for (server in overpass_servers) {
      osmdata::set_overpass_url(server)
      osm_roads <- tryCatch(
        opq(bbox = study_bbox, timeout = 300) %>%
          add_osm_feature(key = "highway", value = c(road_classes, link_classes)) %>%
          osmdata_sf(),
        error = function(e) {
          message("Overpass request failed (", server, "): ", conditionMessage(e))
          NULL
        })
      if (!is.null(osm_roads)) break
    }
    if (!is.null(osm_roads)) break
    Sys.sleep(30)  # wait before trying all servers again
  }
  if (is.null(osm_roads)) stop("All Overpass servers failed; try again later.")
  saveRDS(osm_roads, osm_cache)
}

roads <- osm_roads$osm_lines %>%
  filter(highway %in% c(road_classes, link_classes)) %>%
  select(highway)

# Build and clean the network, preserving road classes.
road_network <- roads %>%
  as_sfnetwork(directed = FALSE) %>%       # undirected: one-way rules ignored
  activate("edges") %>%
  filter(!edge_is_multiple()) %>%           # drop duplicate edges
  filter(!edge_is_loop()) %>%               # drop self-loops
  convert(to_spatial_subdivision) %>%       # connect lines at shared interior points
  convert(to_spatial_smooth,                # remove pseudo-nodes, but...
          require_equal = "highway",        # ...never merge different road classes
          summarise_attributes = "first") %>%
  activate("nodes") %>%
  filter(group_components() == 1)           # keep the largest connected component

# Calculate edge lengths, assign speeds, and calculate travel times.
road_network <- road_network %>%
  activate("edges") %>%
  mutate(length = as.numeric(edge_length()),   # meters
         speed_mph = case_when(
           highway == "motorway"  ~ 55,
           highway == "trunk"     ~ 45,
           highway == "primary"   ~ 35,
           highway == "secondary" ~ 30,
           highway == "tertiary"  ~ 25,
           highway %in% link_classes ~ 20,
           TRUE ~ NA_real_),
         # meters -> miles, divided by mph = hours; x 60 = minutes
         travel_time_min = (length / 1609.344) / speed_mph * 60)

# Every edge must have a listed road class and a valid speed.
stopifnot(!any(is.na(road_network %>% st_as_sf("edges") %>% pull(speed_mph))))

# Checkpoint: these checks should run without errors.
stopifnot(inherits(road_network, "sfnetwork"))
stopifnot(all(c("length", "speed_mph", "travel_time_min") %in%
                names(st_as_sf(road_network, "edges"))))
stopifnot(nrow(road_network %>% st_as_sf("nodes")) > 0)
stopifnot(nrow(road_network %>% st_as_sf("edges")) > 0)

4. Test One Park-and-Ride Trip

Before writing a function, test the complete workflow with one row from origins. Follow the simulation assumptions above. Pass weights argument explicitly to st_network_paths() to find the fastest paths. Use the comments below to organize your test.

Do not treat every empty edge path as a disconnected destination. An empty path can also occur when the origin and destination snap to the same network node.

Create test_trip, a one-row data.frame containing at least:

• GEOID;
• selected_station (the eligible station reached fastest by car);
• drive_time_min: estimated driving time from the tract centroid to the selected station, in minutes;
• transit_time_min: minimum scheduled journey duration from the selected station to North Avenue during the assigned morning window, in minutes; and
• total_time_min: sum of driving and transit travel times, in minutes.

# Choose one origin tract for testing.
test_origin <- origins[1, ]

# Nodes and edge travel times used for routing.
net_nodes <- road_network %>% st_as_sf("nodes")
edge_minutes <- road_network %>% st_as_sf("edges") %>% pull(travel_time_min)

# Snap the origin and stations to their nearest network nodes.
origin_node  <- st_nearest_feature(test_origin, net_nodes)
station_node <- st_nearest_feature(park_ride_stations, net_nodes)
target_nodes <- unique(station_node)  # two stations may share a node

# Calculate travel-time-weighted paths to all eligible parking stations.
paths <- st_network_paths(road_network,
                          from = origin_node,
                          to = target_nodes,
                          weights = "travel_time_min")

# Identify reachable stations and sum the edge travel times for each path.
# (The edge path column is `edge_paths` in older sfnetworks, `edge_path` in 1.0+.)
edge_path_list <- if ("edge_paths" %in% names(paths)) paths$edge_paths else paths$edge_path

node_minutes <- map2_dbl(edge_path_list, target_nodes, function(ep, node) {
  if (node == origin_node) return(0)   # same node: empty path, 0 minutes
  if (length(ep) == 0) return(NA)      # empty path to another node: unreachable
  sum(edge_minutes[ep])
})

drive_options <- park_ride_stations %>%
  st_drop_geometry() %>%
  mutate(drive_time_min = node_minutes[match(station_node, target_nodes)]) %>%
  filter(!is.na(drive_time_min))

# Select the station with the shortest driving time.
best_station <- drive_options %>%
  slice_min(drive_time_min, n = 1, with_ties = FALSE)

# Calculate transit journeys from that station to North Avenue station during the assigned morning window; retain the minimum journey duration.
test_transit <- travel_times(filtered_stop_times = morning_stop_times,
                             stop_name = best_station$stop_name,
                             time_range = 3 * 3600,   # full 7-10 AM window
                             arrival = FALSE,
                             max_transfers = 1)

transit_min <- test_transit %>%
  filter(to_stop_name == destination_stop$stop_name[1]) %>%
  summarise(m = min(as.numeric(travel_time)) / 60) %>%
  pull(m)

# Create test_trip with the origin GEOID, selected station, driving time, transit time, and total time.
test_trip <- data.frame(
  GEOID = test_origin$GEOID,
  selected_station = best_station$stop_name,
  drive_time_min = best_station$drive_time_min,
  transit_time_min = transit_min,
  total_time_min = best_station$drive_time_min + transit_min
)

test_trip

##         GEOID                selected_station drive_time_min transit_time_min
## 1 13089020100 INMAN PARK-REYNOLDSTOWN STATION       6.744808               14
##   total_time_min
## 1       20.74481

# Checkpoint: these checks should run without errors.
required_trip_fields <- c(
  "GEOID", "selected_station", "drive_time_min",
  "transit_time_min", "total_time_min"
)
stopifnot(inherits(test_trip, "data.frame"))
stopifnot(nrow(test_trip) == 1)
stopifnot(all(required_trip_fields %in% names(test_trip)))

The test tract (GEOID 13089020100) reaches Inman Park-Reynoldstown Station fastest, in about 6.7 minutes by car. The fastest scheduled train journey from there to North Avenue between 7 and 10 AM takes 14 minutes, for a total of about 20.7 minutes. This one-row result has the required fields, so the same steps are turned into a function below.

5. Write a Reusable Function

Convert the tested workflow into a function. Keep the function focused on one origin point and return a one-row data.frame with the same fields as test_trip. You may add arguments if your implementation requires them, but keep the arguments shown below and return the required fields.

calculate_park_ride <- function(
    origin_point,
    road_network,
    park_ride_stations,
    morning_stop_times,
    destination_stop,
    max_transfers) {

  # Implement the workflow tested in the previous section.

  # Network nodes and edge travel times (minutes).
  net_nodes <- road_network %>% st_as_sf("nodes")
  edge_minutes <- road_network %>% st_as_sf("edges") %>% pull(travel_time_min)

  # Snap the origin and the eligible stations to their nearest nodes.
  origin_node  <- st_nearest_feature(origin_point, net_nodes)
  station_node <- st_nearest_feature(park_ride_stations, net_nodes)
  target_nodes <- unique(station_node)

  # Fastest paths, weighted by travel time.
  paths <- st_network_paths(road_network,
                            from = origin_node,
                            to = target_nodes,
                            weights = "travel_time_min")
  edge_path_list <- if ("edge_paths" %in% names(paths)) paths$edge_paths else paths$edge_path

  # Driving minutes per target node: 0 if the origin snaps to the same node,
  # NA if no path exists, otherwise the sum of edge travel times.
  node_minutes <- map2_dbl(edge_path_list, target_nodes, function(ep, node) {
    if (node == origin_node) return(0)
    if (length(ep) == 0) return(NA)
    sum(edge_minutes[ep])
  })

  drive_options <- park_ride_stations %>%
    st_drop_geometry() %>%
    mutate(drive_time_min = node_minutes[match(station_node, target_nodes)]) %>%
    filter(!is.na(drive_time_min))

  if (nrow(drive_options) == 0) stop("no eligible station is reachable by road")

  # Station with the shortest driving time (not distance, not total time).
  best_station <- drive_options %>%
    slice_min(drive_time_min, n = 1, with_ties = FALSE)

  # Minimum scheduled journey to North Avenue, 7-10 AM, at most max_transfers.
  transit <- travel_times(filtered_stop_times = morning_stop_times,
                          stop_name = best_station$stop_name,
                          time_range = 3 * 3600,
                          arrival = FALSE,
                          max_transfers = max_transfers)

  transit_min <- transit %>%
    filter(to_stop_name == destination_stop$stop_name[1]) %>%
    summarise(m = min(as.numeric(travel_time)) / 60) %>%
    pull(m)

  if (length(transit_min) == 0 || !is.finite(transit_min)) {
    stop("no transit journey to North Avenue in the morning window")
  }

  data.frame(
    GEOID = origin_point$GEOID,
    selected_station = best_station$stop_name,
    drive_time_min = best_station$drive_time_min,
    transit_time_min = transit_min,
    total_time_min = best_station$drive_time_min + transit_min
  )
}

6. Apply the Function to All Origins

Apply calculate_park_ride() to every row in origins. The iteration code is provided below. If you retain the required object and function names, you do not need to modify it. No additional coding is required in this section; run the provided chunks, inspect the results, and explain any unsuccessful trips.

Create trip_results, with one row per attempted origin and all required fields. The provided loop records an NA row if a calculation raises an error, then continues to the next origin. Test one trip successfully first; the loop also catches programming errors, which should be corrected.

# Provided iteration structure
trip_results_list <- vector("list", nrow(origins))

for (i in seq_len(nrow(origins))) {
  trip_results_list[[i]] <- tryCatch(
    calculate_park_ride(
      origin_point = origins[i, ],
      road_network = road_network,
      park_ride_stations = park_ride_stations,
      morning_stop_times = morning_stop_times,
      destination_stop = destination_stop,
      max_transfers = 1
    ),
    error = function(e) {
      message("Tract ", origins$GEOID[i], ": ", conditionMessage(e))
      data.frame(
        GEOID = origins$GEOID[i],
        selected_station = NA_character_,
        drive_time_min = NA_real_,
        transit_time_min = NA_real_,
        total_time_min = NA_real_
      )
    }
  )
}

trip_results <- bind_rows(trip_results_list)

# Checkpoint: inspect failures before creating the final visualizations.
stopifnot(inherits(trip_results, "data.frame"))
stopifnot(all(required_trip_fields %in% names(trip_results)))
stopifnot(nrow(trip_results) == nrow(origins))

sum(is.na(trip_results$total_time_min))

## [1] 0

If there were any unsuccessful trips, briefly explain why.

Explanation: There were no unsuccessful trips. All 124 origins returned a total travel time, and the checkpoint above reports 0 missing values. Every tract centroid snapped to the connected road network, at least one eligible station was reachable by car from each origin, and each selected station had a scheduled journey to North Avenue between 7 and 10 AM with at most one transfer.

7. Visualize and Interpret the Results

Join trip_results to the Census attributes and tract polygons to prepare the maps and scatterplots.

Create the following four visualizations:

1. An interactive map comparison of household income and total travel time.
2. An interactive map comparison of the proportion of minority population and total travel time.
3. A scatterplot of household income and total travel time, including a fitted trend line.
4. A scatterplot of the proportion of minority population and total travel time, including a fitted trend line.

For each map comparison, you may use two layers in one map, two synchronized maps side by side, or two separate interactive maps. Make both variables easy to identify and compare. If using separate maps, use the same geographic extent and zoom level.

Each visualization should have a clear title and clearly labeled units. All maps must include legends. Maps should make the geographic pattern legible without obscuring the Census geography.

# Map household income and total park-and-ride travel time.
results_map <- census_tracts %>%
  left_join(trip_results, by = "GEOID")

time_pal   <- colorNumeric("YlOrRd", domain = results_map$total_time_min, na.color = "#d9d9d9")
income_pal <- colorNumeric("Blues",  domain = results_map$hhinc,          na.color = "#d9d9d9")

tract_popup <- ~paste0("<b>Tract ", GEOID, "</b><br>",
                       "Median household income: ", dollar(hhinc), "<br>",
                       "Minority population: ", percent(pct_minority, 0.1), "<br>",
                       "Total travel time: ", round(total_time_min, 1), " min<br>",
                       "Parking station: ", selected_station)

map_income <- leaflet(results_map) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(fillColor = ~income_pal(hhinc), fillOpacity = 0.7,
              color = "#555555", weight = 0.6, popup = tract_popup) %>%
  addLegend("bottomright", pal = income_pal, values = ~hhinc,
            title = "Median household<br>income (USD)",
            labFormat = labelFormat(prefix = "$", big.mark = ",")) %>%
  addControl("<b>Median household income</b>", position = "topright")

map_time <- leaflet(results_map) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(fillColor = ~time_pal(total_time_min), fillOpacity = 0.7,
              color = "#555555", weight = 0.6, popup = tract_popup) %>%
  addCircleMarkers(data = park_ride_stations, radius = 4, color = "#08306b",
                   fillOpacity = 1, stroke = FALSE, label = ~stop_name) %>%
  addCircleMarkers(data = destination_stop, radius = 6, color = "#000000",
                   fillOpacity = 1, stroke = FALSE, label = "North Avenue Station") %>%
  addLegend("bottomright", pal = time_pal, values = ~total_time_min,
            title = "Total park-and-ride<br>travel time (minutes)") %>%
  addControl("<b>Total travel time to North Avenue Station</b>", position = "topright")

# Two synchronized maps: same extent and zoom, panning one moves the other.
leafsync::sync(map_income, map_time)

# Map the proportion of minority population and total park-and-ride travel time.
minority_pal <- colorNumeric("Purples", domain = c(0, 1), na.color = "#d9d9d9")

map_minority <- leaflet(results_map) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(fillColor = ~minority_pal(pct_minority), fillOpacity = 0.7,
              color = "#555555", weight = 0.6, popup = tract_popup) %>%
  addLegend("bottomright", pal = minority_pal, values = c(0, 1),
            title = "Minority population<br>(proportion, 0-1)") %>%
  addControl("<b>Proportion of minority population</b>", position = "topright")

leafsync::sync(map_minority, map_time)

# Plot household income against total park-and-ride travel time.
results_plot <- results_map %>% st_drop_geometry()

# Seven tracts have no median income estimate in the ACS, so they are left out
# of this plot and its correlation (117 of 124 tracts are shown).
ggplot(filter(results_plot, !is.na(hhinc)), aes(x = hhinc, y = total_time_min)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE, color = "red") +
  scale_x_continuous(labels = dollar) +
  labs(title = "Median household income and park-and-ride travel time",
       x = "Median household income (USD)",
       y = "Total travel time to North Avenue Station (minutes)") +
  theme_bw()

# Plot the proportion of minority population against total park-and-ride travel time.
ggplot(results_plot, aes(x = pct_minority, y = total_time_min)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE, color = "red") +
  scale_x_continuous(labels = percent, limits = c(0, 1)) +
  labs(title = "Minority population and park-and-ride travel time",
       x = "Minority population (share of tract residents)",
       y = "Total travel time to North Avenue Station (minutes)") +
  theme_bw()

# Correlations to support the findings (pairwise complete tracts only).
cor.test(results_plot$hhinc, results_plot$total_time_min)

## 
##  Pearson's product-moment correlation
## 
## data:  results_plot$hhinc and results_plot$total_time_min
## t = 2.5052, df = 115, p-value = 0.01364
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.04793207 0.39279758
## sample estimates:
##       cor 
## 0.2274856

cor.test(results_plot$pct_minority, results_plot$total_time_min)

## 
##  Pearson's product-moment correlation
## 
## data:  results_plot$pct_minority and results_plot$total_time_min
## t = -3.4969, df = 122, p-value = 0.0006573
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4539854 -0.1325670
## sample estimates:
##        cor 
## -0.3018288

Findings

In no more than 200 words, summarize the most important spatial and statistical patterns. Discuss whether the results suggest disparities in park-and-ride accessibility and identify at least two limitations of the simulation. Separate what the analysis shows from explanations that would require additional evidence.

Findings: Explanations Simulated park-and-ride times to North Avenue range from 7.2 to 26.5 minutes (mean 15.0) across 124 tracts, and every origin returned a result. Times are shortest in the north near Lindbergh Center (about 7 to 9 minutes) and across the southwest, where West End, Oakland City, West Lake and Ashby are 9 to 12 minutes from North Avenue by train. They are longest on the east side and at the western edge, where trips rely on East Lake, Inman Park-Reynoldstown or Hamilton E Holmes (20 to 26 minutes).

Travel time rises slightly with income (r = 0.23, p = 0.014; 117 tracts with income data) and falls with the minority share (r = -0.30, p < 0.001). In this simulation, lower-income and majority-minority tracts, concentrated in the southwest, do not face longer park-and-ride trips. These correlations are weak and show association only.

The pattern reflects where free-parking stations sit; explaining it would need evidence on station siting and travel behavior. Limitations: assumed speeds ignore congestion, signals and one-way streets; tract centroids stand in for residents; drive and train times are added without timing the connection; and the method assumes car access, which many lower-income households lack.

Your analysis and written responses end here. The remaining sections display a result summary, selected trip results, and session information. Run them and include their output in your HTML submission. Do not edit or delete the code chunks below.

Results Summary

# Provided output: do not edit or delete.
result_summary <- data.frame(
  rail_stations = nrow(rail_stations),
  eligible_parking_stations = nrow(park_ride_stations),
  origins = nrow(origins),
  result_rows = nrow(trip_results),
  duplicate_result_geoids = sum(duplicated(trip_results$GEOID)),
  missing_total_times = sum(is.na(trip_results$total_time_min))
)

knitr::kable(result_summary)

rail_stations	eligible_parking_stations	origins	result_rows	duplicate_result_geoids	missing_total_times
21	9	124	124	0	0

# Provided output: do not edit or delete.
travel_time_summary <- trip_results %>%
  dplyr::summarise(
    min_min = min(total_time_min, na.rm = TRUE),
    max_min = max(total_time_min, na.rm = TRUE),
    mean_min = mean(total_time_min, na.rm = TRUE)
  )

knitr::kable(travel_time_summary, digits = 2, caption = "Simulated total travel time (minutes)")

Simulated total travel time (minutes)

min_min	max_min	mean_min
7.16	26.45	14.99

# Provided output: do not edit or delete.
eligible_station_table <- park_ride_stations %>%
  sf::st_drop_geometry() %>%
  dplyr::select(stop_id, stop_name) %>%
  dplyr::arrange(stop_name)

knitr::kable(eligible_station_table)

stop_id	stop_name
510002	ASHBY STATION
510012	EAST LAKE STATION
510014	EDGEWOOD-CANDLER PARK STATION
510017	HAMILTON E HOLMES STATION
510019	INMAN PARK-REYNOLDSTOWN STATION
510024	LINDBERGH CENTER STATION
510029	OAKLAND CITY STATION
510033	WEST END STATION
510034	WEST LAKE STATION

# Provided output: do not edit or delete.
station_selections <- trip_results %>%
  dplyr::count(selected_station, name = "origin_count", sort = TRUE)

knitr::kable(station_selections)

selected_station	origin_count
INMAN PARK-REYNOLDSTOWN STATION	35
LINDBERGH CENTER STATION	27
ASHBY STATION	24
WEST END STATION	15
EDGEWOOD-CANDLER PARK STATION	5
OAKLAND CITY STATION	5
WEST LAKE STATION	5
EAST LAKE STATION	4
HAMILTON E HOLMES STATION	4

# Provided output: do not edit or delete.
set.seed(2026)
selected_geoids <- origins %>%
  sf::st_drop_geometry() %>%
  dplyr::distinct(GEOID) %>%
  dplyr::arrange(GEOID) %>%
  dplyr::slice_sample(n = 20) %>%
  dplyr::arrange(GEOID)

selected_results <- selected_geoids %>%
  dplyr::left_join(trip_results, by = "GEOID")

knitr::kable(selected_results, digits = 2)

GEOID	selected_station	drive_time_min	transit_time_min	total_time_min
13089020500	EDGEWOOD-CANDLER PARK STATION	1.39	16	17.39
13089020902	INMAN PARK-REYNOLDSTOWN STATION	5.57	14	19.57
13089022405	EAST LAKE STATION	7.45	19	26.45
13121000700	ASHBY STATION	3.87	9	12.87
13121001203	INMAN PARK-REYNOLDSTOWN STATION	4.77	14	18.77
13121001205	INMAN PARK-REYNOLDSTOWN STATION	5.08	14	19.08
13121001502	INMAN PARK-REYNOLDSTOWN STATION	3.88	14	17.88
13121001600	INMAN PARK-REYNOLDSTOWN STATION	2.31	14	16.31
13121001801	INMAN PARK-REYNOLDSTOWN STATION	5.32	14	19.32
13121002802	INMAN PARK-REYNOLDSTOWN STATION	3.88	14	17.88
13121005504	OAKLAND CITY STATION	5.34	11	16.34
13121006900	INMAN PARK-REYNOLDSTOWN STATION	5.03	14	19.03
13121008301	WEST LAKE STATION	2.53	12	14.53
13121008400	ASHBY STATION	2.99	9	11.99
13121008702	ASHBY STATION	6.72	9	15.72
13121008801	HAMILTON E HOLMES STATION	10.94	15	25.94
13121009103	LINDBERGH CENTER STATION	5.78	7	12.78
13121009106	LINDBERGH CENTER STATION	3.71	7	10.71
13121011802	ASHBY STATION	3.65	9	12.65
13121012000	WEST END STATION	2.66	9	11.66

Reproducibility

sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.6.2
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] scales_1.4.0      leafsync_0.1.0    leaflet_2.2.3     tidygraph_1.3.1  
##  [5] sfnetworks_0.6.6  osmdata_0.4.1     tidycensus_1.8.1  gtfsrouter_0.1.4 
##  [9] tidytransit_1.8.0 sf_1.1-3          lubridate_1.9.5   forcats_1.0.1    
## [13] stringr_1.6.0     dplyr_1.2.1       purrr_1.2.2       readr_2.2.0      
## [17] tidyr_1.3.2       tibble_3.3.1      ggplot2_4.0.3     tidyverse_2.0.0  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1        farver_2.1.2            S7_0.2.2               
##  [4] fastmap_1.2.0           digest_0.6.39           timechange_0.4.0       
##  [7] lifecycle_1.0.5         magrittr_2.0.5          compiler_4.6.1         
## [10] rlang_1.3.0             sass_0.4.10             tools_4.6.1            
## [13] igraph_2.3.4            data.table_1.18.6.1     knitr_1.52             
## [16] labeling_0.4.3          htmlwidgets_1.6.4       classInt_0.4-11        
## [19] curl_8.0.0              xml2_1.6.0              RColorBrewer_1.1-3     
## [22] gtfsio_1.2.2            KernSmooth_2.23-27      withr_3.0.3            
## [25] grid_4.6.1              e1071_1.7-17            cli_3.6.6              
## [28] rmarkdown_2.32          crayon_1.5.3            generics_0.1.4         
## [31] otel_0.2.0              rstudioapi_0.19.0       httr_1.4.9             
## [34] tzdb_0.5.0              DBI_1.3.0               cachem_1.1.0           
## [37] proxy_0.4-29            splines_4.6.1           rvest_1.0.5            
## [40] s2_1.1.11               vctrs_0.7.3             tigris_2.2.1           
## [43] Matrix_1.7-6            jsonlite_2.0.0          hms_1.1.4              
## [46] crosstalk_1.2.2         jquerylib_0.1.4         units_1.0-1            
## [49] glue_1.8.1              lwgeom_0.2-17           leaflet.providers_3.0.0
## [52] stringi_1.8.9           gtable_0.3.6            pillar_1.11.1          
## [55] rappdirs_0.3.4          htmltools_0.5.9         R6_2.6.1               
## [58] httr2_1.3.0             wk_0.9.5                sfheaders_0.4.5        
## [61] evaluate_1.0.5          lattice_0.23-1          bslib_0.12.0           
## [64] class_7.3-24            Rcpp_1.1.2              geodist_0.1.1          
## [67] zip_3.0.2               uuid_1.2-2              nlme_3.1-171           
## [70] mgcv_1.9-4              xfun_0.61               fs_2.1.0               
## [73] pkgconfig_2.0.3