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")

service_date <- as.Date("2026-11-18")
time_range <- c(7 * 3600, 10 * 3600)

# 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_stations <- gtfs$routes |> 
  filter(route_type == 1) |> 
  select(route_id) |> 
  inner_join(gtfs$trips, by="route_id") |> # Performing series of joins as needed
  select(trip_id) |> 
  inner_join(gtfs$stop_times, by ="trip_id") |> 
  select(stop_id) |> 
  inner_join(gtfs$stops, by ="stop_id") |> 
  select(stop_id, parent_station, stop_name, stop_lat, stop_lon) |> 
  distinct(parent_station, .keep_all = TRUE) |> 
  st_as_sf(coords = c("stop_lon", "stop_lat"), crs = st_crs(4326)) |> 
  st_filter(st_as_sfc(study_bbox)) # Filtering to study area

park_ride_stations <- rail_stations |> 
  filter(stop_name %in% free_parking_station_names) # Filtering acceptable stations

destination_stop <- rail_stations |> 
  filter(stop_name == "NORTH AVE STATION") # Declaring destination stop for the study

# Prepare transfers and stop times for the assigned date and morning period.
gtfs <- gtfs |> 
  gtfs_transfer_table(
  d_limit = 200,
  min_transfer_time = 120
) # Using the gtfs_transfer_table to prepare a table of transit transfers

morning_stop_times <- filter_stop_times(
  gtfs_obj = gtfs,
  extract_date = service_date,
  min_departure_time = time_range[1],
  max_arrival_time = time_range[2]
) # Filtering stopping time in the 7AM to 10AM range

transit_times <- morning_stop_times |> 
  travel_times(
    stop_name = destination_stop$stop_name,
    time_range = time_range,
    arrival = TRUE,
    max_transfers = 1,
    stop_dist_check = FALSE) # Calculating transit times from the selected stations during the given time range, allowing a max of 1 transfers

# 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.

Your explanation:

For this study, a simple calculation of minority population as the non-white population is sufficient. Metro-Atlanta is composed of diverse social groups, and building up minority population as a metric by adding group populations is likely to exclude some populations, or get dominated by only a few groups. Hence, I used the difference between the total population and the non-hispanic white population as the basis of numerically defining minority population.

Store pct_minority as a proportion from 0 to 1. Keep the tract polygons for mapping and create one centroid per tract for routing. Calculate centroids in a suitable projected CRS, then retain only origins whose centroids fall 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.

census_tracts <- get_acs(
  geography = "tract",
  state = "GA",
  county = c("Fulton", "DeKalb", "Clayton"),
  year = 2024,
  survey = "acs5",
  variables = c(
    hhinc = "B19013_001",
    total_pop = "B03002_001",
    white_nh = "B03002_003"
  ),
  geometry = TRUE,
  output = "wide"
) |>
  st_transform(crs = 4326) |>
  mutate(
    hhinc = hhincE,
    pct_minority = 1 - white_nhE / total_popE # Calculating minority population, as described above
  ) |>
  select(GEOID, hhinc, pct_minority) # Filtering the columns needed for the study

# Create tract centroids and retain origins within the study area.

origins <- census_tracts |>
  st_transform(crs = 26916) |>
  st_centroid() |>
  st_transform(4326) |>
  st_filter(st_as_sfc(study_bbox)) |>
  st_transform(26916)

# Keep the corresponding tract polygons and attributes for mapping.

# 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.

osm_data <- opq(
  bbox = study_bbox,
  timeout = 180
) |> 
  add_osm_feature(
    key = "highway",
    value = c("motorway", "motorway_link",
                            "trunk", "trunk_link", 
                            "primary", "primary_link",
                            "secondary", "secondary_link",
                            "tertiary", "tertiary_link")
  ) |> 
  osmdata_sf() |> 
  osm_poly2line()

# Build and clean the network, preserving road classes.
roads <- osm_data$osm_lines |> 
  st_transform(26916) |> 
  select(osm_id, name, highway, geometry) |> 
  mutate(
    speed = case_when(
      grepl("_link$", highway) ~ 20,  # Assigning speed values as defined above
      highway == "motorway" ~ 55,
      highway == "trunk" ~ 45,
      highway == "primary" ~ 35,
      highway == "secondary" ~ 30,
      highway == "tertiary" ~ 25))
# Cleaning data:
road_network <- roads |> 
  as_sfnetwork(directed = FALSE) |> 
  activate("edges") |> 
  filter(!edge_is_multiple()) |> # Filtering out duplicates
  filter(!edge_is_loop()) |> # Filtering out edges with identical origin and destination nodes
  convert(to_spatial_subdivision) |> # Creating missing nodes at intersections 
  convert(to_spatial_smooth, require_equal = "highway", summarise_attributes = "first") # Removing consecutive links of the same class while preserving information

# Calculate edge lengths, assign speeds, and calculate travel times.

road_network <- road_network |> 
  activate("edges") |> 
  mutate(len = edge_length()) |> 
  mutate(
  speed_mph = speed,
  length = as.numeric(edge_length()),
  travel_time_min = (len / 1609.344) / speed_mph * 60) # Converting length to miles, dividing by speed, and converting time to minutes

# 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)

# Extract cleaned network edges for mapping
road_edges <- road_network |>
  activate("edges") |>
  st_as_sf()

#Creating sanity-check map:
tmap_mode("view")

tm_shape(road_edges) +
  tm_lines(
    col = "travel_time_min",
    lwd = 2,
    col.legend = tm_legend("travel_time_min"),
    popup.vars = c(
      "Road Name" = "name",
      "Road Class" = "highway",
      "Speed (mph)" = "speed_mph",
      "Length (m)" = "length",
      "Travel Time (min)" = "travel_time_min"
    )
  )

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.

selected_tract <- origins |> 
  filter(GEOID == "13121005300")

# Calculate travel-time-weighted paths to all eligible parking stations.

paths_to_stations <- st_network_paths(road_network, from = selected_tract, to = st_transform(park_ride_stations, 26916), weights = "travel_time_min", type = "shortest")

# Identify reachable stations and sum the edge travel times for each path.

# Creating reusable variible for edge travel times
edge_times <- road_network |>
  activate("edges") |>
  pull(travel_time_min)

station_paths <- park_ride_stations |> 
  select(stop_name) |> 
  bind_cols(paths_to_stations) |> # Connecting station data to paths to corresponding stations
  mutate(
    drive_time = mapply(
      function(edges, nodes) {
        if (length(edges) > 0) { # Summing up edge travel times when one or more edges are preseny in the path
          sum(edge_times[edges])
        } else if (length(nodes) == 1) {0} # If the number of nodes is zero, the tract likely snapped to the station directly, so the driving time is 0
      else {NA_real_}}, edge_paths, node_paths # Disconnected or missing values are handled as NA
    )
  ) |> 
  filter(!is.na(drive_time)) # Filtering out NA data

# Select the station with the shortest driving time.
nearest_station <- station_paths |> 
  arrange(drive_time) |> 
  slice(1)

# Calculate transit journeys from that station to North Avenue station during the assigned morning window; retain the minimum journey duration.
drive_time_min <- nearest_station$drive_time

transit_time_from_nearest_seconds <- transit_times |> 
  filter(from_stop_name == nearest_station$stop_name) |> 
  pull(travel_time)

transit_time_min <- transit_time_from_nearest_seconds / 60.0

total_time_min <- transit_time_min + drive_time_min # Adding the two times

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

# 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)))

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) {
  
  transit_times <- morning_stop_times |> 
  travel_times(
    stop_name = destination_stop$stop_name,
    time_range = c(25200, 36000),
    arrival = TRUE,
    max_transfers = 1,
    stop_dist_check = FALSE)
  
  paths_to_stations <- st_network_paths(road_network, from = origin_point, to = st_transform(
    park_ride_stations, 26916), weights = "travel_time_min", type = "shortest")
  
  edge_times <- road_network |>
  activate("edges") |>
  pull(travel_time_min)
  
  station_paths <- park_ride_stations |> 
  select(stop_name) |> 
  bind_cols(paths_to_stations) |> 
  mutate(
    drive_time = mapply(
      function(edges, nodes) {
        if (length(edges) > 0) {
          sum(edge_times[edges])
        } else if (length(nodes) == 1) {0}
      else {NA_character_}}, edge_paths, node_paths
    )
  ) |> 
  filter(!is.na(drive_time))
  
  nearest_station <- station_paths |> 
  arrange(drive_time) |> 
  slice(1)

  drive_time_min <- nearest_station$drive_time
  
  transit_time_from_nearest_seconds <- transit_times |> 
    filter(from_stop_name == nearest_station$stop_name) |> 
    pull(travel_time)
  
  transit_time_min <- transit_time_from_nearest_seconds / 60.0
  
  total_time_min <- transit_time_min + drive_time_min
  
  # Creating the required data frame as the output
  trip <- data.frame(
    GEOID = origin_point$GEOID,
    selected_station = nearest_station$stop_name,
    drive_time_min = drive_time_min,
    transit_time_min = transit_time_min,
    total_time_min = total_time_min
  )
  
  return(trip)

}

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] 2

# Mapping tracts with missing data:
study_area <- st_as_sfc(study_bbox)

na_tracts <- census_tracts |>
  semi_join(
    trip_results |> filter(is.na(total_time_min)),
    by = "GEOID"
  ) |>
  st_transform(4326)

road_edges <- road_network |>
  activate("edges") |>
  st_as_sf() |>
  st_transform(4326) |>
  st_crop(study_bbox)

map_origins <- origins |>
  st_transform(4326) |>
  st_filter(study_area)

# Creating a mappable transformed shape object for stations, which is also used later
map_stations <- park_ride_stations |>
  st_transform(4326) |>
  st_filter(study_area)

tmap_mode("view")

tm_basemap("OpenStreetMap.HOT") +
  tm_shape(study_area) +
    tm_borders()+
  tm_shape(road_edges) +
  tm_lines(
    col = "steelblue",
    col_alpha = 1,
    lwd = 0.8
  ) +
  tm_shape(na_tracts) +
  tm_polygons(
    fill = "red",
    fill_alpha = 0.2,
    col = "darkred",
    lwd = 1,
    popup.vars = c("GEOID")) +
  tm_shape(map_origins) +
  tm_dots(
    fill = "gold",
    col = "black",
    size = 0.5,
    popup.vars = c("GEOID")
  ) +
  tm_shape(map_stations) +
  tm_symbols(
    fill = "limegreen",
    col = "black",
    shape = 21,
    size = 0.8,
    popup.vars = c("Station" = "stop_name")
  )

If there were any unsuccessful trips, briefly explain why.

Two Census tracts near Emory University returned missing travel times because no eligible park-and-ride stations were reachable through the constructed OSM road network. Zooming into the map, it can be seen that the in-network street near the centroids for these two tracts gets cutoff from the rest of the network by the bounding box. Hence, no feasible paths to any MARTA stations are found. This could either be addressed by letting the centroid snap to thee second-closest node that is not disconnected from the network, or by widening the buffer slightly.

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.

tmap_mode("view")

trip_census <- trip_results |> 
  left_join(census_tracts, by = "GEOID") |> 
  st_as_sf(crs = 4326)

trip_census <- trip_census |> 
  mutate(hhinc_10k = (hhinc/10000)) |> 
  mutate(pct_minority_100 = pct_minority*100)

income_map <- tm_shape(trip_census) +
  tm_polygons(
    fill = "hhinc_10k",
    fill.scale = tm_scale_continuous(values = "Blues"),
    fill.legend = tm_legend(
      "Median Household Income (10K$)",
      text.size = 0.7,
      title.size = 0.85,
      width = 20,
      height = 0.5,
      orientation = "landscape"
    ),
    col = "gray35",
    lwd = 0.5,
    popup.vars = c(
                  "Census Tract" = "GEOID",
                  "Household Income ($)" = "hhinc",
                  "Minority Population (%)" = "pct_minority_100",
                  "Total Travel Time (min)" = "total_time_min",
                  "Selected Station" = "selected_station"
                )
  ) +
  tm_shape(map_stations) +
  tm_symbols(
    fill = "yellow",
    col = "black",
    shape = 21,
    size = 0.8,
    popup.vars = c("Station" = "stop_name")
  ) + tm_shape(destination_stop) + tm_symbols(fill = "red",
    col = "black",
    shape = 21,
    size = 0.8,
    popup = "stop_name") +
  tm_layout(
    legend.position = c("left", "bottom")
  )
travel_map <- tm_shape(trip_census) +
  tm_polygons(
    fill = "total_time_min",
    fill.scale = tm_scale_continuous(values = "Oranges"),
    fill.legend = tm_legend(
      "Total Travel Time (min)",
      text.size = 0.7,
      title.size = 0.85,
      width = 20,
      height = 0.5,
      orientation = "landscape"
    ),
    col = "gray35",
    lwd = 0.5,
    popup.vars = c(
  "Census Tract" = "GEOID",
  "Household Income ($)" = "hhinc",
  "Minority Population (%)" = "pct_minority_100",
  "Total Travel Time (min)" = "total_time_min",
  "Selected Station" = "selected_station"
)
  ) +
  tm_shape(map_stations) +
  tm_symbols(
    fill = "yellow",
    col = "black",
    shape = 21,
    size = 0.8,
    popup.vars = c("Station" = "stop_name")
  ) + tm_shape(destination_stop) + tm_symbols(fill = "red",
    col = "black",
    shape = 21,
    size = 0.8,
    popup = "stop_name") +
  tm_layout(
    legend.position = c("left", "bottom")
  )

tmap_arrange(income_map, travel_map, ncol = 2, sync = TRUE)

# Map the proportion of minority population and total park-and-ride travel time.

tmap_mode("view")

minority_map <- tm_shape(trip_census) +
  tm_polygons(
    fill = "pct_minority_100",
    fill.scale = tm_scale_continuous(values = "Greens"),
    fill.legend = tm_legend(
      "Minority Population (%)",
      text.size = 0.7,
      title.size = 0.85,
      width = 20,
      height = 0.5,
      orientation = "landscape"
    ),
    col = "gray35",
    lwd = 0.5,
    popup.vars = c(
  "Census Tract" = "GEOID",
  "Household Income ($)" = "hhinc",
  "Minority Population (%)" = "pct_minority_100",
  "Total Travel Time (min)" = "total_time_min",
  "Selected Station" = "selected_station"
)
  ) +
  tm_shape(map_stations) +
  tm_symbols(
    fill = "yellow",
    col = "black",
    shape = 21,
    size = 0.8,
    popup.vars = c("Station" = "stop_name")
  ) + tm_shape(destination_stop) + tm_symbols(fill = "red",
    col = "black",
    shape = 21,
    size = 0.8,
    popup = "stop_name") +
  tm_layout(
    legend.position = c("left", "bottom")
  )

tmap_arrange(minority_map, travel_map, ncol = 2, sync = TRUE)

# Plot household income against total park-and-ride travel time.

ggplot(trip_census, aes(x = hhinc, y = total_time_min))+
  geom_point(alpha = 0.9, size = 1) +
  geom_smooth(method = "lm", se = TRUE)+
  labs(
    x = "Median Household Income ($)",
    y = "Total Park-and-Ride Travel Time (min)",
    title = "Median Household Income vs Total Park-and-Ride Travel Time by Tract"
  ) +
  theme_minimal()+
  stat_cor(r.accuracy = 0.01, p.accuracy = 0.001)

# Plot the proportion of minority population against total park-and-ride travel time.

ggplot(trip_census, aes(x = pct_minority_100, y = total_time_min))+
  geom_point(alpha = 0.9, size = 1) +
  geom_smooth(method = "lm", se = TRUE)+
  labs(
    x = "Percentage of Minority Population",
    y = "Total Park-and-Ride Travel Time (min)",
    title = "Minority Population Proportion vs Total Park-and-Ride Travel Time by Tract"
  ) +
  theme_minimal()+
  stat_cor(r.accuracy = 0.01, p.accuracy = 0.001)

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.

Your findings: There is a weak positive correlation between median household income and total travel time and a weak negative correlation between minority population percentage and total travel time. Minority dominated and lower-income tracts in the southwest area of the study seem to have the most efficient connectivity. The study indicates general spatial trends, but the two demographic variables studied do not reliably explain the spatial distribution and efficiency of park-and-ride accessibility. The most significant limitation of this methodology is the travel time calculation. The simulation requires each node in the network to first connect to a station with parking, and then add the transit time. This causes tracts that are close to North Avenue to require going farther out to a MARTA station, and then transit. This assumption is not realistic, as it doesn’t account for direct walking or driving, and skews the statistical relationships revealed by the plots. The analysis also relies on simplified driving speeds, scheduled transit durations, and a geographically constrained road network. Given that the constrained set of tracts is already close to North Avenue, the results do not really reveal insight that would be useful in assessing equitable access based on actual transit demand to North Avenue.

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	2

# 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.51	25.86	14.96

# 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
908963	ASHBY STATION
908566	EAST LAKE STATION
908699	EDGEWOOD-CANDLER PARK STATION
907933	HAMILTON E HOLMES STATION
908437	INMAN PARK-REYNOLDSTOWN STATION
908055	LINDBERGH CENTER STATION
908301	OAKLAND CITY STATION
908704	WEST END STATION
908023	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
ASHBY STATION	25
LINDBERGH CENTER STATION	25
WEST END STATION	15
EDGEWOOD-CANDLER PARK STATION	5
OAKLAND CITY STATION	5
WEST LAKE STATION	5
HAMILTON E HOLMES STATION	4
EAST LAKE STATION	3
NA	2

# 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	NA	NA	NA	NA
13121000700	ASHBY STATION	3.87	9	12.87
13121001203	INMAN PARK-REYNOLDSTOWN STATION	4.78	14	18.78
13121001205	INMAN PARK-REYNOLDSTOWN STATION	5.09	14	19.09
13121001502	INMAN PARK-REYNOLDSTOWN STATION	3.87	14	17.87
13121001600	INMAN PARK-REYNOLDSTOWN STATION	2.30	14	16.30
13121001801	INMAN PARK-REYNOLDSTOWN STATION	5.32	14	19.32
13121002802	INMAN PARK-REYNOLDSTOWN STATION	3.89	14	17.89
13121005504	OAKLAND CITY STATION	5.35	11	16.35
13121006900	INMAN PARK-REYNOLDSTOWN STATION	5.02	14	19.02
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.86	15	25.86
13121009103	LINDBERGH CENTER STATION	5.92	7	12.92
13121009106	LINDBERGH CENTER STATION	3.86	7	10.86
13121011802	ASHBY STATION	3.65	9	12.65
13121012000	WEST END STATION	2.67	9	11.67

Reproducibility

sessionInfo()

## R version 4.6.0 (2026-04-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] knitr_1.51        units_1.0-1       gtfsrouter_0.1.4  tidytransit_1.8.0
##  [5] tidygraph_1.3.1   sfnetworks_0.6.6  osmdata_0.4.1     ggpubr_1.0.0     
##  [9] tmap_4.4-1        tigris_2.2.1      tidycensus_1.8.1  sf_1.1-2         
## [13] lubridate_1.9.5   forcats_1.0.1     stringr_1.6.0     dplyr_1.2.1      
## [17] purrr_1.2.2       readr_2.2.0       tidyr_1.3.2       tibble_3.3.1     
## [21] ggplot2_4.0.3     tidyverse_2.0.0  
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.3.0               httr2_1.3.0             tmaptools_3.3          
##   [4] s2_1.1.11               logger_0.4.3            rlang_1.3.0            
##   [7] magrittr_2.0.5          otel_0.2.0              e1071_1.7-17           
##  [10] compiler_4.6.0          mgcv_1.9-4              png_0.1-9              
##  [13] sfheaders_0.4.5         vctrs_0.7.3             rvest_1.0.5            
##  [16] pkgconfig_2.0.3         wk_0.9.5                crayon_1.5.3           
##  [19] fastmap_1.2.0           backports_1.5.1         labeling_0.4.3         
##  [22] lwgeom_0.2-17           leafem_0.2.5            rmarkdown_2.31         
##  [25] tzdb_0.5.0              spacesXYZ_1.6-0         xfun_0.59              
##  [28] cachem_1.1.0            jsonlite_2.0.0          yyjsonr_0.1.22         
##  [31] uuid_1.2-2              gtfsio_1.2.2            terra_1.9-46           
##  [34] broom_1.0.13            parallel_4.6.0          R6_2.6.1               
##  [37] bslib_0.11.0            stringi_1.8.7           RColorBrewer_1.1-3     
##  [40] car_3.1-5               jquerylib_0.1.4         stars_0.7-3            
##  [43] Rcpp_1.1.2              base64enc_0.1-6         leaflet.providers_3.0.0
##  [46] Matrix_1.7-5            splines_4.6.0           igraph_2.3.4           
##  [49] timechange_0.4.0        tidyselect_1.2.1        rstudioapi_0.19.0      
##  [52] abind_1.4-8             leafgl_0.2.4            maptiles_0.12.0        
##  [55] codetools_0.2-20        curl_7.1.0              lattice_0.22-9         
##  [58] leafsync_0.1.0          withr_3.0.3             S7_0.2.2               
##  [61] evaluate_1.0.5          geodist_0.1.1           proxy_0.4-29           
##  [64] zip_3.0.2               xml2_1.6.0              pillar_1.11.1          
##  [67] carData_3.0-6           KernSmooth_2.23-26      generics_0.1.4         
##  [70] sp_2.2-3                hms_1.1.4               scales_1.4.0           
##  [73] class_7.3-23            glue_1.8.1              tools_4.6.0            
##  [76] leaflegend_1.2.8        data.table_1.18.4       ggsignif_0.6.4         
##  [79] fs_2.1.0                XML_3.99-0.24           grid_4.6.0             
##  [82] crosstalk_1.2.2         colorspace_2.1-3        nlme_3.1-169           
##  [85] cols4all_0.10           raster_3.6-32           Formula_1.2-6          
##  [88] cli_3.6.6               rappdirs_0.3.4          gtable_0.3.6           
##  [91] rstatix_1.1.0           sass_0.4.10             digest_0.6.39          
##  [94] classInt_0.4-11         htmlwidgets_1.6.4       farver_2.1.2           
##  [97] htmltools_0.5.9         lifecycle_1.0.5         leaflet_2.2.3          
## [100] httr_1.4.8