GTFS Data

gtfs_cta <- read_gtfs("https://www.transitchicago.com/downloads/sch_data/google_transit.zip")

The first and most important data for us to retrieve is the GTFS data, which includes the subway stations, routes, and their geographical locations. There is publicly available GTFS data that we can directly retrieve from CTA’s official website.

After downloading the data, what we had to do was to filter out the subway routes (route_type == 1) & (location_type == 1). We joined the data from the stops, routes, trips and shapes columns to get the shape sf objects for the routes and the stations.

Here is the map that shows the CTA “L” Train network that we retrieved from the GTFS data. We can see that all the subway lines starts from a loop in the downtown region and goes radial into the suburban area. The subway loop is the namesake of Chicago’s downtown or CBD region, “the Loop”.

Ridership Data

The official ridership data can be found at CTA’s website as well as the Chicago Data Portal. It includes daily ridership for every station dating back from 2001 til this date (2024). We loaded the dataset in R and did data cleaning procedures. We wanted to gfilter out time range between the years 2015 and 2019, as these are the last years before subway ridership is affected by the COVID event.

We also filtered out only the weekday riderships of each station, as we think it may be most representative of the daily commuting patterns of the subway. Then we grouped then ridership by station and calculated the daily means of each station.

# Use forward slashes
ridership <- read.csv("CTA_-_Ridership_-__L__Station_Entries_-_Daily_Totals_20240922.csv")
ridership <- ridership %>%
  mutate(station_id = as.character(station_id)) %>%
  mutate(date = trimws(date), date = mdy(date)) %>%
  filter(date >= as.Date("2015-01-01") & date <= as.Date("2019-12-31"))

avg_rides <- ridership %>% 
  filter(daytype == "W") %>%
  mutate(year = year(date)) %>% 
  group_by(station_id, stop_name) %>%
  #group_by(station_id, stop_name, year) %>% 
  summarise(avg_rides = mean(rides, na.rm = TRUE), .groups = 'drop') 

stations_ridership <- station_sf %>%
  left_join(avg_rides, by = c("stop_id" = "station_id"))

Here is the map showing the riderships of each station, with larger circle indicating higher ridership. Unfortunately the size-based legend is not available, but you can view ech station’s ridership by clicking it and it’ll show in the pop-up.

Census Demographic Data

We also can fetch the data from the Census Bureau, which is quite convenient with a census api key using the tidycensus package.

The data extraction uses the get_acs function to fetch block group-level data from the American Community Survey (ACS) for Cook and Lake Counties in Illinois. Important variables include household income, population, vehicle ownership, and senior population by age and gender. These variables are processed to find total counts (like total senior population) and densities (like car density, senior density, and population density) by dividing counts by the area of each spatial unit. The data is transformed to the same coordinate system and filtered to remove empty geometries. This helps combine demographic and spatial details for analysis.

Bike Stations Data

There might be a pattern for subway passengers to use bike to reach their destinations within a mile from subway stations. So we used GBFS data from Divvy and calculated the density of bike stations for each block group.

OSM Street Network and Network Analysis on Service Area

We would like to fetch the road network data for the Chicago region and the area served by the “L” Train network. We did this in the osmdata package and created a custom bounding box to download the data. Ideally we should do all analysis in r markdown, but due to the large size of road network data(~800MB), we had to use ArcGIS for the network analysis.

Here is the bounding box we created to download the road network.

Here is the code we used to fetch the road network data.

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

After we downloaded the road network, we converted it to an sf network object, and activated the edges. To clean out the road network, we first filtered out all edges which are same pairs ( filter(!edge_is_multiple()) ) and edges which starts and ends at the same node ( filter(!edge_is_loop()) ). There are some road intersections that are not marked as nodes in this road network, and we used sfnetworks::to_spatial_subdivision to create nodes for those intersections. There are also “pseudo” nodes which only connects to two edges which is not essentially a network node so we deleted it using sfnetworks::to_spatial_smooth.

From this point on, the capacity in R simple features class does not allow us to do network analysis on the massive network. We exported the sfnetwork to a shapefile and did the analysis in ArcGIS Pro.

In ArcGIS, we imported the osm road network and the CTA Station point locations. In the osm_edges shapefile, we created two new columns, one is “length” that calculates the geographical length of all road sections, and the other is “minutes” that calculates the estimated walking time (in minutes) with 5 km per hour speed using the “length” column.

We built a Network Dataset using the osm_edges and osm_nodes shapefiles, and set the traveling mode to minutes using the column we created earlier. We are able to run a 5-10-15 minute walking distance service area analysis. We exported the result to a shapefile and imported it back to R.

## Reading layer `cta_polygons' from data source 
##   `C:\Users\benso\Documents\Georgia Tech\CP8883\Final Project\shp\cta_polygons.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -87.90694 ymin: 41.71177 xmax: -87.59051 ymax: 42.08323
## Geodetic CRS:  WGS 84

We decided to use the 15-minute walking distance region. Here’s what it looks like:

We assigned each census block a binary value based on whether it intersects with the 15-min walking distance polygon.

Til this point, all the data we needed are ready for us to do a regression analysis on the different variables.

Distribution of Car Density

The distribution of car density is highly skewed, with most areas having very low vehicle density and a few areas exhibiting extremely high densities, forming a long tail. This suggests vehicles are concentrated in certain urban regions, while many areas may be less car-reliant or less populated. The skewness reflects varied land-use patterns, ranging from dense urban cores to sparsely populated suburban zones.

Car Density vs Intersects

The scatterplot of car density vs. intersects shows two clear clusters due to the binary nature of intersects: areas with CTA coverage (1) and without coverage (0). The red logistic regression curve indicates a positive relationship, where higher car density increases the likelihood of CTA coverage, especially in moderately dense areas. However, the curve flattens at very high car densities, suggesting that beyond a certain threshold, more cars do not significantly affect the probability of coverage. This reveals that CTA focuses on areas with moderate car density, balancing transit accessibility with urban vehicle use.

Distribution of Household Income

The distribution of household income has a single peak around the $50,000–$100,000 range, representing median income levels in the area. The curve steadily declines as income increases, with fewer households earning above $150,000. This distribution highlights the region’s socioeconomic diversity, with many middle-income households and fewer at the extremes. Most blocks are neither extremely wealthy nor impoverished, reflecting common urban and suburban trends.

Household Income vs Intersects

The scatter plot of household income vs. CTA station coverage shows a near-flat regression line, indicating little correlation between income and proximity to a CTA station. However, the points reveal that higher-income households are more often found in areas without CTA coverage (intersects = 0). This suggests wealthier neighborhoods may depend less on public transit, pointing to a possible disparity in transit usage across income levels.

Distribution of Population Density

The distribution of population density is highly skewed, with most areas having densities below 10,000 people per square kilometer. The number of areas drops sharply as density increases, indicating that very dense blocks are uncommon and concentrated in specific regions. This pattern reflects urban planning, where high densities are usually found near city centers or transit hubs.

Population Density vs. Intersects

The scatterplot of population density vs. CTA coverage (intersects) reveals a non-linear relationship. The red curve indicates that higher population densities increase the likelihood of CTA coverage, but the effect levels off beyond a certain density. This suggests the CTA focuses on moderately dense urban areas rather than exclusively targeting the most densely populated blocks.

Distribution of Bike Density

The distribution of bike density shows a heavily skewed pattern, with most areas reporting very low values close to zero. Only a small number of areas have higher bike density, indicating that bike infrastructure is concentrated in specific locations, leaving much of the region with minimal resources.

Bike Density vs Intersects

The scatterplot of bike density vs. CTA coverage (intersects) shows an upward trend, with higher bike density linked to a greater likelihood of being within 15 minutes of a CTA station. However, distinct patterns emerge:

1. High Bike Density Distribution:

• Unlike previous variables (e.g., car density or population density), where areas with high values tend to cluster around intersects = 1 (covered by CTA), here we observe a more balanced distribution of high bike density points across both intersects = 0 (not covered) and intersects = 1 (covered).
• This suggests that areas with higher bike density are not exclusively reliant on CTA coverage, indicating that bike infrastructure may serve as a complementary transit option, independent of CTA presence.

2. Low Bike Density Distribution:

• For areas with low bike density, the points are more evenly distributed between intersects = 0 and intersects = 1, reinforcing the idea that bike density alone does not determine CTA coverage in these regions.
• This highlights the potential role of other factors, such as population density or urban planning priorities.

3. Unique Insights:

• The balanced presence of high bike density points across both intersects = 0 and intersects = 1 suggests that bike infrastructure might play an independent role in transit accessibility.
• This differs from variables like car or population density, where high values more strongly align with CTA coverage, indicating that bike density is less predictive of CTA coverage in isolation.
• This result points to the need for a combined transit strategy, where bikes complement rather than compete with CTA services, particularly in areas with well-established bike infrastructure.