Grids and Segmentations - Geospatial Data Cookbook

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

# 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

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

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

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

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

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

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

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

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

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