Road Features

import numpy as np
import geopandas as gpd

es_index = gpd.read_file('./data/demo_index/demo_index_elementary.geojson')
ms_index = gpd.read_file('./data/demo_index/demo_index_middle.geojson')
hs_index = gpd.read_file('./data/demo_index/demo_index_high.geojson')

Road Feature Analysis

The second component of this analysis is calculating a road segment-based index based on road-level variables. This index will be based on Daily Vehicle Miles Traveled (VMT) and K-Factor value, which is a congestion metric representing the fraction of average annual daily traffic (AADT) that occurs during the 30th busiest hour of the year.

Road data sourced from PennDot’s Traffic Information Repository (TIRe) provides these variables. However, this information is not consolidated in a single shapefile and is instead broken up based on both road ownership (state vs. city/local roads) and value bins (i.e. 0-5000 AADT), and was aggregated manually.

volume = gpd.read_file('./data/aadt_combined_phl.geojson')

# select only necessary variables to generate the index
volume_sel = volume[['DLY_VMT', 'K_FACTOR', 'geometry']][volume.is_valid]

# create a natural log DLY_VMT column for visualization purposes
volume_sel = volume_sel.assign(
    log_DLY_VMT = np.log1p(volume_sel['DLY_VMT'])
)

Vehicle volumes are significantly higher - often by a factor of 10 or more - on state roads in Philadelphia, namely I-95, I-676, and I-76. Taking the natural log of vehicle volume allows fluctuations across the county to be visualized.

volume_sel.explore(
    column='log_DLY_VMT',
    cmap='OrRd',
    legend=True,
    tooltip=True,
    tiles='CartoDB Positron',
    legend_kwds={'caption': 'Current K-Factor'}
)

Make this Notebook Trusted to load map: File -> Trust Notebook

volume_sel.explore(
    column='K_FACTOR',
    cmap='OrRd',
    legend=True,
    tooltip=True,
    tiles='CartoDB Positron',
    legend_kwds={'caption': 'Current K-Factor'}
)

Make this Notebook Trusted to load map: File -> Trust Notebook

VMT and K-Factor values were indexed following the same min-max normalization method. The wide range of values for VMT seen on a handful of road segments (highway roads) and the dominance of many short road segments with middling VMT values (local roads) produce a heavily left-skewed distribution that would otherwise obscure very low-volume roads that would be prime candidates for SRTS improvements. Therefore, the log-transformed VMT variable, which amplifies the difference of volume in roads with lower VMT values, will be indexed instead of raw VMT values.

def volume_index(df):
    df = df.copy()
    
    # create a normalized index value for each column where the higher the value, the lower the volume or k-factor
    df['vol_idx'] = 1 - ((df['log_DLY_VMT'] - df['log_DLY_VMT'].min()) / (df['log_DLY_VMT'].max() - df['log_DLY_VMT'].min()))
    df['kfact_idx'] = 1 - ((df['K_FACTOR'] - df['K_FACTOR'].min()) / (df['K_FACTOR'].max() - df['K_FACTOR'].min()))
    
    return df

volume_idx = volume_index(volume_sel)

Demographic variable values will eventually be joined to each road segment, though this is complicated by how not every road segment is fully contained within each school catchment zone. To circumvent this, the catchment zone that each road segment primarily intersects with - based on the length of each road segment that falls within each catchment area - will be associated with the latter.

# define function to determine which catchment area each road segment predominantly fall within
def assign_catchment_id(roads_gdf, catchment_gdf, id_col):
    
    # Create copy to avoid modifying original
    roads = roads_gdf.copy()
    roads['road_idx'] = roads.index
    
    # Perform spatial intersection
    intersected = gpd.overlay(roads, catchment_gdf[[id_col, 'geometry']], how='intersection')
    
    # Calculate intersection length
    intersected['intersection_length'] = intersected.geometry.length
    
    # Find catchment zone with maximum intersection length for each road segment
    idx_max = intersected.groupby('road_idx')['intersection_length'].idxmax()
    max_intersections = intersected.loc[idx_max]
    
    # Merge the catchment ID back to original roads
    roads = roads.merge(
        max_intersections[['road_idx', id_col]], 
        on='road_idx', 
        how='left'
    )
    
    # Drop temporary index column
    roads = roads.drop(columns='road_idx')
    
    return roads

# assign catchment IDs to volume index values and remove road segments that do not intersect any catchment area
es_volume_idx = assign_catchment_id(volume_idx, es_index, 'ES_ID').dropna(subset=['ES_ID'])
ms_volume_idx = assign_catchment_id(volume_idx, ms_index, 'MS_ID').dropna(subset=['MS_ID'])
hs_volume_idx = assign_catchment_id(volume_idx, hs_index, 'HS_ID').dropna(subset=['HS_ID'])

The final index generated from this analysis relies on how each is independently normalized through the same method, so that they can then be added together to form a final index value. No weighting of these values was done; however, future iterations of this approach could implement weights on specific variables to scale their impact on this index value (e.x. multiplying normalized VMT values )

# join demographic index values to volume index values based on catchment area IDs
es_final = es_volume_idx.merge(es_index[['ES_ID', 'demo_index']], on='ES_ID', how='left')
ms_final = ms_volume_idx.merge(ms_index[['MS_ID', 'demo_index']], on='MS_ID', how='left')
hs_final = hs_volume_idx.merge(hs_index[['HS_ID', 'demo_index']], on='HS_ID', how='left')

# calculate a combined index value based on demographic and volume indices
def combined_index(df):
    df = df.copy()
    df['combined_index'] = df['demo_index'] + df['vol_idx'] + df['kfact_idx']
    return df

es_final_idx = combined_index(es_final)
ms_final_idx = combined_index(ms_final)
hs_final_idx = combined_index(hs_final)