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("tidytransit")
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(gtfsrouter)## Warning: package 'gtfsrouter' was built under R version 4.5.2## Registered S3 method overwritten by 'gtfsrouter':
## method from
## summary.gtfs gtfsiolibrary(sf)## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUElibrary(tidycensus)
library(mapview)
library(leafsync)
library(leaflet)
library(dbscan)##
## Attaching package: 'dbscan'## The following object is masked from 'package:stats':
##
## as.dendrogramlibrary(sfnetworks)
library(tigris)## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.library(tidygraph)##
## Attaching package: 'tidygraph'## The following object is masked from 'package:stats':
##
## filterlibrary(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(osmdata)## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyrightlibrary(here)## here() starts at C:/Users/HojungYu/OneDrive/GT_Semester3/Urban_Analytics/major_1library(magrittr)
library(purrr)##
## Attaching package: 'purrr'## The following object is masked from 'package:magrittr':
##
## set_nameslibrary(tmap)
library(ggplot2)
library(units)## udunits database from C:/Users/HojungYu/AppData/Local/R/win-library/4.5/units/share/udunits/udunits2.xmlStep 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 //////////////////////////////////////////////////////////////////////
head(gtfs)## $trips
## # A tibble: 58,695 × 10
## route_id service_id trip_id trip_headsign trip_short_name direction_id
## <chr> <chr> <chr> <chr> <chr> <int>
## 1 26898 2 10843793 BLUE EASTBOUND TO … "" 0
## 2 26898 2 10843794 BLUE EASTBOUND TO … "" 0
## 3 26898 2 10843795 BLUE EASTBOUND TO … "" 0
## 4 26898 2 10843796 BLUE EASTBOUND TO … "" 0
## 5 26898 2 10843843 BLUE EASTBOUND TO … "" 0
## 6 26898 2 10843797 BLUE EASTBOUND TO … "" 0
## 7 26898 2 10843798 BLUE EASTBOUND TO … "" 0
## 8 26898 2 10843799 BLUE EASTBOUND TO … "" 0
## 9 26898 2 10843800 BLUE EASTBOUND TO … "" 0
## 10 26898 2 10843801 BLUE EASTBOUND TO … "" 0
## # ℹ 58,685 more rows
## # ℹ 4 more variables: block_id <chr>, shape_id <chr>,
## # wheelchair_accessible <int>, bikes_allowed <int>
##
## $agency
## # A tibble: 1 × 7
## agency_id agency_name agency_url agency_timezone agency_lang agency_phone
## <chr> <chr> <chr> <chr> <chr> <chr>
## 1 MARTA Metropolitan At… https://w… America/New_Yo… en 404-848-5000
## # ℹ 1 more variable: agency_fare_url <chr>
##
## $calendar
## # A tibble: 71 × 10
## service_id monday tuesday wednesday thursday friday saturday sunday
## <chr> <int> <int> <int> <int> <int> <int> <int>
## 1 20 0 0 0 0 0 0 0
## 2 28 0 0 0 0 0 0 0
## 3 39 0 0 0 0 0 0 0
## 4 40 0 0 0 0 0 0 0
## 5 41 0 0 0 0 0 0 0
## 6 42 0 0 0 0 0 0 0
## 7 24 0 0 0 0 0 0 0
## 8 25 0 0 0 0 0 0 0
## 9 26 0 0 0 0 0 0 0
## 10 27 0 0 0 0 0 0 0
## # ℹ 61 more rows
## # ℹ 2 more variables: start_date <date>, end_date <date>
##
## $calendar_dates
## # A tibble: 336 × 3
## service_id date exception_type
## <chr> <date> <int>
## 1 51903 2025-08-23 1
## 2 50903 2025-08-23 1
## 3 51003 2025-08-23 1
## 4 51103 2025-08-23 1
## 5 51203 2025-08-23 1
## 6 51303 2025-08-23 1
## 7 51403 2025-08-23 1
## 8 51904 2025-08-24 1
## 9 50904 2025-08-24 1
## 10 51004 2025-08-24 1
## # ℹ 326 more rows
##
## $routes
## # A tibble: 118 × 9
## route_id agency_id route_short_name route_long_name route_desc route_type
## <chr> <chr> <chr> <chr> <chr> <int>
## 1 26773 MARTA 1 Marietta Blvd/Jose… "" 3
## 2 26774 MARTA 2 Ponce de Leon Aven… "" 3
## 3 26775 MARTA 3 Martin Luther King… "" 3
## 4 26776 MARTA 4 Moreland Avenue "" 3
## 5 26777 MARTA 5 Piedmont Road / Sa… "" 3
## 6 27008 MARTA 6 Clifton Road / Emo… "" 3
## 7 26779 MARTA 8 North Druid Hills … "" 3
## 8 26780 MARTA 9 Boulevard / Tilson… "" 3
## 9 26781 MARTA 12 Howell Mill Road /… "" 3
## 10 26782 MARTA 14 14th Street / Blan… "" 3
## # ℹ 108 more rows
## # ℹ 3 more variables: route_url <chr>, route_color <chr>,
## # route_text_color <chr>
##
## $shapes
## # A tibble: 417,271 × 5
## shape_id shape_pt_lat shape_pt_lon shape_pt_sequence shape_dist_traveled
## <chr> <dbl> <dbl> <int> <dbl>
## 1 135656 33.8 -84.5 1 0
## 2 135656 33.8 -84.5 2 0.0496
## 3 135656 33.8 -84.5 3 0.101
## 4 135656 33.8 -84.5 4 0.149
## 5 135656 33.8 -84.5 5 0.172
## 6 135656 33.8 -84.5 6 0.178
## 7 135656 33.8 -84.5 7 0.184
## 8 135656 33.8 -84.5 8 0.203
## 9 135656 33.8 -84.5 9 0.257
## 10 135656 33.8 -84.5 10 0.265
## # ℹ 417,261 more rows# =========== 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"))## To install your API key for use in future sessions, run this function with `install = TRUE`.# //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.
# //TASK //////////////////////////////////////////////////////////////////////
# TASK ////////////////////////////////////////////////////////////////////////
# Convert the CRS of `census` to GCS
census <- suppressMessages(
get_acs(
geography="tract",
state="GA",
county=c("Fulton", "Dekalb", "Clayton"),
variables = c(
hhinc = "B19013_001E",
non_hispanic_white = "B03002_003E",
race_ethnic_total = "B03002_001E"
),
year = 2022,
survey = "acs5",
geometry = TRUE,
output = "wide"
)
)## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%census <- census %>%
mutate(
minority_pct = (race_ethnic_total - non_hispanic_white) / race_ethnic_total * 100,
)
# //TASK //////////////////////////////////////////////////////////////////////
# TASK ////////////////////////////////////////////////////////////////////////
# get centroids of `census` polygons,
# extract centroids that fall inside the `bbox`,
# and assign it to `home` object.
cen_central <- st_centroid(census)## Warning: st_centroid assumes attributes are constant over geometriesif (is.na(st_crs(bbox))) st_crs(bbox) <- st_crs(cen_central)
bbox <- st_transform(bbox, st_crs(cen_central))
# Keep only centroids inside the box
home <- cen_central[bbox, ]
mapview(home)# //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'
bbox_transformed <- st_transform(bbox, st_crs(osm_road$osm_line))
osm <- osm_road$osm_line[bbox_transformed,] %>% select(osm_id, highway)
net <- sfnetworks::as_sfnetwork(osm, directed = FALSE)
# clean networks
simple_net <- net %>%
activate("edges") %>%
filter(!edge_is_multiple()) %>% # remove multiple edges
filter(!edge_is_loop()) # remove loops
# Using spatial morpher
subdiv_net <- convert(simple_net, sfnetworks::to_spatial_subdivision)## Warning: to_spatial_subdivision assumes attributes are constant over geometriessubdiv_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))
smoothed_net <- convert(subdiv_net, sfnetworks::to_spatial_smooth)
# Below is done for visualization purpose
# Extract removed points for mapping purposes
removed <- setdiff(subdiv_net %>% st_as_sf('nodes') %>% pull(custom_id),
smoothed_net %>% st_as_sf('nodes') %>% pull(custom_id))
smooth_map <- leaflet() %>%
addProviderTiles(providers$CartoDB.DarkMatter) %>%
setView( -84.3854, 33.7618, zoom = 13) %>%
addPolylines(data = smoothed_net %>% st_as_sf('edges'), col = "grey", weight = 3, opacity = 0.9, popup = smoothed_net %>% st_as_sf('edges') %>% rownames()) %>%
addCircles(data = smoothed_net %>% st_as_sf('nodes'), fillColor = 'yellow', stroke = F, radius = 20, fillOpacity = 0.7, group = "kept") %>%
addCircles(data = subdiv_net %>% st_as_sf("nodes") %>% filter(custom_id %in% removed),
fillColor = "red", stroke = F, radius = 15, fillOpacity = 0.4, group = "removed") %>%
addControl(html = htmltools::HTML("<b>Red points denote removed nodes</b>"), position = "bottomright") %>%
addLayersControl(overlayGroups = c("kept", "removed"))
smooth_maposm <- smoothed_net
# **YOUR CODE HERE..**
# ...
# //TASK //////////////////////////////////////////////////////////////////////
# TASK ////////////////////////////////////////////////////////////////////////
# Add a new column named 'length' to the edges part of the object `osm`.
osm <- osm %>%
mutate(length = st_length(.))
# **YOUR CODE HERE..**
# ...
# //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
osm <- osm %>%
st_transform(crs = 26916)
home_1 <- home_1 %>%
st_transform(crs = 26916)
station <- station %>%
st_transform(crs = 26916)
paths <- st_network_paths(osm, from = home_1, to = station, type='shortest')
# //TASK //////////////////////////////////////////////////////////////////////
# TASK ////////////////////////////////////////////////////////////////////////
# Find the closest station and assign the value to `closest_station`
# Measure the route lengths
edges_sf <- st_as_sf(osm, "edges") %>%
mutate(len_m = set_units(st_length(geometry), "m"))
dist_all <- map_dbl(paths$edge_paths, function(eids) {
sum(as.numeric(edges_sf$len_m[eids]))
})
closest_idx <- which.min(dist_all)
# Get the index number for the three closest stations.
# Extract paths and stations
closest_path <- map(closest_idx, ~ osm %>%
activate("nodes") %>%
slice(paths$node_paths[[.x]]))
closest_station <- map(closest_idx, ~ station %>% slice(.x))
closest_station## [[1]]
## Simple feature collection with 1 feature and 12 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 739222.9 ymin: 3738153 xmax: 739222.9 ymax: 3738153
## Projected CRS: NAD83 / UTM zone 16N
## # A tibble: 1 × 13
## stop_id stop_code stop_name stop_desc zone_id stop_url location_type
## <chr> <chr> <chr> <chr> <chr> <chr> <int>
## 1 485 908963 ASHBY STATION 65 JOSEPH E LO… "" "" NA
## # ℹ 6 more variables: parent_station <chr>, stop_timezone <chr>,
## # wheelchair_boarding <int>, geometry <POINT [m]>, stop_lat <dbl>,
## # stop_lon <dbl># //TASK //////////////////////////////////////////////////////////////////////
# TASK ////////////////////////////////////////////////////////////////////////
# Get the distance to the closest station.
# index of the closest
closest_dist <- dist_all[closest_idx]
# //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.
m_per_mile <- 1609.344
speed_mph <- 20
trvt_osm_m <- (closest_dist / m_per_mile) / speed_mph * 60
# //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.
start_stop_name <- closest_station[[1]]$stop_name
trvt <- travel_times(filtered_stop_times = am_stop_time,
stop_name = start_stop_name,
time_range = 3600,
arrival = FALSE, # if TRUE, all journeys end at Midtown station.
max_transfers = 1,
return_coords = TRUE) %>%
filter(to_stop_name == "MIDTOWN STATION")
# //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,
date = "2025-11-10",
start_hour = 7,
end_hour = 10,
speed_mph = 20,
crs = 26916,
max_transfers = 1){
stopifnot(nrow(home) >= 1)
home_1 <- home[1, ]
osm <- st_transform(osm, crs)
home_1 <- st_transform(home_1, crs)
station <- st_transform(station, crs)
paths <- sfnetworks::st_network_paths(osm, from = home_1, to = station, type = "shortest")
# Precompute edge lengths
edges_sf <- st_as_sf(osm, "edges") %>%
mutate(len_m = set_units(st_length(geometry), "m"))
# Sum per-path lengths (meters)
dist_all <- map_dbl(paths$edge_paths, function(eids) {
sum(as.numeric(edges_sf$len_m[eids]))
})
# Closest station by network distance
closest_idx <- which.min(dist_all)
closest_station <- station[closest_idx, , drop = FALSE]
closest_dist_m <- dist_all[closest_idx]
# Home -> station time (minutes) at `speed_mph` ---
m_per_mile <- 1609.344
trvt_osm_m <- (closest_dist_m / m_per_mile) / speed_mph * 60
# Station -> Midtown via GTFS (seconds) in time window, allow ≤ max_transfers
day_use <- as.Date(date)
am_stop_time <- filter_stop_times(gtfs_obj = gtfs,
extract_date = day_use, # can't go back in time
min_departure_time = 3600*start_hour,
max_arrival_time = 3600*end_hour)
start_stop_name <- closest_station$stop_name
midtown_stop_name <- midtown$stop_name
trvt <- travel_times(filtered_stop_times = am_stop_time,
stop_name = start_stop_name,
time_range = 3600,
arrival = FALSE, # if TRUE, all journeys end at Midtown station.
max_transfers = max_transfers,
return_coords = TRUE) %>%
filter(to_stop_name == midtown_stop_name)
trvt_gtfs_m <- trvt$travel_time/60
total_trvt <- trvt_osm_m + trvt_gtfs_m
# 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)
# **YOUR CODE HERE..**
# //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
# TASK ////////////////////////////////////////
# Create an interactive map displaying `census` (polygons) and `home` (points), effectively visualizing household income and travel time to Midtown Station, respectively.
# **YOUR CODE HERE..**
census_wgs84 <- st_transform(census, 4326)
home_wgs84 <- st_transform(home, 4326)
pal_census <- colorNumeric(
palette = "YlOrRd",
domain = census_wgs84$hhinc
)
pal_home <- colorNumeric(
palette = "viridis",
domain = home_wgs84$trvt
)
m <- leaflet() %>%
addProviderTiles(providers$CartoDB.Positron, group = "Basemap") %>%
addPolygons(
data = census_wgs84,
group = "Census (Household Income)",
fillColor = ~pal_census(hhinc),
fillOpacity = 0.7,
weight = 1, # Border thickness
color = "#BDBDC3", # Light grey border
popup = ~paste("<b>Household Income:</b>", round(hhinc, 2))
) %>%
addCircleMarkers(
data = home_wgs84,
group = "Home (Travel Time)",
color = ~pal_home(trvt),
radius = 5,
stroke = FALSE, # No border on circles
fillOpacity = 0.8,
popup = ~paste("<b>Travel Time:</b>", round(trvt, 1), "min")
) %>%
addLegend(
pal = pal_census,
values = census_wgs84$hhinc,
opacity = 0.7,
title = "Household Income",
position = "bottomright",
group = "Census"
) %>%
addLegend(
pal = pal_home,
values = home_wgs84$trvt,
opacity = 0.8,
title = "Travel Time (min)",
position = "bottomleft",
group = "Home"
) %>%
addLayersControl(
overlayGroups = c("Census (Household Income)", "Home (Travel Time)"),
baseGroups = c("Basemap"),
options = layersControlOptions(collapsed = FALSE) # Keep controls open
)
m # Create an interactive map displaying `census` (polygons) and `home` (points) effectively visualizing the percentage of minority population and travel time to Midtown Station, respectively.
# **YOUR CODE HERE..**
pal_minor <- colorNumeric(
palette = "inferno",
domain = census_wgs84$minority_pct
)
y <- leaflet() %>%
addProviderTiles(providers$CartoDB.Positron, group = "Basemap") %>%
addPolygons(
data = census_wgs84,
group = "Census (Minority)",
fillColor = ~pal_minor(minority_pct),
fillOpacity = 0.7,
weight = 1, # Border thickness
color = "#BDBDC3", # Light grey border
popup = ~paste("<b>Name:</b>", NAME, "<br>",
"<b>Minority Pct:</b>", round(minority_pct, 2), "%")
) %>%
addCircleMarkers(
data = home_wgs84,
group = "Home (Travel Time)",
color = ~pal_home(trvt),
radius = 5,
stroke = FALSE, # No border on circles
fillOpacity = 0.8,
popup = ~paste("<b>Travel Time:</b>", round(trvt, 1), "min")
) %>%
addLegend(
pal = pal_minor,
values = census_wgs84$minority_pct,
opacity = 0.7,
title = "Minority Percentage",
position = "bottomright",
group = "Census"
) %>%
addLegend(
pal = pal_home,
values = home_wgs84$trvt,
opacity = 0.8,
title = "Travel Time (min)",
position = "bottomleft",
group = "Home"
) %>%
addLayersControl(
overlayGroups = c("Census (Minority %)", "Home (Travel Time)"),
baseGroups = c("Basemap"),
options = layersControlOptions(collapsed = FALSE) # Keep controls open
)
y # Create a scatter plot with a trend line showing the relationship between household income (x-axis) and travel time to Midtown Station (y-axis).
# **YOUR CODE HERE..**
ggplot(data = home, aes(x = hhinc, y = trvt)) +
geom_point(alpha = 0.6, color = "gray40") +
geom_smooth(method = "lm", formula = y ~ x, color = "blue", se = FALSE) +
labs(
title = "Income vs. Travel Time to Midtown Station",
x = "Household Income ($)",
y = "Travel Time (minutes)"
) +
theme_minimal()## Warning: Removed 6 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).
# 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).
# **YOUR CODE HERE..**
ggplot(data = home, aes(x = minority_pct, y = trvt)) +
geom_point(alpha = 0.6, color = "gray40") +
geom_smooth(method = "lm", formula = y ~ x, color = "blue", se = FALSE) +
labs(
title = "Minority Percentage vs. Travel Time to Midtown Station",
x = "Minority Percentage (%)",
y = "Travel Time (minutes)"
) +
theme_minimal()
# //TASK //////////////////////////////////////The Atlanta household income map shows income disparities across Atlanta. Northern areas show wealthier populations, whereas southern regions show lower household incomes. The percentage of minority ethnic groups also shows similar spatial patterns to the household income map.
Travel time was calculated by adding the driving time to the nearest transit and the shortest path from the nearest transit to ‘Midtown’ transit. Neighborhoods near Orme Park (Census Tract 2.02) show a longer travel time of 30 minutes than nearby neighborhoods, even though the Euclidean distance to the midtown station was actually similar to census tracts nearby. This represents the Euclidean distance is not linearly proportional to travel time, indicating the importance of measuring accessibility with travel time or travel distance, rather than using Euclidean distance.
The household income and travel time scatterplot shows a positive relationship between household income and travel time. This means that household income might be positively correlated with travel time, indicating that the shortest travel time using transit is not considered an important factor for people who live in wealthy neighborhoods. Rather, people with higher household incomes prefer to live away from transit. This seems plausible, considering the car-oriented culture of Atlanta. On the other hand, the percentage of minority ethnic groups does not show any relationship.