import multiprocessing as mp
import time
import pandas as pd
import geopandas as gpd
import osmnx as ox
import networkx as nx
import shapely
def fmt_geoid( in_str ):
# Last digits of GEOD (remove the US prefix)
return in_str[9:]
trip_times = 600 # 10 mins
bus_stop_buffer = 100 # meters
"""
Premise:
We want to determine how accessible public transit is in Columbus OH/Franklin County vs. Indianapolis IN/Marion County. This
will be achieved by creating a model to calculate how many census blocks are within "reasonable walkability" to a bus stop
(X minute walk at Y kph) or any block within Z meters of a bus stop.
"""
class City:
# WIP class for studying transit accessibility in major metro areas. This serves as the roadmap on how the
# transit analysis was done.
# In order to not pass the entire graph to the multiprocessing pool, we can establish it here to work as
# a sort of global variable.
worker_graph = None
# These work like global variables used as a folium styling function.
covered_series = None
stop_on_loc_series = None
def __init__( self, name, utm_crs ):
# When creating a city, input which UTM zone it exists in so that all data
# added later is projected to this CRS.
self.name = name # string
self.crs = utm_crs # pyproj CRS
self.coverage_gdf = None
self.units_per_block = None
self.stops = None
self.stop_network = None
self.analysis = None
self.walk_graph = None
self.walk_nodes = None
self.walk_edges = None
def read_coverage_file( self, fname ):
# When a transit coverage polygon is created for the metro area, add it here.
in_gdf = gpd.read_file( fname )
self.coverage_gdf = in_gdf.set_crs( self.crs, allow_override=True )
# Create unique identifier for join operation in analysis. Used to style webmap in Folium
self.coverage_gdf['covered'] = True
print( f'{self.name} transit coverage polygon loaded' )
# In the future, I will re-engineer the creation of the coverage polygon as a method to this class.
# def create_coverage_file( ... )
def read_census_blocks( self, hu_fname, blocks_fname ):
# Inputs, file names for the housing units data (from US census) and blocks geodata.
hu_block = pd.read_csv( hu_fname ).drop(0) # First row is like a column metadata we don't want
hu_col = 'H1_001N'
# Now run this on the whole column.
hu_block[ 'GEOID_FIX' ] = hu_block['GEO_ID'].apply( fmt_geoid )
# Convert dtype of the housing units to int
hu_block[ hu_col ] = hu_block[ hu_col ].astype( int )
gdf = gpd.read_file( blocks_fname )
# Now merge spatial data with housing unit count data
self.units_per_block = gdf.merge(
hu_block,
left_on='GEOID',
right_on='GEOID_FIX',
how='inner'
)
self.units_per_block.to_crs( self.crs, inplace=True )
print( f'{self.name} census blocks loaded' )
def read_metro_stops( self, stops_file ):
# Input is any file format that can be read using geopandas, that contains point bus stop data. Any input
# crs is supported because it is immediately projected.
self.stops = gpd.read_file( stops_file )
self.stops.to_crs( self.crs, inplace=True )
print( f'{self.name} stops loaded' )
def analyze_transit_accessibility( self, bus_stop_distance_threshold=100 ):
# Creates city.analysis instance variable, which is based off the city.units_per_block dataframe.
# Adds true/false columns for the 2 "coverage" criteria, e.g. has a bus stop within
# <bus_stop_distance_threshold> meters or is within an n minute walk on
# the walking network.
p = 'within'
p2 = 'intersects'
if self.stops is None:
print( 'bus stops have not been loaded.' )
return
if self.units_per_block is None:
print( 'housing units data has not been loaded.' )
return
# Determine whether the census block is within a simple 100m buffer of each bus stop.
stop_network = shapely.ops.unary_union( self.stops['geometry'].tolist() )
stop_network = stop_network.buffer( bus_stop_distance_threshold )
self.stop_network = gpd.GeoDataFrame( geometry=[stop_network], crs=self.crs )
# Create unique identifier for join operation in analysis. Used to style webmap in Folium
self.stop_network['bus_stop_onsite'] = True
# Perform a simple intersect. We will assume that a person can reasonably walk to the
# end of their block to begin the route to their bus stop. This is much more lenient than
# a typical transit analysis would be, because most people will prefer to take their car
# unless transit is undeniably more convenient.
analysis = self.units_per_block.sjoin( self.stop_network, how='left', predicate=p2 )
# Remove unnecessary columns (including the index_right columns that result in error
# from later spatial joins.
cols = ['BLOCK', 'GEOID', 'H1_001N', 'geometry', 'bus_stop_onsite']
analysis = analysis[ cols ]
# Determine whether census block is within the walkability network of bus stops.
analysis = analysis.sjoin( self.coverage_gdf, how='left', predicate=p )
# Remove unnecessary columns (including the index_right columns that result in error
# from later spatial joins.
cols = ['BLOCK', 'GEOID', 'H1_001N', 'geometry', 'covered', 'bus_stop_onsite']
analysis = analysis[ cols ]
analysis['covered'] = analysis['covered'].fillna( False )
analysis['bus_stop_onsite'] = analysis['bus_stop_onsite'].fillna( False )
print( analysis.head() )
self.analysis = analysis
def set_series_vars( self ):
# For folium plotting, set the series global variables to the dataframe from the city class. Used due to
# limitations with the style_function parameter in the folium.GeoJson function.
global covered_series
global stop_on_loc_series
covered_series = self.analysis.set_index( 'GEOID' )['covered']
stop_on_loc_series = self.analysis.set_index( 'GEOID' )['bus_stop_onsite']
def build_walking_network(self, place, walking_speed_kph=5):
#download and prepare the pedestrian street network.
start = time.time()
if place is None:
place = self.name
print(f"Downloading walking network for {place}...")
G = ox.graph_from_place( place, network_type="walk" )
G = ox.project_graph(G)
nx.set_edge_attributes( G, values=walking_speed_kph, name="speed_kph" )
G = ox.add_edge_travel_times(G)
self.walk_graph = G
self.walk_nodes, self.walk_edges = ox.graph_to_gdfs(G)
# Match every layer to the graph CRS
self.crs = self.walk_nodes.crs
if self.stops is not None:
self.stops = self.stops.to_crs(self.crs)
if self.units_per_block is not None:
self.units_per_block = ( self.units_per_block.to_crs(self.crs) )
end = time.time()
runtime = end-start
print(
f"Walking network loaded "
f"({len(self.walk_nodes)} nodes, "
f"{len(self.walk_edges)} edges) "
f"runtime: {runtime} s"
)
# Static method allows for class-related workflows that do not rely on the
# instance of the class
@staticmethod
def _make_iso_polygon(args):
#Worker function for multiprocessing.
graph = city._worker_graph
center_node, trip_time, edge_buffer = args
travel_time = nx.single_source_dijkstra_path_length( graph,
center_node, weight="travel_time" )
subnodes = [ n for n, t in travel_time.items() if t <= trip_time ]
subgraph = graph.subgraph( subnodes )
if subgraph.number_of_edges() == 0:
return None
_, subedges = ox.graph_to_gdfs( subgraph )
return subedges.buffer( edge_buffer ).unary_union
def build_coverage_polygon( self, trip_time=600, edge_buffer=80, processes=4 ):
# Generate a transit coverage polygon by constructing a walking isochrone around
# every transit stop.
if self.walk_graph is None:
print( "Walking network has not been built." )
return
if self.stops is None:
print( "Bus stops have not been loaded.")
return
print("Finding nearest network nodes...")
# Find nearest nodes from OSMnx network over each bus stop
stop_nodes = ox.distance.nearest_nodes( self.walk_graph,
X=self.stops.geometry.x, Y=self.stops.geometry.y )
# For benchmark analysis, uncomment
#stop_nodes = stop_nodes[:10]
self.stop_nodes = stop_nodes
print( f"Generating {len(stop_nodes)} walking isochrones...")
args = [( node, trip_time, edge_buffer ) for node in stop_nodes ]
city._worker_graph = self.walk_graph
start = time.time()
with mp.Pool(processes) as pool:
polygons = pool.map( city._make_iso_polygon, args )
end = time.time()
runtime = end-start
print( f'Isochrone runtime: {runtime}' )
polygons = [ p for p in polygons if p is not None ]
print( f"Merging {len(polygons)} polygons..." )
coverage_polygon = shapely.ops.unary_union( polygons )
self.coverage_polygon = coverage_polygon
self.coverage_gdf = gpd.GeoDataFrame({"covered": [True]},
geometry=[coverage_polygon],
crs=self.crs,
)
print("Coverage polygon created.")
def save_coverage(self, filename):
#Save the generated transit coverage polygon.
if self.coverage_gdf is None:
print( "Coverage polygon has not been generated." )
return
self.coverage_gdf.to_file( filename, driver="GeoJSON", )
print(f"Coverage saved to '{filename}'")
def make_styling_function( self ):
# Create a function to use in folium when creating the web app. This is because
# the styles are being determined by the dataframe indexes.
covered_series = self.analysis.set_index( 'GEOID' )['covered']
stop_on_loc_series = self.analysis.set_index( 'GEOID' )['bus_stop_onsite']
def styling_function( feature ):
covered = covered_series.get(int(feature["id"][-5:]), None)
stop_on_loc = stop_on_loc_series.get(int(feature["id"][-5:]), None)
return {
"fillOpacity": 0.25,
"weight": 0.75,
"fillColor": "#00ff00" if covered or stop_on_loc else "#ff0000"
}
return styling_function
def assessment( self ):
df = self.analysis
total_hu_count = df['H1_001N'].sum()
hu_covered = df.loc[ df['covered'] | df['bus_stop_onsite'] ]['H1_001N'].sum()
p = hu_covered / total_hu_count
ds = [{
'Total Housing Units': total_hu_count,
'Housing Units in Covered Areas:': hu_covered,
'Proportion': p
}]
outdf = pd.DataFrame( ds )
print( outdf.head() )
self.assessment = outdf
1) Download H1 dataset with H1_001N (# of housing units per census block) for the metro area.
Link: https://data.census.gov/table/DECENNIALSF12000.H001?q=Housing+Units&g=050XX00US18097$1000000
** Ensure the filter is set for the right area before downloading.
2) Download the US Census Blocks geographic data ** Link: https://hub.arcgis.com/datasets/d795eaa6ee7a40bdb2efeb2d001bf823_0/
3) Obtain transit data for the metro, ensure there is latitude/longitude data so that it can be used by OSMnx.
4) Update the 'cfg' variable with a new dict entry, and set each variable to the file name you used. If any directories were added, ensure they are specified as well.
5) Use the methods to create and optionally save a "transit coverage area" for the metro area.
6) Use the .analyze_transit_accessiblity() method to generate the city.analysis dataframe, so that it can be visualized.
# First start by configuring each city. Add your own cities as entries to the configuration!
# Possible extension: create programmatric structures to organize city data by state.
# Configure the input datasets for each city.
cfg = {
'Marion County, Indiana, USA':
{
'CRS':'EPSG:26916',
'housing_units_file': r'Data/Indianapolis/IN_DECENNIALDHC2020.H1-Data.csv',
'census_block_geodata': r'Data/Indianapolis/US_Census_Blocks_Marion_County_IN.geojson',
'metro_geodata': r'Data/Indianapolis/IndyGo_Bus_Stops.geojson',
'coverage_file': r'Data/Indianapolis/indy_covg_area.geojson'
},
'Franklin County, Ohio, USA':
{
'CRS':'EPSG:32617',
'housing_units_file': r'Data/Columbus/OH_DECENNIALDHC2020.H1-Data.csv',
'census_block_geodata': r'Data/Columbus/US_Census_Blocks_Franklin_County_OH.geojson',
'metro_geodata': r'Data/Columbus/COTA_Stops_Current_full.geojson',
'coverage_file': r'Data/Columbus/cbus_covg_area.geojson'
}
}
WARNING: For large metro areas (Columbus, Indianapolis) this can take HOURS
Careful! It is easy to overwrite the 'datasets' variable
#datasets = {}
for name, config in cfg.items():
print( f'**{name}**' )
# Initialize and establish the metro area to study.
#datasets[ name ] = city( name, config['CRS'] )
# Download the walking network from OSMnx
datasets[ name ].build_walking_network( name )
datasets[ name ].read_metro_stops( config['metro_geodata'] )
# Create the coverage polygon from the walking network.
datasets[ name ].build_coverage_polygon( processes=6 )
# Use this line to save the coverage dataset as needed.
#datasets[ name ].save_coverage( 'file.geojson' )
Visualize Inputs (merged)
datasets['Marion County, Indiana, USA'].coverage_gdf.geometry[0]
datasets['Franklin County, Ohio, USA'].coverage_gdf.geometry[0]
datasets = {}
for name, config in cfg.items():
print( f'**{name}**' )
# Initialize and establish the metro area to study.
datasets[ name ] = City( name, config['CRS'] )
# Read in the H1 dataset for the metro area from data.census.gov (2020 census must be used), as well as the
# geographic data for census blocks
datasets[ name ].read_census_blocks( config['housing_units_file'], config['census_block_geodata'] )
# Assuming you have already created a transit coverage area polygon dataset using the .build_coverage_polygon() and
# saved it using the .save_coverage() methods, you can read it in using this method to save hours of re-processing it
# each time.
datasets[ name ].read_coverage_file( config['coverage_file'] )
# GTFS transit feed data is preferred here.
datasets[ name ].read_metro_stops( config['metro_geodata'] )
# Performs the analysis, saves as a dataframe that can be accessed using city.analysis.
datasets[ name ].analyze_transit_accessibility( bus_stop_buffer )
# Display assessment (get simple proportion of housing units in service area)
for obj in datasets.values():
obj.assessment()
import folium
m = folium.Map(location=[39.75, -84.5], zoom_start=8) # Centered between Cbus and Indy on start
radius = 4
# Add each city to web map:
for name, obj in datasets.items():
# Adds the census blocks, color coded by transit accessibility to web map.
print( f'**{name}**' )
obj.set_series_vars()
folium.GeoJson( obj.analysis, name=f'{name} Census Blocks', style_function=obj.make_styling_function() ).add_to( m )
# Opportunity for advancement... ensure description columns are the same in all stop sets
if name == 'Marion County, Indiana, USA':
desc_col = 'DESCRIPTION'
else:
desc_col = 'StopName'
for _, row in obj.stops.to_crs( 'EPSG:4326' ).iterrows():
# Adds each bus stop from the transit service
lon = row.geometry.x
lat = row.geometry.y
folium.CircleMarker(
location=[ lat, lon ],
radius=radius,
color="black",
weight=1,
fill_opacity=0.6,
opacity=1,
fill_color="green",
fill=False, # gets overridden by fill_color
tooltip=row[ desc_col ],
).add_to(m)
folium.LayerControl().add_to(m)
# Run to save map, open in browser.
# The notebook wouldn't save with the app running inside of it, so this was my workaround.
# I had to refresh the page a few times for the HTML to work in my browser, but it does work.
m.save( 'indy_cbus_transit.html' )
Kåresdotter, E., Page, J., & Håkansson, M. (2022). First mile/last mile problems in smart and sustainable cities: A case study in Stockholm County. Journal of Urban Technology, 29(2), 115–137. https://doi.org/10.1080/10630732.2022.2033949
Sarker, R. I., Mailer, M., & Sikder, S. K. (2020). Walking to a public transport station: Empirical evidence on willingness and acceptance in Munich, Germany. Smart and Sustainable Built Environment, 9(1), 38–53. https://doi.org/10.1108/SASBE-07-2017-0031
Lu, Y., Kimpton, A., Prato, C. G., Sipe, N., & Corcoran, J. (2024). First and last mile travel mode choice: A systematic review of the empirical literature. International Journal of Sustainable Transportation, 18(1), 1–14. https://doi.org/10.1080/15568318.2023.2218285
Central Ohio Transit Authority. (n.d.). Bus stops. https://cota1974.maps.arcgis.com/home/item.html?id=c7ce9e3057894e56902a9853bf4c96bb. Retrieved July 28, 2026.
City of Indianapolis & Marion County. (2025). IndyGo bus stops. https://data.indy.gov/datasets/indygo-bus-stops/about. Retrieved July 2, 2026.
U.S. Census Bureau. (2020). 2020 Census Demographic and Housing Characteristics File: Table H1—Housing units. Retrieved July 31, 2026, from https://data.census.gov/table/DECENNIALDHC2020.H1?q=H1:+HOUSING+UNITS&g=050XX00US39049$1000000
U.S. Census Bureau. (2020). 2020 Census Demographic and Housing Characteristics File: Table H1—Housing units. Retrieved July 31, 2026, from https://data.census.gov/table/DECENNIALDHC2020.H1?q=H1:+HOUSING+UNITS&g=050XX00US18097$1000000
U.S. Census Bureau. (2024). American Community Survey 1-year estimates: Table B08134—Means of transportation to work by travel time to work. Retrieved August 3, 2026, from https://data.census.gov/table/ACSDT1Y2024.B08134?q=B08134&g=040XX00US37
U.S. Census Bureau. (2020). U.S. Census Blocks. Federal Geographic Data Committee / ArcGIS Hub. Retrieved July 31, 2026, from https://hub.arcgis.com/datasets/fedmaps::u-s-census-blocks-1/about