Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Grids and Segmentations

Authors
Affiliations
TU Wien
TU Wien
TU Wien
Binder

As you already learned in this cookbook, raw spatial data, such as irregular street segments and scattered intersection points, cannot be fed directly into most standard Machine Learning algorithms. To make this data meaningful for ML models, we need to transform it into a structured format. In this exercise we will do this by dividing the area of interest into regular cells using segmentation and computing aggregate statistics for each of those cells. This creates a standardized matrix of features, made up of our grid cells, that ML models can easily process.

So our next processing step will be to create two different grids for the area inside our bounding box:

  • A hexagonal grid with a side length of 250m

  • A square grid with a side length of 500m

What is Segmentation and why do we need it?

Segmentation refers to dividing a large geographic area into smaller, regularly shaped cells, typically squares or hexagons.

We need segmentation because raw spatial data doesn’t work well with most Machine Learning models. In order to be useful for the ML model, we need to pre-process the data in a way it can easily understand. By placing a grid over the area, we can group (aggregate) the data inside each cell. For example, we can count the number of intersections within each cell. This turns messy spatial data into a clean table where each cell is one row and the calculated values are features.

Square grids are simple and easy to work with, while hexagonal grids often give a more balanced view of neighboring areas. Therefore, the choice of grid type can influence the results which should be taken into consideration when choosing a certain grid.

Another important factor is the grid size. Using different grid sizes (like our 250m hexagons and 500m squares) helps capture patterns at different levels, and while finer grids will catch more details, coarser grids will be better for displaying general trends in the data.

Imports and Configuration

import sys
from pathlib import Path

import branca.colormap as cm
import folium
import geopandas as gpd
import pandas as pd
import requests
from shapely import wkb
from shapely.geometry import box

# ruff: noqa: D205

# Importing custom functions
sys.path.append(str(Path("..").resolve()))
from src.geoai.cookbook_functions import (
    create_hexagonal_grid,
    create_polygon_grid,
    get_base_map,
)

# Configuration
CITY_NAME = "Vienna"
BBOX = "16.335005,48.187854,16.400923,48.209995"
DATA_DIR = "./osm_data"
OUTPUT_DIR = "./output"

# Computing the center of the bounding box
min_lon, min_lat, max_lon, max_lat = map(float, BBOX.split(","))
bounds = [[min_lat, min_lon], [max_lat, max_lon]]
center_lat = (min_lat + max_lat) / 2
center_lon = (min_lon + max_lon) / 2

Now we load both feature tables from GitLab LFS (with a local cache): OSM intersections (intersections.csv) and GTFS stops (gtfs_stops.csv).

# Load OSM intersection features
csv_url = "https://gitlab.tuwien.ac.at/api/v4/projects/15378/repository/files/intersections.csv/raw?ref=main&lfs=true"
PROJECT_ROOT = Path.cwd()

csv_path = PROJECT_ROOT / "output" / "Vienna" / "intersections.csv"
csv_path.parent.mkdir(parents=True, exist_ok=True)

if not csv_path.exists():
    response = requests.get(csv_url, timeout=30)
    response.raise_for_status()
    csv_path.write_bytes(response.content)

df_int = pd.read_csv(csv_path)
df_int["geometry"] = df_int["geom"].apply(lambda x: wkb.loads(bytes.fromhex(x)))
intersections_gdf = gpd.GeoDataFrame(df_int, geometry="geometry", crs="EPSG:4326")

# Load GTFS stop features from GitLab LFS (same pattern as intersections.csv)
gtfs_csv_url = (
    "https://gitlab.tuwien.ac.at/api/v4/projects/15378/repository/files/"
    "gtfs_stops.csv/raw?ref=main&lfs=true"
)
gtfs_csv = Path(OUTPUT_DIR) / CITY_NAME / "gtfs_stops.csv"
gtfs_csv.parent.mkdir(parents=True, exist_ok=True)
if not gtfs_csv.exists():
    response = requests.get(gtfs_csv_url, timeout=60)
    response.raise_for_status()
    gtfs_csv.write_bytes(response.content)

gtfs_df = pd.read_csv(gtfs_csv)
stops_gdf = gpd.GeoDataFrame(
    gtfs_df,
    geometry=gpd.points_from_xy(gtfs_df["stop_lon"], gtfs_df["stop_lat"]),
    crs="EPSG:4326",
)
print(
    f"Loaded {len(intersections_gdf)} OSM intersections and {len(stops_gdf)} GTFS stops"
)
Loaded 8155 OSM intersections and 395 GTFS stops

Grids

In the next cell, we set up our two different grid. As already mentioned, those are:

  • A hexagonal grid with a side length of 250m

  • A square grid with a side length of 500m

# Determine Map Bounds from the intersections (with 100m buffer)
intersection_meters = intersections_gdf.to_crs("EPSG:3857")
minx, miny, maxx, maxy = intersection_meters.total_bounds
total_area = box(minx, miny, maxx, maxy)

# Hex Grid (250m)
hex_grid = create_hexagonal_grid(total_area, side_length=250)
hex_grid.set_crs("EPSG:3857", inplace=True)
hex_grid_geo = hex_grid.to_crs("EPSG:4326")

# Rect Grid (500m)
minx, miny, maxx, maxy = total_area.bounds
rect_grid = create_polygon_grid(
    width=maxx - minx,
    height=maxy - miny,
    cell_size=(500, 500),
    origin=(minx, miny),
    crs="EPSG:3857",
)
rect_grid_geo = rect_grid.to_crs("EPSG:4326")

Now we can plot the two grids we just created on the map to have a look at them and see how they divide the area.

# Creating the hexagonal grid
m_hex_grid = get_base_map(center_lat, center_lon)

folium.GeoJson(hex_grid_geo).add_to(m_hex_grid)

m_hex_grid
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...
# Creating the rectangular grid
m_rect_grid = get_base_map(center_lat, center_lon)

folium.GeoJson(rect_grid_geo).add_to(m_rect_grid)

m_rect_grid
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...

Aggregation

Now that we have created our grids, we bring them together with both point layers from feature engineering.

We need to figure out which points fall within which grid cell. This process is called a “Spatial Join”. Once we know which points belong to which cell, we can summarize the data for that specific area.

For every cell we will calculate:

From OSM intersections

  • Density: total number of intersections inside the cell

  • Average irregularity: average “skewness” (delta) of those intersections

From GTFS stops

  • Stop density: number of stops inside the cell

  • Average route connectivity: mean route_count of stops in the cell

Since we want the same calculation for hexagonal and rectangular grids, we define small helpers and apply them to both grids.

def aggregate_intersection_stats(
    grid_gdf: gpd.GeoDataFrame, intersections_gdf: gpd.GeoDataFrame
) -> gpd.GeoDataFrame:
    """Link intersections to grid cells using a Spatial Join."""
    joined = gpd.sjoin(grid_gdf, intersections_gdf, how="left", predicate="contains")
    stats = (
        joined.groupby(joined.index)
        .agg({"osm_id": "count", "delta": "mean"})
        .rename(columns={"osm_id": "total_count", "delta": "avg_delta"})
    )
    if "i_tpe" in joined.columns:
        t_counts = joined.pivot_table(
            index=joined.index, columns="i_tpe", aggfunc="size", fill_value=0
        )
        t_counts.columns = [f"count_{c}" for c in t_counts.columns]
        stats = stats.join(t_counts, how="left")
    grid_out = grid_gdf.merge(stats, left_index=True, right_index=True)
    cnt_cols = ["total_count"] + [c for c in grid_out.columns if "count_" in c]
    grid_out[cnt_cols] = grid_out[cnt_cols].fillna(0)
    return grid_out


def aggregate_stop_stats(
    grid_gdf: gpd.GeoDataFrame, stops: gpd.GeoDataFrame
) -> gpd.GeoDataFrame:
    """Link GTFS stops to grid cells using a Spatial Join."""
    joined = gpd.sjoin(grid_gdf, stops, how="left", predicate="contains")
    stats = joined.groupby(joined.index).agg(
        stop_count=("stop_id", "count"), avg_route_count=("route_count", "mean")
    )
    out = grid_gdf.merge(stats, left_index=True, right_index=True, how="left")
    out["stop_count"] = out["stop_count"].fillna(0).astype(int)
    out["avg_route_count"] = out["avg_route_count"].fillna(0.0)
    return out


hex_data = aggregate_intersection_stats(hex_grid_geo, intersections_gdf)
rect_data = aggregate_intersection_stats(rect_grid_geo, intersections_gdf)
hex_gtfs = aggregate_stop_stats(hex_grid_geo, stops_gdf)
rect_gtfs = aggregate_stop_stats(rect_grid_geo, stops_gdf)

Visualisation

We plot the aggregated values onto the hexagonal and rectangular grids. For each view we show the OSM metric and the matching GTFS metric on the same grid type, so you can compare street-network structure and transit coverage side by side.

def add_grid(
    data: pd.DataFrame,
    col: str,
    name: str,
    cmap_colors: list[str],
    m: folium.Map,
    tooltip_fields: list[str] | None = None,
    tooltip_aliases: list[str] | None = None,
) -> None:
    """Add a Choropleth layer to the map."""
    actual_colors = cmap_colors
    show = True

    # Create Legend
    vmin = data[col].min()
    vmax = data[col].max()

    # Handle completely empty data
    if pd.isna(vmin) or pd.isna(vmax):
        vmin, vmax = 0, 1

    colormap = cm.LinearColormap(colors=actual_colors, vmin=vmin, vmax=vmax)
    colormap.caption = name

    m.add_child(colormap)

    fields = tooltip_fields or ["total_count", "avg_delta"]
    aliases = tooltip_aliases or ["Count:", "Avg Delta:"]
    fields = [f for f in fields if f in data.columns]
    aliases = aliases[: len(fields)]

    # Create Layer
    geojson_layer = folium.GeoJson(
        data,
        style_function=lambda feature: {
            "fillColor": colormap(feature["properties"][col])
            if feature["properties"][col] is not None
            else "transparent",
            "color": "gray",
            "weight": 0.5,
            "fillOpacity": 0.7 if feature["properties"][col] is not None else 0,
        },
        tooltip=folium.GeoJsonTooltip(fields=fields, aliases=aliases)
        if fields
        else None,
    )

    fg = folium.FeatureGroup(name=name, show=show)
    geojson_layer.add_to(fg)
    fg.add_to(m)

We also keep the raw OSM intersection points available as an overlay on the choropleths.

# Create a feature group for all intersection points
fg_pts = folium.FeatureGroup(name="Points", show=True)

for _, row in intersections_gdf.iterrows():
    # Color logic
    color = "blue"
    if row.get("i_tpe") == "Car":
        color = "crimson"
    elif row.get("i_tpe") == "Path":
        color = "green"

    # Popup
    popup_html = f"""
    <div style="font-family: sans-serif; font-size: 12px;">
        <b>ID:</b> {row.get("osm_id")}<br>
        <b>Type:</b> {row.get("i_tpe")}<br>
        <b>Arms:</b> {row.get("num_ways")}<br>
        <b>Delta:</b> {row.get("delta"):.2f}<br>
    </div>
    """

    folium.CircleMarker(
        [row.geometry.y, row.geometry.x],
        radius=1,
        color=color,
        fill=True,
        popup=folium.Popup(popup_html, max_width=200),
    ).add_to(fg_pts)

Density on the hexagonal grid

Intersection density (OSM) and stop density (GTFS) on the 250 m hexagons.

m_hex_osm = get_base_map(center_lat, center_lon)
add_grid(
    data=hex_data,
    col="total_count",
    name="OSM intersection density",
    cmap_colors=["#f2f0f7", "#54278f"],
    m=m_hex_osm,
)
fg_pts.add_to(m_hex_osm)
m_hex_osm
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...
m_hex_gtfs = get_base_map(center_lat, center_lon)
add_grid(
    data=hex_gtfs,
    col="stop_count",
    name="GTFS stop density",
    cmap_colors=["#fff5eb", "#d94801"],
    m=m_hex_gtfs,
    tooltip_fields=["stop_count", "avg_route_count"],
    tooltip_aliases=["Stops:", "Avg routes:"],
)
m_hex_gtfs
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...

Density on the rectangular grid

The same comparison on 500 m squares.

m_rect_osm = get_base_map(center_lat, center_lon)
add_grid(
    data=rect_data,
    col="total_count",
    name="OSM intersection density",
    cmap_colors=["#f2f0f7", "#54278f"],
    m=m_rect_osm,
)
fg_pts.add_to(m_rect_osm)
m_rect_osm
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...
m_rect_gtfs = get_base_map(center_lat, center_lon)
add_grid(
    data=rect_gtfs,
    col="stop_count",
    name="GTFS stop density",
    cmap_colors=["#fff5eb", "#d94801"],
    m=m_rect_gtfs,
    tooltip_fields=["stop_count", "avg_route_count"],
    tooltip_aliases=["Stops:", "Avg routes:"],
)
m_rect_gtfs
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...

Connectivity / irregularity on the hexagonal grid

OSM average intersection irregularity (delta) next to GTFS average routes per stop.

m_hex_delta = get_base_map(center_lat, center_lon)
add_grid(
    data=hex_data,
    col="avg_delta",
    name="OSM avg delta",
    cmap_colors=["#f2f0f7", "#54278f"],
    m=m_hex_delta,
)
fg_pts.add_to(m_hex_delta)
m_hex_delta
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...
m_hex_routes = get_base_map(center_lat, center_lon)
add_grid(
    data=hex_gtfs,
    col="avg_route_count",
    name="GTFS avg routes/stop",
    cmap_colors=["#f7fbff", "#08306b"],
    m=m_hex_routes,
    tooltip_fields=["stop_count", "avg_route_count"],
    tooltip_aliases=["Stops:", "Avg routes:"],
)
m_hex_routes
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...

Connectivity / irregularity on the rectangular grid

m_rect_delta = get_base_map(center_lat, center_lon)
add_grid(
    data=rect_data,
    col="avg_delta",
    name="OSM avg delta",
    cmap_colors=["#f2f0f7", "#54278f"],
    m=m_rect_delta,
)
fg_pts.add_to(m_rect_delta)
m_rect_delta
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...
m_rect_routes = get_base_map(center_lat, center_lon)
add_grid(
    data=rect_gtfs,
    col="avg_route_count",
    name="GTFS avg routes/stop",
    cmap_colors=["#f7fbff", "#08306b"],
    m=m_rect_routes,
    tooltip_fields=["stop_count", "avg_route_count"],
    tooltip_aliases=["Stops:", "Avg routes:"],
)
m_rect_routes
/opt/conda/lib/python3.12/site-packages/folium/raster_layers.py:130: UserWarning: CartoDB tiles now require an API key. Please provide one to continue using the tiles. You can request the key at https://carto.com/basemaps/apikey/.
  tiles = tiles.build_url(fill_subdomain=False, scale_factor="{r}")  # type: ignore
Loading...

Summary

In this notebook you learned:

  • What Segmentation is and why it is necessary

  • Types of grids and how to create them

  • How to spatial-join and aggregate OSM and GTFS point features onto the same grids

  • How to compare those aggregates side by side on Folium choropleths

Introduction
Feature Engineering
Introduction
Visualisation