How to Create a Morphological Tessellation in Python

Problem statement

Many urban-form measures need a unit of land around each building: the plot. Cadastral plots are often unavailable, inconsistent between countries, or subdivided in ways that have nothing to do with form. Morphological tessellation replaces them with a geometric proxy: divide all the space in an area among its buildings, giving each building the land that is closer to it than to any other.

momepy builds it in one call. What takes care is the extent — without a limit, cells at the edge of the study area run off to infinity — and the choice between a free tessellation and one bounded by streets, which changes what a cell means.

Quick answer

Build a limit, then tessellate; or tessellate within street-bounded enclosures:

import numpy as np, geopandas as gpd, osmnx as ox, momepy
from shapely.geometry import box

b = gpd.read_parquet("denhaag_bag.parquet")                       # footprints guide, EPSG:28992
cx, cy = 81000, 455000                                             # central Den Haag, 1 km square
aoi = box(cx - 500, cy - 500, cx + 500, cy + 500)
b = b[b.representative_point().within(aoi) & b.is_valid].reset_index(drop=True)
b = b[b.area > 10].reset_index(drop=True)

limit = momepy.buffered_limit(b, buffer=100)                      # land within 100 m of a building
cells = momepy.morphological_tessellation(b, clip=limit)
print(f"{len(b):,} buildings → {len(cells):,} cells")
print(f"cells cover {cells.area.sum() / 1e6:.3f} km², buildings {b.area.sum() / 1e6:.3f} km²")
print(f"median cell {cells.area.median():.0f} m², median building {b.area.median():.0f} m²")

1,990 buildings → 1,990 cells
cells cover 1.450 km², buildings 0.479 km²
median cell 223 m², median building 114 m²

One cell per building, partitioning all the land within the limit. The ratio of the two areas — building footprint over cell — is the coverage ratio of each building's "plot".

Clean footprints, choose a limit, tessellate, verify

1. Clean the footprints

Tessellation needs valid, non-overlapping polygons. Overlaps and duplicates produce empty or split cells; tiny slivers produce tiny cells. Fix or drop them first.

2. Choose the limit

A buffered limit (the union of buffers around buildings) keeps cells from extending into fields and water. Its distance is a parameter: 100 m is common in dense cities, larger in sparse suburbs.

3. Tessellate

momepy.morphological_tessellation erodes each footprint slightly, samples points along the boundaries and builds Voronoi cells from them, then dissolves them per building.

4. Or bound it by streets

momepy.enclosures splits the area by streets (and optionally rivers and railways); momepy.enclosed_tessellation then tessellates within each enclosure, so no cell crosses a street. The result is much closer to a plot.

5. Verify

Check that every building has exactly one cell, that cells do not overlap, and that their union matches the limit. momepy.verify_tessellation reports the failures.

Enclosed tessellation and coverage per cell

The data used on this page

Building footprints with LiDAR heights from the 3D BAG (the Dutch building register with heights from AHN), fetched from its WFS for central Den Haag — 4.285–4.330° E, 52.065–52.095° N — exactly as in How to download building footprints from OSM and Overture, and saved as denhaag_bag.parquet in the Dutch national grid (EPSG:28992). The OpenStreetMap answers are read as the Overpass mirror gave them on 10 October 2026.

import io
import os
import zipfile

import geopandas as gpd
import requests
from pyproj import Transformer

# OpenStreetMap as Overpass returned it on 10 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-morphological-tessellation-python-2026-10-10.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

if not os.path.exists("denhaag_bag.parquet"):
    W, S, E, N = 4.285, 52.065, 4.330, 52.095
    to_rd = Transformer.from_crs(4326, 28992, always_xy=True)
    (x0, y0), (x1, y1) = to_rd.transform(W, S), to_rd.transform(E, N)
    feats, start = [], 0
    while True:
        fc = requests.get("https://data.3dbag.nl/api/BAG3D/wfs", params={
            "service": "WFS", "version": "2.0.0", "request": "GetFeature", "typeNames": "BAG3D:lod12",
            "count": 5000, "startIndex": start, "outputFormat": "application/json",
            "bbox": f"{x0},{y0},{x1},{y1},EPSG:28992"}, timeout=300).json()
        feats += fc["features"]
        if len(fc["features"]) < 5000:
            break
        start += 5000
    bag = gpd.GeoDataFrame.from_features(feats, crs=28992)
    bag["height"] = bag.b3_h_70p - bag.b3_h_maaiveld
    bag[["identificatie", "b3_h_70p", "b3_h_maaiveld", "height", "b3_opp_grond", "geometry"]] \
        .to_parquet("denhaag_bag.parquet")
print(len(gpd.read_parquet("denhaag_bag.parquet")), "buildings")

23929 buildings

Example 1 — enclosed tessellation

ox.settings.http_user_agent = "spatialworkflow-examples/1.0"
ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"   # mirror; overpass-api.de refused
ox.settings.overpass_rate_limit = False
streets = ox.graph_to_gdfs(ox.graph_from_point((52.0799, 4.3113), dist=1200,
                                                network_type="drive"), nodes=False).to_crs(28992)
encl = momepy.enclosures(streets, limit=gpd.GeoSeries([aoi], crs=28992))
etess = momepy.enclosed_tessellation(b, encl, n_jobs=1)
built = etess[etess.index >= 0]
print(f"{len(encl)} enclosures; {len(built):,} building cells, "
      f"{(etess.index < 0).sum()} empty enclosures kept as cells")
print(f"median cell {built.area.median():.0f} m² (free tessellation: {cells.area.median():.0f} m²)")

491 enclosures; 2,058 building cells, 418 empty enclosures kept as cells
median cell 205 m² (free tessellation: 223 m²)

Bounding by streets shrinks the cells along every street frontage: the street space no longer belongs to the nearest building. There are also 68 more building cells than buildings, because a footprint that crosses an enclosure line — a building over a passage, a street centre-line drawn through a courtyard — gets a cell on each side. Enclosures without buildings — squares, parks — survive as cells with no building attached.

Example 2 — coverage per cell

cells["building_area"] = b.area.reindex(cells.index)
cells["coverage"] = cells.building_area / cells.area
print(cells.coverage.describe(percentiles=[.1, .5, .9]).round(2).to_string())

count    1990.00
mean        0.53
std         0.22
min         0.00
10%         0.20
50%         0.55
90%         0.81
max         1.01

Half of the cells are more than 55% built on, and a tenth more than 81%. The maximum a hair above 1 is geometry, not data: the cell is built from a slightly simplified footprint. Coverage per cell is the most basic density measure and the one a tessellation makes possible without cadastral plots. Its distribution — not just its mean — distinguishes a perimeter-block district from one of towers in open space.

Example 3 — verify the tessellation

problems = momepy.verify_tessellation(cells, b)
print(f"buildings without a cell: {len(set(b.index) - set(cells.index))}")
print(f"overlapping cell area: {cells.union_all().area - cells.area.sum():.1f} m² "
      f"(negative would mean overlap)")
print(f"cells ≠ limit: {abs(cells.union_all().area - limit.area):.1f} m²")

buildings without a cell: 0
overlapping cell area: 0.0 m² (negative would mean overlap)
cells ≠ limit: 0.0 m²

A tessellation that passes these checks is a partition: every square metre inside the limit belongs to exactly one building.

Voronoi around footprints, and limits as choices

Why Voronoi cells around footprints, not centroids

Voronoi cells around centroids give a big building the same share of land as a small one next to it. Sampling points along each footprint's boundary gives each building the land nearest its walls, which is what "its plot" intuitively means. Eroding the footprints slightly first keeps touching buildings from sharing boundary points.

Why the limit is a modelling choice

Without a limit, cells at the edge extend to the edge of the study area — or to infinity. With a buffered limit, the distance decides how much surrounding open land a building "owns". Density metrics per cell change with that distance, so it belongs in the method description.

Why enclosures are closer to plots

Plots are bounded by streets; a free tessellation ignores them and hands half of every street to the buildings on either side. Enclosed tessellation restores the street as a boundary. It needs a clean street network — gaps in the network merge enclosures across streets.

Performance, multipart buildings and barriers

• Performance: tessellation of tens of thousands of buildings takes minutes; n_jobs parallelises enclosed tessellation but uses every core by default.
• shrink and segment control erosion and sampling density; defaults of 0.4 m and 0.5 m suit metric data.
• Multipart buildings become one cell per row; explode or dissolve consistently first.
• Buildings touching the AOI edge: select by representative point, and treat edge cells as incomplete.
• Non-building barriers (rivers, railways) go into enclosures as additional barriers.

Internal links

• Morphological tessellation explained — what a cell stands for
• Tessellation fails or produces odd cells — the fix guide
• How to calculate coverage and floor area ratio — density per cell
• How to extract urban blocks from a street network — the street-bounded unit
• How to build Voronoi service areas — the underlying geometry
• How to download building footprints from OSM and Overture — the input

FAQ

What is a morphological tessellation?

A partition of the land around buildings into cells, each containing the space closer to one building's footprint than to any other.

How do I create one in Python?

momepy.morphological_tessellation(buildings, clip=limit), with a limit such as momepy.buffered_limit(buildings, buffer=100).

What is an enclosed tessellation?

A tessellation computed within street-bounded enclosures, so no cell crosses a street.

What buffer should the limit use?

Around 100 m in dense urban areas; larger for suburbs. Report it.

Why do some buildings have no cell?

Invalid or overlapping footprints, or duplicates. Clean the input and check with momepy.verify_tessellation.

Is a tessellation cell the same as a cadastral plot?

No — it is a geometric proxy, closest to a plot when bounded by streets.