Tutorial#

This notebook demonstrates the core functionality of Ferrobus, a high-performance transit routing library for Python.

Setup#

First, let’s import the necessary libraries:

import datetime
import json
import time

import geopandas as gpd
from shapely import wkt

import ferrobus
import warnings

# Suppress all warnings
warnings.filterwarnings("ignore")

1. Creating a Transit Model#

The transit model is the core component that represents the transportation network. It combines GTFS transit data with OpenStreetMap road network data.

# Create a transit model with both GTFS and OSM data
model = ferrobus.create_transit_model(
    max_transfer_time=1800,  # 30 minutes as max transfer time
    osm_path="your_data.osm.pbf",
    gtfs_dirs=[
        "/feed_1",
        "/feed_2",
    ],
    date=datetime.date(2024, 4, 1),  # date for service filtering
)
# Display model info
print(model)
TransitModel with 8691 stops, 73775 routes and 2020348 trips

2. Creating Transit Points#

Transit points represent geographic coordinates in the model, which can be used as origins or destinations for routing. Those pre-computed points allow for extremly fast routing queries with RAPTOR without the need to compute same walking paths from point to transit stops over and over again.

# Create transit points with different parameters
origin = ferrobus.create_transit_point(
    lat=60.04,
    lon=30.34,
    transit_model=model,
    max_walking_time=1200,  # 20 minutes maximum walking time
    max_nearest_stops=10,  # Consider up to 10 nearest transit stops
)

# Create a simpler transit point with default parameters
destination = ferrobus.create_transit_point(lat=59.93, lon=30.37, transit_model=model)

3. Finding Routes#

The core functionality of Ferrobus is finding optimal routes between points.

# Find a route with basic parameters
start_time = time.perf_counter()
route = ferrobus.find_route(
    transit_model=model,
    start_point=origin,
    end_point=destination,
    departure_time=43200,  # 12:00 noon in seconds since midnight
    max_transfers=3,  # Allow up to 3 transfers
)
end_time = time.perf_counter()

# Display route information
print(f"Route found in {end_time - start_time:.3f} seconds")
print(f"Travel time: {route['travel_time_seconds'] / 60:.1f} minutes")
print(f"Transit time: {route['transit_time_seconds'] / 60:.1f} minutes")
print(f"Walking time: {route['walking_time_seconds'] / 60:.1f} minutes")
print(f"Number of transfers: {route['transfers']}")
Route found in 0.021 seconds
Travel time: 55.6 minutes
Transit time: 50.2 minutes
Walking time: 5.3 minutes
Number of transfers: 3

4. Detailed Journey Visualization#

For a more detailed view, we can get the full journey with all legs (walking, transit) as GeoJSON for visualization.

# Get detailed journey information as GeoJSON
journey = ferrobus.detailed_journey(
    transit_model=model,
    start_point=origin,
    end_point=destination,
    departure_time=43200,  # 12:00 noon
    max_transfers=3,
)

# Convert to GeoDataFrame for visualization
journey_gdf = gpd.GeoDataFrame.from_features(json.loads(journey), crs=4326)

# Display journey information
print(f"Journey has {len(journey_gdf)} legs")

# Visualize the journey
journey_gdf.explore(
    column="leg_type",
    cmap="Set1",
    legend=True,
    tiles="CartoDB Dark Matter",
    legend_kwds={"colorbar": False, "loc": "upper right"},
)
Journey has 9 legs
Make this Notebook Trusted to load map: File -> Trust Notebook

5. Travel Time Matrix#

Calculate travel times between multiple points.

# Show part of the matrix (first 3×3 submatrix)
from statistics import mean

import numpy as np
import pandas as pd

grid = gpd.read_file("/SPB_grid.gpkg")

points = [
    ferrobus.create_transit_point(point.y, point.x, model)
    for point in grid.geometry.centroid.to_list()
]

# Calculate travel time matrix
matrix = ferrobus.travel_time_matrix(
    transit_model=model,
    points=points,
    departure_time=43200,  # 12:00 noon
    max_transfers=3,
)

# Display matrix dimensions
print(f"Matrix size: {len(matrix)} × {len(matrix[0])}")

# Convert to DataFrame for better visualization
matrix_np = np.array(matrix)
matrix_df = pd.DataFrame(matrix_np[:3, :3])
print("Travel times (seconds):")
print(matrix_df)


def compute_mean_if_enough_data(values, threshold=600):
    # Filter out None values
    valid_values = [v for v in values if v is not None]
    # Only compute the mean if we have enough valid values
    if len(valid_values) > threshold:
        return mean(valid_values)
    return None


# Compute the mean for each row in the matrix and store in grid["mean_travel_time"]
grid["mean_travel_time"] = [compute_mean_if_enough_data(row) for row in matrix]

grid.explore(
    column="mean_travel_time",
    cmap="RdYlGn_r",
    legend=True,
    legend_title="Mean Travel Time (seconds)",
    scheme="quantiles",
    legend_kwds={"colorbar": False, "loc": "upper right"},
    k=8,
)
Matrix size: 963 × 963
Travel times (seconds):
      0     1     2
0     0  2369  5399
1  2631     0  5560
2  5916  5969     0
Make this Notebook Trusted to load map: File -> Trust Notebook

6. Isochrones#

Calculate isochrones (areas reachable within a given time) from a point.

# Create an isochrone index to speed up calculations
# The bounding box is defined as a WKT polygon string
bounding_box = grid.unary_union.convex_hull.wkt
index = ferrobus.create_isochrone_index(
    model, bounding_box, 9
)  # 9 represents the grid precision

# Calculate isochrones for different time thresholds
time_thresholds = [3600, 2400, 1800, 900]  # 30, 40, 50 minutes in seconds
isochrones = []

for threshold in time_thresholds:
    isochrone = ferrobus.calculate_isochrone(
        model,
        start=origin,
        departure_time=43200,  # 12:00 noon
        max_transfers=3,
        cutoff=threshold,
        index=index,
    )
    isochrones.append(isochrone)

# Convert to GeoDataFrame for visualization
isochrone_gdf = gpd.GeoDataFrame(
    {
        "time_threshold": [t / 60 for t in time_thresholds],  # Convert to minutes
        "geometry": [wkt.loads(poly) for poly in isochrones],
    },
    crs=4326,
)

# Visualize the isochrones
isochrone_gdf.explore(
    column="time_threshold",
    cmap="viridis_r",
    legend=True,
    legend_title="Reachable within (minutes)",
)
Snapped 9370 of 10459
Make this Notebook Trusted to load map: File -> Trust Notebook
pa_iso = ferrobus.calculate_percent_access_isochrone(
    model,
    start=origin,
    departure_range=(79200, 82800),  # From 22:00 to 23:00 in seconds
    sample_interval=60,
    max_transfers=3,
    cutoff=3600,
    index=index,
)

pa_iso = (
    gpd.GeoDataFrame.from_features(json.loads(pa_iso), crs=4326)
    .dissolve("percent_access")
    .reset_index()
)

pa_iso.explore(
    column="percent_access",
    cmap="RdYlGn",
    legend=True,
    legend_title="Percent chance of reaching this",
    scheme="natural_breaks",
    k=5,
    opacity=0.3,
    legend_kwds={"colorbar": False, "loc": "upper right"},
)
Make this Notebook Trusted to load map: File -> Trust Notebook