GTFS to Public Transit Graphs¶

# Load GTFS data into DuckDB (current API)
sample_gtfs_path = Path("./data/itm_london_gtfs.zip")

print("Loading London Transport GTFS data into DuckDB...")
gtfs_con = c2g.load_gtfs(sample_gtfs_path)
print("GTFS database ready.")

# Inspect loaded GTFS tables from DuckDB
tables_df = gtfs_con.execute("SHOW TABLES").df()
table_names = sorted(tables_df["name"].tolist())

print(f"Found {len(table_names)} data tables")
print(", ".join(table_names))

for table in ["stops", "routes", "stop_times"]:
    if table in table_names:
        n = gtfs_con.execute(f"SELECT COUNT(*) AS n FROM {table}").df().loc[0, 'n']
        print(f"Total {table}: {n:,}")

# Explore GTFS tables available in DuckDB
print("Available GTFS tables:")
for i, table_name in enumerate(table_names, 1):
    num_records = gtfs_con.execute(
        f"SELECT COUNT(*) AS n FROM {table_name}"
    ).df().loc[0, 'n']
    print(f"  {i}. {table_name}: {num_records:,} records")

print("\n" + "=" * 50)
print("GTFS Table Descriptions:")
print("=" * 50)
print("agency         - Transit operators (TfL, etc.)")
print("calendar       - Service patterns (weekdays/weekends)")
print("routes         - Transit lines")
print("stops          - Physical stop locations (with coordinates)")
print("stop_times     - Scheduled arrivals/departures")
print("trips          - Individual vehicle journeys")
print("calendar_dates - Service exceptions (holidays, etc.)")

table_names

# Agency information - Who operates the transit services?
print("Transit Agencies:")
print("This table contains information about transportation operators")
gtfs_con.execute("SELECT * FROM agency LIMIT 5").df()

gtfs_con.execute("SELECT * FROM calendar LIMIT 5").df()

gtfs_con.execute("SELECT * FROM calendar_dates LIMIT 5").df()

gtfs_con.execute("SELECT * FROM routes LIMIT 5").df()

# Stops - The spatial foundation of transit networks
stops_df = gtfs_con.execute(
    """
    SELECT * EXCLUDE (geometry), ST_AsText(geometry) AS geometry_wkt
    FROM stops
    LIMIT 5
    """
).df()
stops_preview = gpd.GeoDataFrame(
    stops_df.drop(columns=["geometry_wkt"]),
    geometry=gpd.GeoSeries.from_wkt(stops_df["geometry_wkt"]),
    crs="EPSG:4326",
)

print("Transit Stops (with Spatial Coordinates):")
print("Notice how city2graph stores stop geometry in DuckDB and we can convert to GeoPandas")
print(f"Coordinate Reference System: {stops_preview.crs}")
print(f"Geometry type: {stops_preview.geometry.geom_type.iloc[0]}")

stops_preview.head()

gtfs_con.execute("SELECT * FROM stop_times LIMIT 5").df()

gtfs_con.execute("SELECT * FROM trips LIMIT 5").df()

# Build a full stops GeoDataFrame from DuckDB for quick map preview
stops_all_df = gtfs_con.execute(
    """
    SELECT * EXCLUDE (geometry), ST_AsText(geometry) AS geometry_wkt
    FROM stops
    """
).df()
stops_gdf = gpd.GeoDataFrame(
    stops_all_df.drop(columns=["geometry_wkt"]),
    geometry=gpd.GeoSeries.from_wkt(stops_all_df["geometry_wkt"]),
    crs="EPSG:4326",
)

# Reproject to British National Grid for accurate map display
stops_gdf_bng = stops_gdf.to_crs(epsg=27700)

print("Visualizing London's Transit Stops")
ax = c2g.plot_graph(nodes=stops_gdf_bng)
ctx.add_basemap(ax, crs=stops_gdf_bng.crs, source=ctx.providers.CartoDB.DarkMatter)
plt.show()

# Transform GTFS schedules into a travel network graph
travel_summary_nodes, travel_summary_edges = c2g.travel_summary_graph(
    gtfs_con,
    calendar_start="20250601",  # Analyze services for June 1, 2025
    calendar_end="20250601",     # Single day analysis for demonstration
)

travel_summary_nodes = travel_summary_nodes.to_crs(epsg=27700)
travel_summary_edges = travel_summary_edges.to_crs(epsg=27700)

print("Created graph with:")
print(f"   • {len(travel_summary_nodes):,} nodes (stops with connections)")
print(f"   • {len(travel_summary_edges):,} edges (stop-to-stop connections)")
print("   • Each edge contains travel time (seconds) and frequency data")

travel_summary_nodes.head()

# Examine the edges (connections) in our travel network
print("Network Edges (Transit Connections):")
print("Each row represents a direct connection between two stops")
print("\nKey metrics city2graph calculated:")
print("   • travel_time_sec: Average time to travel between stops (seconds)")
print("   • frequency: Number of services per day on this connection")
print("   • geometry: LineString for mapping and spatial analysis")

print(f"\nPerformance insight:")
print(f"   Fastest connection: {travel_summary_edges['travel_time_sec'].min():.0f} seconds")
print(f"   Busiest connection: {travel_summary_edges['frequency'].max():.0f} services/day")
print(f"   Average travel time: {travel_summary_edges['travel_time_sec'].mean():.0f} seconds")

travel_summary_edges.head()

# Examine the edges (connections) in our travel network
print("Network Edges (Transit Connections):")
print("Each row represents a direct connection between two stops")
print("\nKey metrics city2graph calculated:")
print("   • travel_time_sec: Average time to travel between stops (seconds)")
print("   • frequency: Number of services per day on this connection")
print("   • geometry: LineString for mapping and spatial analysis")

print(f"\nPerformance insight:")
print(f"   Fastest connection: {travel_summary_edges['travel_time_sec'].min():.0f} seconds")
print(f"   Busiest connection: {travel_summary_edges['frequency'].max():.0f} services/day")
print(f"   Average travel time: {travel_summary_edges['travel_time_sec'].mean():.0f} seconds")

travel_summary_edges.head()

# Get Greater London boundary using city2graph
london_boundary = c2g.get_boundaries("Greater London, UK").to_crs(epsg=27700)

# Ensure projected CRS for spatial operations
travel_summary_nodes = travel_summary_nodes.to_crs(epsg=27700)
travel_summary_edges = travel_summary_edges.to_crs(epsg=27700)

# Clip graph to Greater London extent
print("Clipping nodes and edges to Greater London boundary...")
travel_summary_nodes, travel_summary_edges = c2g.clip_graph(
    (travel_summary_nodes, travel_summary_edges),
    london_boundary,
    keep_outer_neighbors=False,
 )

print("Boundary clipping complete:")
print(f"   Nodes within London: {len(travel_summary_nodes):,}")
print(f"   Edges within London: {len(travel_summary_edges):,}")

# Quick check of London-wide network
print("Greater London network summary")
print(f"Nodes: {len(travel_summary_nodes):,}")
print(f"Edges: {len(travel_summary_edges):,}")
travel_summary_nodes.head()

# Plot the network with dual encoding (color + width) using city2graph
c2g.plot_graph(
    nodes=travel_summary_nodes,
    edges=travel_summary_edges,
    edge_color='travel_time_sec',
    edge_linewidth=travel_summary_edges['frequency'] / 250,
    edge_alpha=0.8,
    )
plt.show()

travel_summary_graph = c2g.gdf_to_nx(
    travel_summary_nodes,
    travel_summary_edges,
)

travel_rx_graph = c2g.nx_to_rx(travel_summary_graph)

betweenness_centrality = rx.betweenness_centrality(travel_rx_graph)

# Set the betweenness centrality as a node attribute
nx.set_node_attributes(
    travel_summary_graph,
    betweenness_centrality, 
    "betweenness_centrality"
    )

travel_nodes, travel_edges = c2g.nx_to_gdf(travel_summary_graph)

# Calculate quantiles for node size to handle skew in centrality values
# Drop duplicate bin edges to avoid qcut errors when many nodes share the same centrality
centrality_bins = pd.qcut(
    travel_nodes['betweenness_centrality'],
    q=10,
    labels=False,
    duplicates="drop"
)

travel_nodes['centrality_quantile'] = centrality_bins.fillna(0).astype(int) + 1
travel_nodes['node_size_visual'] = travel_nodes['centrality_quantile'] * 15

fig, ax = plt.subplots(figsize=(16, 16))

# Use city2graph to plot the network structure
# Node color represents centrality (skews handled by quantiles)
# Edge width represents frequency
c2g.plot_graph(
    nodes=travel_nodes,
    edges=travel_edges,
    markersize=8,
    node_color='centrality_quantile',
    node_alpha=0.9,
    edge_color='#bbbbbb',
    edge_linewidth=travel_edges['frequency'] / 200,
    edge_alpha=0.7,
    figsize=(16, 16),
    title="Central London Transit Network\nBetweenness Centrality of Stops",
    legend=False,  # remove colorbar
    ax=ax
)

# Add legend for betweenness centrality ranges
cmap = plt.cm.viridis
norm = plt.Normalize(
    vmin=travel_nodes['centrality_quantile'].min(),
    vmax=travel_nodes['centrality_quantile'].max()
)

quantile_ranges = (
    travel_nodes
    .groupby('centrality_quantile')['betweenness_centrality']
    .agg(['min', 'max'])
    .reset_index()
)

handles = [
    mlines.Line2D(
        [],
        [],
        color=cmap(norm(row['centrality_quantile'])),
        marker='o',
        linestyle='',
        markersize=8,
        label=f"{row['min']:.6f} – {row['max']:.6f}"
    )
    for _, row in quantile_ranges.iterrows()
]
ax.legend(handles=handles, title='Betweenness Centrality', loc='lower right')

# Add basemap with a clean, modern style
ctx.add_basemap(ax, crs=travel_nodes.crs, source=ctx.providers.CartoDB.DarkMatter)

plt.show()

# Calculate quantiles for node size to handle skew in centrality values
# We bin the data into 5 quantiles and scale the size for visibility
filtered_travel_nodes['centrality_quantile'] = pd.qcut(
    filtered_travel_nodes['betweenness_centrality'], 
    q=5, 
    labels=False
) + 1
filtered_travel_nodes['node_size_visual'] = filtered_travel_nodes['centrality_quantile'] * 15

fig, ax = plt.subplots(figsize=(16, 16))

# Use city2graph to plot the network structure
# Node size represents centrality (skews handled by quantiles)
# Edge width represents frequency
c2g.plot_graph(
    nodes=filtered_travel_nodes,
    edges=filtered_travel_edges,
    node_color='centrality_quantile',
    markersize=100,
    node_alpha=0.9,
    edge_color='#bbbbbb',
    edge_linewidth=filtered_travel_edges['frequency'] / 200,
    edge_alpha=0.7,
    ax=ax
)

# Add basemap with a clean, modern style
ctx.add_basemap(ax, crs=filtered_travel_nodes.crs, source=ctx.providers.CartoDB.DarkMatter)

# Add legend for betweenness centrality ranges
cmap = plt.cm.viridis
norm = plt.Normalize(
    vmin=filtered_travel_nodes['centrality_quantile'].min(),
    vmax=filtered_travel_nodes['centrality_quantile'].max()
)

quantile_ranges = (
    filtered_travel_nodes
    .groupby('centrality_quantile')['betweenness_centrality']
    .agg(['min', 'max'])
    .reset_index()
)

handles = [
    mlines.Line2D(
        [],
        [],
        color=cmap(norm(row['centrality_quantile'])),
        marker='o',
        linestyle='',
        markersize=8,
        label=f"{row['min']:.6f} – {row['max']:.6f}"
    )
    for _, row in quantile_ranges.iterrows()
]
ax.legend(handles=handles, title='Betweenness Centrality', loc='lower right')

plt.show()

# Keep walking-network download limited for performance
analysis_radius_m = 10_000  # 10 km radius around Piccadilly Circus

# Piccadilly Circus (WGS84 input coordinate)
piccadilly_circus_lon, piccadilly_circus_lat = -0.1337, 51.5101
piccadilly_circus_point = gpd.GeoSeries.from_xy(
    [piccadilly_circus_lon],
    [piccadilly_circus_lat],
    crs="EPSG:4326",
).to_crs(epsg=27700).iloc[0]

analysis_area = gpd.GeoDataFrame(
    {"name": ["Piccadilly Circus 10km"]},
    geometry=[piccadilly_circus_point.buffer(analysis_radius_m)],
    crs="EPSG:27700",
)
analysis_poly_wgs84 = analysis_area.to_crs(epsg=4326).geometry.iloc[0]

# Configure OSMnx for reliable downloads
ox.settings.use_cache = True
ox.settings.cache_folder = "./cache/osmnx"
ox.settings.requests_timeout = 180

print("Fetching Central London walk network from OSM...")
print(f"Query polygon area: {analysis_area.geometry.iloc[0].area / 1e6:.1f} km^2")

street_graph = ox.graph_from_polygon(
    analysis_poly_wgs84,
    network_type="walk",
    simplify=True,
    retain_all=True,
    truncate_by_edge=True,
)

if street_graph.number_of_nodes() == 0 or street_graph.number_of_edges() == 0:
    raise RuntimeError(
        "OSM download returned an empty graph. "
        "Try a larger area or check connectivity/API availability."
    )

print(
    f"OSM network downloaded successfully: "
    f"{street_graph.number_of_nodes():,} nodes, {street_graph.number_of_edges():,} edges"
)

# Convert to GeoDataFrames and reproject to British National Grid
street_nodes, street_edges = c2g.nx_to_gdf(street_graph)
street_nodes = street_nodes.to_crs(epsg=27700)
street_edges = street_edges.to_crs(epsg=27700)

# Snap isochrone center to the nearest street node around Piccadilly Circus
nearest_node_idx = street_nodes.geometry.distance(piccadilly_circus_point).idxmin()
center_point = street_nodes.loc[nearest_node_idx, "geometry"]
snap_distance_m = piccadilly_circus_point.distance(center_point)

# Compute walking travel time (seconds) assuming 4.5 km/h = 4500 m/h
street_edges["length"] = street_edges.geometry.length
street_edges["travel_time_sec"] = street_edges["length"] / 4500 * 3600

print(f"Isochrone origin snapped to nearest street node ({snap_distance_m:.1f} m from Piccadilly Circus)")
print(f"Street nodes (limited area): {len(street_nodes):,}")
print(f"Street edges (limited area): {len(street_edges):,}")
display(street_edges.head())

walking_labels = ["5 minutes", "10 minutes", "15 minutes"]
walking_threshold_seconds = [300, 600, 900]
# Compute layered walking isochrones in one call
walking_isochrones = c2g.create_isochrone(
    nodes=street_nodes,
    edges=street_edges,
    center_point=center_point,
    threshold=walking_threshold_seconds,
    edge_attr="travel_time_sec",
    method="concave_hull_knn",
    k=50
 )

walking_isochrones["label"] = walking_labels

for _, row in walking_isochrones.iterrows():
    area_km2 = row.geometry.area / 1e6
    print(f"{row['label']:12} -> {area_km2:6.2f} km^2 reachable area")

fig, ax = plt.subplots(figsize=(16, 14))

color_by_label = {
    "5 minutes": "#ff6b6b",
    "10 minutes": "#48dbfb",
    "15 minutes": "#1dd1a1",
}

# Plot larger thresholds first for layering
for _, row in walking_isochrones.sort_values("threshold", ascending=False).iterrows():
    row_gdf = gpd.GeoDataFrame([row], geometry="geometry", crs=walking_isochrones.crs)
    row_gdf.plot(
        ax=ax,
        alpha=0.35,
        color=color_by_label.get(row["label"], "#1dd1a1"),
        edgecolor="#111111",
        linewidth=1.25,
        label=row["label"],
        zorder=2,
    )

# Origin marker
ax.plot(center_point.x, center_point.y, "r*", markersize=28, markeredgecolor="white", label="Origin", zorder=4)

# Basemap for context
ctx.add_basemap(ax, crs=street_nodes.crs, source=ctx.providers.CartoDB.Positron, alpha=0.6)

ax.set_title("Walking Accessibility Isochrones\nWalk Travel Times from Center of London", fontsize=16, color="black", pad=18)
ax.legend(loc="upper right", framealpha=0.9)
ax.axis("off")
plt.tight_layout()
plt.show()

multimodal_labels = ["5 minutes", "10 minutes", "15 minutes"]
multimodal_threshold_seconds = [300, 600, 900]

# Compute layered multimodal isochrones in one call
multimodal_isochrones = c2g.create_isochrone(
    nodes=hetero_nodes,
    edges=hetero_edges,
    center_point=center_point,
    threshold=multimodal_threshold_seconds,
    edge_attr="travel_time_sec",
    cut_edge_types=[("bus_station", "is_next_to", "bus_station")],  # Only cut walking edges for multimodal isochrone
    method="concave_hull_knn",
    k=50
 )

multimodal_isochrones["label"] = multimodal_labels

for _, row in multimodal_isochrones.iterrows():
    area_km2 = row.geometry.area / 1e6
    print(f"{row['label']:12} -> {area_km2:6.2f} km^2 reachable area")

fig, ax = plt.subplots(figsize=(16, 14))

color_by_label = {
    "5 minutes": "#ff6b6b",
    "10 minutes": "#48dbfb",
    "15 minutes": "#1dd1a1",
}

# Plot larger thresholds first for layering
for _, row in multimodal_isochrones.sort_values("threshold", ascending=False).iterrows():
    row_gdf = gpd.GeoDataFrame([row], geometry="geometry", crs=multimodal_isochrones.crs)
    row_gdf.plot(
        ax=ax,
        alpha=0.35,
        color=color_by_label.get(row["label"], "#1dd1a1"),
        edgecolor="#111111",
        linewidth=1.25,
        label=row["label"],
        zorder=2,
    )

# Origin marker
ax.plot(center_point.x, center_point.y, "r*", markersize=28, markeredgecolor="white", label="Origin", zorder=4)

# Basemap for context
ctx.add_basemap(ax, crs=street_nodes.crs, source=ctx.providers.CartoDB.Positron, alpha=0.6)

ax.set_title("Multimodal Accessibility Isochrones\nWalk + Transit Travel Times from Center of London", fontsize=16, color="black", pad=18)
ax.legend(loc="upper right", framealpha=0.9)
ax.axis("off")
plt.tight_layout()
plt.show()

# Use meaningful London locations as origins (WGS84 lon/lat)
london_center_specs = [
    ("King's Cross St Pancras", -0.1238, 51.5308),
    ("Waterloo Station", -0.1132, 51.5033),
    ("Canary Wharf", -0.0195, 51.5054),
]

center_labels = [name for name, _, _ in london_center_specs]
center_lons = [lon for _, lon, _ in london_center_specs]
center_lats = [lat for _, _, lat in london_center_specs]

candidate_centers = gpd.GeoSeries.from_xy(
    center_lons,
    center_lats,
    crs="EPSG:4326",
).to_crs(street_nodes.crs)

# Snap each candidate center to nearest street node
multi_centers = []
for pt in candidate_centers:
    nearest_idx = street_nodes.geometry.distance(pt).idxmin()
    multi_centers.append(street_nodes.loc[nearest_idx, "geometry"])

# Compute one isochrone per center, then combine
iso_parts = []
for label, pt in zip(center_labels, multi_centers, strict=False):
    iso = c2g.create_isochrone(
        nodes=street_nodes,
        edges=street_edges,
        center_point=pt,
        threshold=900,  # 15 minutes
        edge_attr="travel_time_sec",
        method="concave_hull_knn",
        k=50,
    ).reset_index(drop=True)
    iso["center_label"] = label
    iso_parts.append(iso)

multi_center_isochrones = gpd.GeoDataFrame(
    pd.concat(iso_parts, ignore_index=True),
    geometry="geometry",
    crs=iso_parts[0].crs,
)

for _, row in multi_center_isochrones.iterrows():
    area_km2 = row.geometry.area / 1e6
    print(f"{row['center_label']}: {area_km2:.2f} km^2 reachable in 15 minutes")

fig, ax = plt.subplots(figsize=(16, 14))

multi_colors = ["#f97316", "#0ea5e9", "#22c55e"]

for i, row in multi_center_isochrones.iterrows():
    row_gdf = gpd.GeoDataFrame([row], geometry="geometry", crs=multi_center_isochrones.crs)
    row_gdf.plot(
        ax=ax,
        alpha=0.30,
        color=multi_colors[i % len(multi_colors)],
        edgecolor="#111111",
        linewidth=1.2,
        label=row["center_label"],
        zorder=2,
    )

for i, (pt, label) in enumerate(zip(multi_centers, center_labels, strict=False)):
    ax.plot(
        pt.x,
        pt.y,
        marker="*",
        markersize=18,
        color=multi_colors[i % len(multi_colors)],
        markeredgecolor="white",
        markeredgewidth=1.0,
        zorder=4,
    )
    ax.text(
        pt.x + 120,
        pt.y + 120,
        label,
        fontsize=10,
        color="#111111",
        zorder=5,
    )

ctx.add_basemap(ax, crs=street_nodes.crs, source=ctx.providers.CartoDB.Positron, alpha=0.6)
ax.set_title("Multiple-Center Isochrones (15-minute walk)\nCentral London", fontsize=16, color="black", pad=18)
ax.legend(loc="upper right", framealpha=0.9)
ax.axis("off")
plt.tight_layout()
plt.show()