Bus Intersection Major Transit Stops, Spatial Pathway - HQ Transit Corridor & Major Bus Stop Explorer

import os
os.environ["CALITP_BQ_MAX_BYTES"] = str(200_000_000_000)

import intake
from calitp_data_analysis.gcs_pandas import GCSPandas
from calitp_data_analysis.gcs_geopandas import GCSGeoPandas
import calitp_data_analysis.magics

from shared_utils import webmap_utils, gtfs_utils_v2, rt_dates
import branca
import lookback_wrappers

from update_vars import analysis_date, GCS_FILE_PATH, INTERSECTION_BUFFER_METERS, MS_TRANSIT_THRESHOLD, SHARED_STOP_THRESHOLD
import datetime as dt

gcspd = GCSPandas()
gcsgp = GCSGeoPandas()

catalog = intake.open_catalog('catalog.yml')

%%capture_parameters
human_date = dt.date.fromisoformat(analysis_date).strftime('%B %d %Y (%A)')
human_date

Current Analysis Date¶

January 14 2026 (Wednesday)

If we are missing data on that date for a particular operator, we will patch in data from the previous three months. Currently patching in:

lookback_wrappers.read_published_operators(analysis_date)

{'2025-12-17': ['Bay Area 511 South San Francisco Shuttle Schedule',
  'eTrans Schedule Remix',
  'Lassen Schedule',
  'Long Beach Schedule',
  'Baldwin Park Schedule',
  'Mountain Transit Schedule',
  'Redwood Coast Schedule',
  'Lake Schedule',
  'Mendocino Schedule',
  'Redding Schedule'],
 '2025-10-15': ['Bay Area 511 Angel Island-Tiburon Ferry Schedule',
  'Nevada County Schedule',
  'El Monte Schedule'],
 '2025-11-05': ['Trinity Schedule',
  'Mariposa Grove Shuttle Schedule',
  'UCSC Schedule',
  'Amtrak San Joaquins Schedule']}

Analysis Segments and Key Stops¶

We use 1,250 meter analysis segments cut from GTFS shapes.
In each segment, we identify the stop with the highest frequency and use it to assign frequency to the segment.

hqta_segments = catalog.hqta_segments.read()

path = f'{GCS_FILE_PATH}all_bus.parquet'

max_arrivals_by_stop = gcspd.read_parquet(f"{GCS_FILE_PATH}max_arrivals_by_stop.parquet")

gdf = gcsgp.read_parquet(path)

stops = gcsgp.read_parquet(f"{GCS_FILE_PATH}stops_with_lookback.parquet")

stops = stops[['stop_id', 'stop_name', 'analysis_date',
      'schedule_gtfs_dataset_key', 'analysis_name', 'geometry']]

stops = stops.rename(columns={'geometry': 'stop_geometry'})

gdf = gdf.merge(stops, on = ['stop_id', 'schedule_gtfs_dataset_key'])
gdf = gdf[~gdf['circuitous_segment']]

map1 = gdf.copy()[['route_id', 'stop_id', 'geometry',
   'fwd_azimuth_360', 'circuitous_segment', 'hq_transit_corr',
   'ms_precursor', 'analysis_name']]

# Source - https://stackoverflow.com/a
# Posted by mkrieger1, modified by community. See post 'Timeline' for change history
# Retrieved 2025-12-08, License - CC BY-SA 4.0

azimuth_cmap = branca.colormap.LinearColormap(
        colors=list(branca.colormap.linear.viridis.colors) + list(reversed(branca.colormap.linear.viridis.colors)),
        vmin=0, vmax=360
)  # this will correctly show 0 and 360 as close together
azimuth_cmap.caption = '360-degree azimuth (heading)'

%%capture

# webmap_utils.export_legend(azimuth_cmap, 'azimuth_viridis_360a.svg', inner_labels=list(range(72, 360, 72)))

segment_state = webmap_utils.set_state_export(map1, filename = 'hqta_segments', subfolder='high_quality_transit_areas/',
                                     map_title='HQTA Segments', overwrite=True, color_col='fwd_azimuth_360',
                                     cmap = azimuth_cmap, legend_url='https://storage.googleapis.com/calitp-map-tiles/azimuth_viridis_360a.svg')

/Users/rae/workspace/cal-itp/data-analyses/_shared_utils/shared_utils/webmap_utils.py:128: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  centroid = (gdf.geometry.centroid.y.mean(), gdf.geometry.centroid.x.mean())

map2 = gdf.copy()[['stop_id', 'stop_name', 'am_max_trips_hr',
                  'pm_max_trips_hr', 'route_dir_count', 'analysis_name',
                  'stop_geometry']].set_geometry('stop_geometry')

map2['color'] = [(0,0,0)] * len(map2)

segment_stop_state = webmap_utils.set_state_export(map2, filename = 'hqta_segment_key_stops', subfolder='high_quality_transit_areas/',
                             existing_state=segment_state, map_title='Key Stops and Segments', overwrite=True,
                                                  manual_centroid=[37.336813156889704, -121.88911054161129])

webmap_utils.render_spa_link(segment_stop_state['spa_link'])

webmap_utils.display_spa_map(segment_stop_state['spa_link'])

pairs = gcspd.read_parquet(f"{GCS_FILE_PATH}pairwise.parquet")

Spatial Intersections¶

We use an azimuth (heading) approach to find intersections, segments are considered intersecting if they diverge at a 45-degree angle or greater. Our goal is to identify locations where riders have access to multiple frequent routes that can take them in different directions.
Intersections are colored in red in the map below.

intersect = gcsgp.read_parquet(f"{GCS_FILE_PATH}all_intersections.parquet")

by_segment = intersect.dissolve(['hqta_segment_id']).reset_index(drop=False)

by_segment['color'] = [(140, 0, 0)] * len(by_segment)

segment_intersect_state = webmap_utils.set_state_export(by_segment, filename = 'hqta_intersection_areas', subfolder='high_quality_transit_areas/',
                             existing_state=segment_state, map_title='Segments with Intersections', overwrite=True,
                                                       manual_centroid=[37.336813156889704, -121.88911054161129])

webmap_utils.render_spa_link(segment_intersect_state['spa_link'])

webmap_utils.display_spa_map(segment_intersect_state['spa_link'])

Intersection Buffers and Stop Groups¶

We use a 500ft buffer around the spatial intersection to find physical stops associated with the intersection.
We consider all of these physical stops to be Major Transit Stops.

by_segment.geometry = by_segment.buffer(INTERSECTION_BUFFER_METERS)

major_bus_spatial = gcsgp.read_parquet(f"{GCS_FILE_PATH}major_stop_bus.parquet")

intersect_buffered_state = webmap_utils.set_state_export(by_segment, filename = 'intersect_buffered', cache_seconds=0,
                           map_title='Intersecton Buffers', overwrite=True,
                                                         manual_centroid=[37.336813156889704, -121.88911054161129])

major_bus_spatial['color'] = [(200, 200, 255)] * len(major_bus_spatial)

intersect_major_state = webmap_utils.set_state_export(major_bus_spatial, filename = 'major_bus_spatial', cache_seconds=0,
                           existing_state=intersect_buffered_state, map_title='Buffered Intersections and Stop Groups', overwrite=True,
                                                         manual_centroid=[37.336813156889704, -121.88911054161129])

webmap_utils.render_spa_link(intersect_major_state['spa_link'])

webmap_utils.display_spa_map(intersect_major_state['spa_link'])