Spatial units

A spatial unit is the polygon everything is measured on and a class is assigned to. Every stage after this one joins on its identifier, unit_id.

Three strategies, all satisfying SpatialUnitStrategy. unit_strategy is config; the default is grid and there is deliberately no auto-selection.

Units must form a partition of the bounding box — they must tile it exactly, with no gaps and no overlaps. Assert it explicitly — a validity test is not a partition test, and enclosure faces formed outside the extent were measured at 222% of the bbox on Berlin and 379% on Rotterdam, silently corrupting the denominator of every area-weighted statistic downstream.

The three differ in scale, which is the thing to know before choosing one:

	median unit
EnclosureUnits, Hong Kong fixture	0.04 ha — a block
GridUnits, 100 m	1.00 ha
PatchUnits	11.69 ha
a WUDAPT polygon, across the study cities	2.2–52 ha
a So2Sat patch	10.24 ha

An enclosure is a city block, bounded by streets; an LCZ patch — the area a class is meant to describe — is a neighbourhood hundreds of metres across. A thinner barrier set does not close that gap — it stops subdividing big faces rather than enlarging small ones — so PatchUnits sets the scale with a merge step instead.

lczkit.units

Spatial-unit generation: EnclosureUnits, GridUnits, and aggregate() between them.

Also home to check_units(), the entry contract every stage keyed on unit_id enforces. It lives with the definition of the unit of exchange rather than with any one consumer.

check_units

check_units(units: GeoDataFrame) -> None

Raise unless units is projected, indexed by unit_id, and uniquely indexed.

Same checks aggregate() and height_metrics() make inline, for the same reason: a geographic CRS makes every area meaningless, and a missing or duplicated unit_id index makes the join silently wrong rather than loud.

Source code in src/lczkit/units/__init__.py

def check_units(units: gpd.GeoDataFrame) -> None:
    """Raise unless `units` is projected, indexed by `unit_id`, and uniquely indexed.

    Same checks `aggregate()` and `height_metrics()` make inline, for the same reason: a geographic
    CRS makes every area meaningless, and a missing or duplicated `unit_id` index makes the join
    silently wrong rather than loud.
    """
    assert_projected_crs(units, "units")
    if units.index.name != "unit_id":
        raise ValueError("units must be indexed by unit_id")
    if not units.index.is_unique:
        raise ValueError("units index must be unique")

Grid

lczkit.units.grid

GridUnits: a 100 m regular grid in the local UTM CRS.

The default, because it is the format every existing LCZ map, validation dataset and WRF workflow uses. lczkit.validation joins against the Demuzere global LCZ map on this same 100 m grid.

GridUnits

GridUnits(cell_size_m: float = DEFAULT_CELL_SIZE_M)

A regular cell_size_m grid, aligned to the local UTM CRS's own coordinate origin.

Aligned to that origin — not to bbox. Two overlapping bboxes therefore assign the same unit_id to the same real-world cell, which is what makes grid unit_ids meaningfully "stable" across runs and across cities sharing a UTM zone, unlike EnclosureUnits' barrier-dependent ids.

Cells are kept whole (never clipped to bbox) and are included if they intersect bbox at all, so the returned grid can extend slightly beyond bbox's edges. barriers is accepted for SpatialUnitStrategy conformance but ignored, per the protocol's own docstring.

Set the cell side in metres, defaulting to the 100 m every LCZ workflow uses.

100 m is what every existing LCZ map, validation dataset and WRF workflow is on, so a run at another size is comparable to nothing. Rejects a non-positive size rather than producing an empty or inverted grid downstream.

Source code in src/lczkit/units/grid.py

def __init__(self, cell_size_m: float = DEFAULT_CELL_SIZE_M) -> None:
    """Set the cell side in metres, defaulting to the 100 m every LCZ workflow uses.

    100 m is what every existing LCZ map, validation dataset and WRF workflow is on, so a
    run at another size is comparable to nothing. Rejects a non-positive size rather than
    producing an empty or inverted grid downstream.
    """
    if cell_size_m <= 0:
        raise ValueError(f"cell_size_m must be positive, got {cell_size_m}")
    self.cell_size_m = cell_size_m

generate

generate(bbox: BBox, barriers: GeoDataFrame | None = None) -> GeoDataFrame

Return the grid cells intersecting bbox, indexed by unit_id.

bbox is lon/lat and the returned frame is in the local UTM CRS, carrying geometry alone — every later stage joins onto it by unit_id. Cells are kept whole and are included if they intersect bbox at all, so the grid can extend slightly past its edges.

barriers is accepted for SpatialUnitStrategy conformance and ignored: a regular grid is defined by the CRS origin and the cell size, and nothing about the city moves it.

Source code in src/lczkit/units/grid.py

def generate(self, bbox: BBox, barriers: gpd.GeoDataFrame | None = None) -> gpd.GeoDataFrame:
    """Return the grid cells intersecting `bbox`, indexed by `unit_id`.

    `bbox` is lon/lat and the returned frame is in the local UTM CRS, carrying geometry
    alone — every later stage joins onto it by `unit_id`. Cells are kept whole and are
    included if they intersect `bbox` at all, so the grid can extend slightly past its edges.

    `barriers` is accepted for `SpatialUnitStrategy` conformance and ignored: a regular grid
    is defined by the CRS origin and the cell size, and nothing about the city moves it.
    """
    del barriers  # grid cells ignore barriers; see SpatialUnitStrategy protocol docstring
    crs = local_utm_crs(bbox)
    bbox_utm = gpd.GeoSeries([box(*bbox)], crs="EPSG:4326").to_crs(crs).iloc[0]
    minx, miny, maxx, maxy = bbox_utm.bounds
    cell = self.cell_size_m

    col_start, col_end = math.floor(minx / cell), math.floor(maxx / cell)
    row_start, row_end = math.floor(miny / cell), math.floor(maxy / cell)

    unit_ids: list[str] = []
    geoms: list[Polygon] = []
    for col in range(col_start, col_end + 1):
        x0 = col * cell
        for row in range(row_start, row_end + 1):
            y0 = row * cell
            cell_geom = box(x0, y0, x0 + cell, y0 + cell)
            if cell_geom.intersects(bbox_utm):
                unit_ids.append(f"grid_{col}_{row}")
                geoms.append(cell_geom)

    return gpd.GeoDataFrame({"unit_id": unit_ids}, geometry=geoms, crs=crs).set_index("unit_id")

Enclosures

lczkit.units.enclosures

EnclosureUnits: momepy.enclosures()-based spatial units, the GeoClimate RSU analogue.

The barrier set is streets, rail, waterbodies and large vegetation patches, and that list is exhaustive rather than illustrative.

Large vegetation patches have no vector source, so assemble_barriers accepts vegetation as an optional layer and omits it until a raster-derived vegetation-patch layer exists to pass in.

The one plausible-looking wrong answer: VectorSource.land_use() returns polygons that include parks and green space, and it is tempting to reach for them as the missing vegetation barrier. Do not. Land use carries functional semantics only — it feeds industrial_fraction — and it is neither a barrier here nor a land-cover source. Rasters own land cover, and the vegetation barrier must come from them.

EnclosureUnits

momepy.enclosures()-based spatial units.

All barrier linework is unioned into one series before polygonization (see momepy.elements.enclosures's source), so momepy's own primary_barriers / additional_barriers split — meaningful elsewhere for RSU semantics — makes no difference to the resulting enclosure geometries here; barriers is passed through as primary_barriers with no additional_barriers.

unit_ids are f"enclosure_{eid}", eid being the sequential integer momepy.enclosures assigns in polygonization order. This is deterministic for a fixed barriers input (same rows, same order, same library versions) but not stable against changes to the set of barriers — unlike GridUnits, whose ids are tied to absolute coordinates, an enclosure's id can shift if upstream barrier data changes. Documented, not fixed: no id scheme survives which regions of the barrier network happen to change.

Enclosures are restricted to bbox. Barrier linework runs past the study area — a rail line entering the extent continues well beyond it — and polygonizing it produces faces lying wholly outside. Left in, those are returned as units: on the Berlin fixture the result covered 222% of the requested extent and on Rotterdam 379%, so the units were not a partition and any area-weighted statistic over them was measured against the wrong denominator. clip=True keeps only the faces whose representative point falls inside limit, which is exact rather than approximate: limit's boundary is part of the noded linework, so every face is already split at the extent's edge and none straddles it. No geometry is modified.

generate

generate(bbox: BBox, barriers: GeoDataFrame | None = None) -> GeoDataFrame

Polygonize barriers into enclosures clipped to bbox, indexed by unit_id.

bbox is lon/lat; barriers must already be in the projected CRS the units come back in, and the returned frame carries geometry alone — every later stage joins onto it by unit_id.

Unlike GridUnits, barriers is required rather than optional: this strategy has no barrier source of its own, and silently returning one face for the whole bbox would be a partition with no information in it. Build the argument with assemble_barriers.

Source code in src/lczkit/units/enclosures.py

def generate(self, bbox: BBox, barriers: gpd.GeoDataFrame | None = None) -> gpd.GeoDataFrame:
    """Polygonize `barriers` into enclosures clipped to `bbox`, indexed by `unit_id`.

    `bbox` is lon/lat; `barriers` must already be in the projected CRS the units come back
    in, and the returned frame carries geometry alone — every later stage joins onto it by
    `unit_id`.

    Unlike `GridUnits`, `barriers` is required rather than optional: this strategy has no
    barrier source of its own, and silently returning one face for the whole bbox would be a
    partition with no information in it. Build the argument with `assemble_barriers`.
    """
    if barriers is None or barriers.empty:
        raise ValueError(
            "EnclosureUnits requires `barriers` (see `assemble_barriers`) — it has no "
            "barrier source of its own, unlike GridUnits, which barriers is optional for."
        )
    assert_projected_crs(barriers, "barriers")
    crs = barriers.crs
    assert crs is not None  # narrows for mypy; assert_projected_crs already guarantees this
    limit = gpd.GeoSeries([box(*bbox)], crs="EPSG:4326").to_crs(crs)

    raw: gpd.GeoDataFrame = momepy.enclosures(barriers, limit=limit, clip=True)
    raw["unit_id"] = "enclosure_" + raw["eID"].astype(str)
    return raw.set_index("unit_id")[["geometry"]]

assemble_barriers

assemble_barriers(streets: GeoDataFrame, waterbodies: GeoDataFrame, *, rail: GeoDataFrame | None = None, vegetation: GeoDataFrame | None = None) -> GeoDataFrame

Combine barrier layers into the single barriers frame EnclosureUnits.generate expects.

That is: one geometry column, no attributes, all layers in the same projected CRS.

rail and vegetation are optional — pass None (the default) when a layer isn't available; enclosures form from whatever barriers are given. All non-empty inputs must already share one CRS (typically the one clean_vectors() reprojected into); this function does not itself reproject anything.

These four layers are the whole eligible barrier set. In particular, do not pass CleanedVectors.land_use as vegetation — see this module's docstring.

Source code in src/lczkit/units/enclosures.py

def assemble_barriers(
    streets: gpd.GeoDataFrame,
    waterbodies: gpd.GeoDataFrame,
    *,
    rail: gpd.GeoDataFrame | None = None,
    vegetation: gpd.GeoDataFrame | None = None,
) -> gpd.GeoDataFrame:
    """Combine barrier layers into the single `barriers` frame `EnclosureUnits.generate` expects.

    That is: one geometry column, no attributes, all layers in the same projected CRS.

    `rail` and `vegetation` are optional — pass `None` (the default) when a layer isn't
    available; enclosures form from whatever barriers are given. All non-empty inputs must
    already share one CRS (typically the one `clean_vectors()` reprojected into); this function
    does not itself reproject anything.

    These four layers are the whole eligible barrier set. In particular, do not pass
    `CleanedVectors.land_use` as `vegetation` — see this module's docstring.
    """
    assert_projected_crs(streets, "streets")
    assert_projected_crs(waterbodies, "waterbodies")
    layers = [streets.geometry, waterbodies.geometry]
    crs = streets.crs
    for extra, name in ((rail, "rail"), (vegetation, "vegetation")):
        if extra is None or extra.empty:
            continue
        assert_projected_crs(extra, name)
        if extra.crs != crs:
            raise ValueError(f"{name}.crs ({extra.crs}) != streets.crs ({crs})")
        layers.append(extra.geometry)
    combined = pd.concat(layers, ignore_index=True)
    return gpd.GeoDataFrame(geometry=gpd.GeoSeries(combined, crs=crs))

Patches

lczkit.units.patches

PatchUnits: street-bounded blocks merged into units at the scale an LCZ patch is drawn at.

The measurement this exists to answer. Enclosures are a block, not a patch. On the Hong Kong fixture momepy.enclosures over the cleaned barrier set returns a median unit of 0.04 ha with 72.9% under 0.1 ha, against WUDAPT polygons that run 2.2-52 ha across the sixteen study cities and a So2Sat patch of 10.24 ha. Berlin shows the same from the other direction: 78% of its enclosures are sub-1000 m2 slivers.

Re-cutting the barrier set does not fix that, which was worth measuring before assuming it:

barriers (HK fixture)	seeds	median	p90	max	% < 0.1 ha
all streets	4095	0.04 ha	0.36 ha	673 ha	72.9%
drop footway/steps/path	769	0.11 ha	2.52 ha	691 ha	48.6%
major + tertiary	509	0.07 ha	3.85 ha	692 ha	56.2%
major only	397	0.07 ha	2.96 ha	692 ha	57.9%

Every barrier set is bimodal — slivers plus a handful of very large faces — because a thinner barrier network does not enlarge the small faces, it only stops subdividing the big ones. The scale is set by a merge step, not by which streets are barriers.

So: seeds, then merge. Stage 1 is EnclosureUnits unchanged, over a barrier set with the pedestrian classes removed — a footpath is not a boundary between two LCZ patches, and it is 50-73% of the network in Berlin, Hong Kong and Milan. Stage 2 merges each seed into the contiguous neighbour it most resembles until the units reach min_area_m2.

This does not replace GridUnits and is not the default. The strategy is config, the default is grid, and there is no auto-selection. Grid cells are what every published LCZ map, validation dataset and WRF workflow uses, and every published figure here is measured on them.

A caveat that must not be lost. The merge reads building surface fraction and mean height, and those are two of the seven dimensions the classifier then scores. Units defined partly by the quantity being measured is the standard shape of a regionalisation (SKATER, AZP and the rest work this way), and it is still a form of circularity: a patch is more homogeneous in BSF than a cell partly because it was built to be. It cannot inflate agreement with an external reference, which is what the validation measures, but it does mean bsf_by_reference_class on patch units is a weaker test than the same table on a grid, so read that table on the grid.

DEFAULT_MIN_AREA_M2 module-attribute

DEFAULT_MIN_AREA_M2 = 50000.0

5 ha - the median WUDAPT polygon across the sixteen study cities, which is the grain the reference is actually drawn at. A 100 m cell is 1 ha and a So2Sat patch 10.24 ha, so this sits between the two objects that are genuinely patch-scale.

A floor, not a centre, which is why it is named for what it does. Merging stops when a unit reaches the minimum, and the merge that gets it there overshoots, so the resulting median lands around twice this value: on the Hong Kong fixture 5 ha gives p10 5.75, median 10.5, p90 24.5 ha. Set it to roughly half the grain wanted.

DEFAULT_MAX_AREA_M2 module-attribute

DEFAULT_MAX_AREA_M2 = 500000.0

50 ha. A ceiling, not a target: it stops a merge chain swallowing a whole industrial estate or a park into one unit, which is the failure mode of merging on size alone.

PEDESTRIAN_CLASSES module-attribute

PEDESTRIAN_CLASSES: frozenset[str] = frozenset({'footway', 'steps', 'path', 'cycleway', 'bridleway'})

Overture road classes excluded from the barrier set by default.

Not a quality judgement about the data - these are real features, correctly mapped. They are excluded because an LCZ patch is a neighbourhood of homogeneous cover and a footpath through a housing estate does not divide one patch from another. They are also the majority of the network where they are mapped at all: 72.7% of Berlin's segments, 72.8% of Hong Kong's, 50.6% of Milan's, against 3.5-7.5% in cities where pedestrian mapping is sparser. Left in, the seed partition is mostly an artefact of how thoroughly a city's footpaths have been surveyed.

pedestrian itself is not here: Overture uses it for plazas and pedestrianised streets, which are genuine breaks in the urban fabric at the width a street has.

MERGE_COLUMNS module-attribute

MERGE_COLUMNS: tuple[str, ...] = ('building_surface_fraction', 'height_of_roughness_elements_m')

The seed-level features the merge compares neighbours on.

Deliberately the two cheap ones. Both come from buildings_area alone with a single overlay - no land cover, no street profile, no second compute_parameters - and between them they carry the two axes LCZ separates built types on: compactness and height band. Adding the land-cover fractions would double the cost of unit generation to refine a decision the classifier makes afterwards anyway.

PatchReport dataclass

PatchReport(n_seeds: int, n_patches: int, n_merges: int, n_isolates: int, n_below_minimum: int, n_blocked_by_max_area: int, n_seeds_split: int, n_above_maximum: int, area_above_maximum: float, min_area_m2: float, max_area_m2: float, seed_area_quantiles: dict[str, float] = dict(), patch_area_quantiles: dict[str, float] = dict())

What the merge did, for the run manifest.

A unit strategy whose output depends on a threshold has to say what that threshold produced, or two runs at different settings are indistinguishable in the archive.

n_isolates instance-attribute

n_isolates: int

Seeds adjacent to nothing - an island, or a face fully enclosed by the study boundary. They stay as their own patch at whatever size they are, since there is nothing to merge them into.

n_below_minimum instance-attribute

n_below_minimum: int

Patches still under min_area_m2 when the merge stopped. Non-zero exactly where isolates or the max_area_m2 ceiling blocked further merging, so it is the number that says whether the target was reachable rather than merely requested.

n_blocked_by_max_area instance-attribute

n_blocked_by_max_area: int

Merges that took the smallest neighbour rather than the most similar one, because every similar neighbour would have breached max_area_m2. High values mean the two thresholds are fighting each other.

n_seeds_split instance-attribute

n_seeds_split: int

Seeds that exceeded max_area_m2 and were subdivided before merging.

Zero on a barrier set that produced no oversized face, and zero whenever max_area_m2 is None. High values say the barrier network left large unbounded faces — usually water, or a periphery Overture does not cover — which is worth seeing rather than inferring from a suspiciously round patch count.

n_above_maximum instance-attribute

n_above_maximum: int

Patches larger than max_area_m2 when the merge stopped.

Non-zero only where a merge could not be avoided: split_oversized cuts every oversized seed before the merge starts, so what remains here is a unit the merge itself pushed over the line because the alternative was leaving a sliver stranded. That trade is deliberate — see n_blocked_by_max_area — and this is the count that says how often it was taken.

area_above_maximum instance-attribute

area_above_maximum: float

Total area in those patches, m². The count alone understates it badly: a handful of units can hold most of the extent.

patch_area_quantiles class-attribute instance-attribute

patch_area_quantiles: dict[str, float] = field(default_factory=dict)

p10/p50/p90 in m², before and after. The pair is the point: it says what the merge moved.

PatchUnits

PatchUnits(*, min_area_m2: float = DEFAULT_MIN_AREA_M2, max_area_m2: float | None = DEFAULT_MAX_AREA_M2, buildings: GeoDataFrame | None = None)

Enclosure seeds merged to LCZ-patch scale, satisfying SpatialUnitStrategy.

buildings is taken at construction rather than passed to generate because the protocol's signature is (bbox, barriers) and widening it for one strategy would put a building layer in the interface every strategy has to accept. Without it the merge runs on size alone, which is supported and worse — MERGE_COLUMNS is what makes a patch homogeneous rather than merely big.

The last PatchReport is kept on the instance so a caller can put it in the manifest without the strategy having to return two things and break the protocol.

Set the area floor the merge works towards, and the building layer it judges by.

min_area_m2 is a floor, not a target: merging stops when a unit reaches it and the merge that got it there overshoots, so 5 ha yields a ~10.5 ha median. max_area_m2 blocks a merge that would overshoot too far. buildings is optional and the merge runs on size alone without it — supported, and worse.

Source code in src/lczkit/units/patches.py

def __init__(
    self,
    *,
    min_area_m2: float = DEFAULT_MIN_AREA_M2,
    max_area_m2: float | None = DEFAULT_MAX_AREA_M2,
    buildings: gpd.GeoDataFrame | None = None,
) -> None:
    """Set the area floor the merge works towards, and the building layer it judges by.

    `min_area_m2` is a floor, not a target: merging stops when a unit *reaches* it and the
    merge that got it there overshoots, so 5 ha yields a ~10.5 ha median. `max_area_m2`
    blocks a merge that would overshoot too far. `buildings` is optional and the merge runs
    on size alone without it — supported, and worse.
    """
    self.min_area_m2 = min_area_m2
    self.max_area_m2 = max_area_m2
    self.buildings = buildings
    self.report: PatchReport | None = None

generate

generate(bbox: BBox, barriers: GeoDataFrame | None = None) -> GeoDataFrame

Seed enclosures over barriers and merge them to patch scale, indexed by unit_id.

Same contract as the other two strategies — bbox lon/lat, geometry-only frame back in the projected CRS — and the same requirement as EnclosureUnits that barriers be supplied, since the seeds are enclosures.

Sets self.report as a side effect. The protocol returns one frame, and a caller that wants the merge outcome in the manifest reads it off the instance afterwards rather than the interface growing a second return value for one strategy.

Source code in src/lczkit/units/patches.py

def generate(self, bbox: BBox, barriers: gpd.GeoDataFrame | None = None) -> gpd.GeoDataFrame:
    """Seed enclosures over `barriers` and merge them to patch scale, indexed by `unit_id`.

    Same contract as the other two strategies — `bbox` lon/lat, geometry-only frame back in
    the projected CRS — and the same requirement as `EnclosureUnits` that `barriers` be
    supplied, since the seeds are enclosures.

    Sets `self.report` as a side effect. The protocol returns one frame, and a caller that
    wants the merge outcome in the manifest reads it off the instance afterwards rather than
    the interface growing a second return value for one strategy.
    """
    seeds = EnclosureUnits().generate(bbox, barriers)
    patches, report = merge_to_patches(
        seeds,
        min_area_m2=self.min_area_m2,
        max_area_m2=self.max_area_m2,
        buildings=self.buildings,
    )
    self.report = report
    return patches

filter_street_barriers

filter_street_barriers(streets: GeoDataFrame, *, drop_classes: frozenset[str] = PEDESTRIAN_CLASSES, class_column: str = 'class') -> GeoDataFrame

streets without the classes that are not boundaries between LCZ patches.

A no-op returning streets unchanged when the frame has no class column, which is the honest behaviour for a VectorSource that does not supply one: silently treating every segment as keepable is what the caller would get anyway, and raising would make the class column a hard requirement of a protocol that does not promise it.

Overture's transportation/segment carries class on 87-99% of rows globally, and OvertureSource already selects it, so on the shipped source this is free.

Source code in src/lczkit/units/patches.py

def filter_street_barriers(
    streets: gpd.GeoDataFrame,
    *,
    drop_classes: frozenset[str] = PEDESTRIAN_CLASSES,
    class_column: str = "class",
) -> gpd.GeoDataFrame:
    """`streets` without the classes that are not boundaries between LCZ patches.

    A no-op returning `streets` unchanged when the frame has no class column, which is the honest
    behaviour for a `VectorSource` that does not supply one: silently treating every segment as
    keepable is what the caller would get anyway, and raising would make the class column a hard
    requirement of a protocol that does not promise it.

    Overture's `transportation/segment` carries `class` on 87-99% of rows globally, and
    `OvertureSource` already selects it, so on the shipped source this is free.
    """
    assert_projected_crs(streets, "streets")
    if class_column not in streets.columns or not drop_classes:
        return streets
    keep = ~streets[class_column].isin(drop_classes)
    return gpd.GeoDataFrame(streets.loc[keep])

seed_features

seed_features(seeds: GeoDataFrame, buildings: GeoDataFrame) -> DataFrame

MERGE_COLUMNS per seed, from buildings alone.

Building surface fraction is footprint area overlaid onto the seed, over seed area — the same definition lczkit.ucp.buildings uses, recomputed here rather than imported because this runs before units exist and compute_parameters takes units as an input.

Height is the unweighted geometric mean, matching Hr, so a seed's value here and its value in the parameter table are the same statistic. Null where the seed holds no building with a height, and the merge handles that as a missing dimension rather than as a zero — a block of unmeasured buildings is not a block of 1 m buildings.

buildings must already carry heights; pass what the cascade returned.

Source code in src/lczkit/units/patches.py

def seed_features(seeds: gpd.GeoDataFrame, buildings: gpd.GeoDataFrame) -> pd.DataFrame:
    """`MERGE_COLUMNS` per seed, from `buildings` alone.

    Building surface fraction is footprint area overlaid onto the seed, over seed area — the same
    definition `lczkit.ucp.buildings` uses, recomputed here rather than imported because this runs
    *before* units exist and `compute_parameters` takes units as an input.

    Height is the unweighted geometric mean, matching `Hr`, so a seed's value here and its value in
    the parameter table are the same statistic. Null where the seed holds no building with a height,
    and the merge handles that as a missing dimension rather than as a zero — a block of unmeasured
    buildings is not a block of 1 m buildings.

    `buildings` must already carry heights; pass what the cascade returned.
    """
    assert_projected_crs(seeds, "seeds")
    assert_projected_crs(buildings, "buildings")
    frame = pd.DataFrame(
        {column: np.nan for column in MERGE_COLUMNS},
        index=seeds.index,
    )
    if buildings.empty:
        frame["building_surface_fraction"] = 0.0
        return frame

    pieces = gpd.overlay(
        seeds.reset_index()[["unit_id", "geometry"]],
        buildings[["geometry", *(c for c in ("height",) if c in buildings.columns)]],
        how="intersection",
        keep_geom_type=True,
    )
    seed_area = seeds.geometry.area
    if pieces.empty:
        frame["building_surface_fraction"] = 0.0
        return frame

    pieces = pieces.assign(piece_area=pieces.geometry.area)
    covered = pieces.groupby("unit_id")["piece_area"].sum()
    frame["building_surface_fraction"] = (
        covered.reindex(seeds.index).fillna(0.0).div(seed_area.where(seed_area > 0)).clip(upper=1.0)
    )
    if "height" in pieces.columns:
        heights = pd.to_numeric(pieces["height"], errors="coerce")
        usable = pieces.loc[heights.gt(0).fillna(False)].assign(log_h=np.log(heights[heights > 0]))
        if not usable.empty:
            frame["height_of_roughness_elements_m"] = np.exp(
                usable.groupby("unit_id")["log_h"].mean()
            ).reindex(seeds.index)
    return frame

split_oversized

split_oversized(seeds: GeoDataFrame, max_area_m2: float | None) -> tuple[GeoDataFrame, int]

Subdivide any seed larger than max_area_m2, returning the seeds and how many were split.

max_area_m2 was a merge guard and not a ceiling, and this is what makes the name true. It refused to combine two seeds into something oversized and had no way to divide a seed that already exceeded it — and enclosure seeds routinely do, because a face bounded by nothing but the study edge is as large as the unmapped ground it covers. On a 4 555 km² Istanbul extent 807 patches exceeded the shipped 50 ha setting and held 72.7% of the total area, the largest being 1 073 km²; one 98 km² unit contained 1 310 buildings and was given a single LCZ label at a uniqueness of 0.12. Benign where the giant faces are sea. Not benign where a city's periphery is simply unsurveyed, which is the case in the cities this package exists to reach.

The cut is a regular grid of side sqrt(max_area_m2) anchored on each seed's own bounds, and it is deliberately the dullest thing that works. It needs no building layer, so it behaves the same on the unmapped hinterland that produces most oversized seeds; it is deterministic, so two runs over the same extent agree; and intersecting a polygon with a grid that covers it preserves the partition exactly — the pieces union back to the seed, and no ground is gained or lost. Pieces that a concave seed leaves disconnected are exploded, so every unit stays a single polygon.

Slivers along the cut lines are expected and are not a problem: merge_to_patches runs next and absorbs anything under min_area_m2 into its most similar neighbour, which is the same treatment enclosure slivers already get.

Piece count is bounded by the seed area over max_area_m2 — that is, by the number of units the caller asked for — so this adds no unbounded operation to the stage.

Source code in src/lczkit/units/patches.py

def split_oversized(
    seeds: gpd.GeoDataFrame, max_area_m2: float | None
) -> tuple[gpd.GeoDataFrame, int]:
    """Subdivide any seed larger than `max_area_m2`, returning the seeds and how many were split.

    **`max_area_m2` was a merge guard and not a ceiling, and this is what makes the name true.**
    It refused to *combine* two seeds into something oversized and had no way to divide a seed that
    already exceeded it — and enclosure seeds routinely do, because a face bounded by nothing but
    the study edge is as large as the unmapped ground it covers. On a 4 555 km² Istanbul extent 807
    patches exceeded the shipped 50 ha setting and held 72.7% of the total area, the largest being
    1 073 km²; one 98 km² unit contained 1 310 buildings and was given a single LCZ label at a
    uniqueness of 0.12. Benign where the giant faces are sea. Not benign where a city's periphery is
    simply unsurveyed, which is the case in the cities this package exists to reach.

    The cut is a regular grid of side `sqrt(max_area_m2)` anchored on each seed's own bounds, and it
    is deliberately the dullest thing that works. It needs no building layer, so it behaves the same
    on the unmapped hinterland that produces most oversized seeds; it is deterministic, so two runs
    over the same extent agree; and intersecting a polygon with a grid that covers it preserves the
    partition exactly — the pieces union back to the seed, and no ground is gained or lost. Pieces
    that a concave seed leaves disconnected are exploded, so every unit stays a single polygon.

    Slivers along the cut lines are expected and are not a problem: `merge_to_patches` runs next and
    absorbs anything under `min_area_m2` into its most similar neighbour, which is the same
    treatment enclosure slivers already get.

    Piece count is bounded by the seed area over `max_area_m2` — that is, by the number of units the
    caller asked for — so this adds no unbounded operation to the stage.
    """
    if max_area_m2 is None:
        return seeds, 0
    oversized = seeds.geometry.area > max_area_m2
    if not bool(oversized.any()):
        return seeds, 0

    side = float(np.sqrt(max_area_m2))
    kept = [seeds.loc[~oversized]]
    for unit_id, geometry in seeds.loc[oversized, "geometry"].items():
        pieces = _grid_pieces(geometry, side)
        kept.append(
            gpd.GeoDataFrame(
                {"unit_id": [f"{unit_id}_s{index:04d}" for index in range(len(pieces))]},
                geometry=gpd.GeoSeries(pieces, crs=seeds.crs),
            ).set_index("unit_id")
        )
    out = pd.concat(kept)
    return gpd.GeoDataFrame(out, geometry="geometry", crs=seeds.crs), int(oversized.sum())

merge_to_patches

merge_to_patches(seeds: GeoDataFrame, features: DataFrame | None = None, *, min_area_m2: float = DEFAULT_MIN_AREA_M2, max_area_m2: float | None = DEFAULT_MAX_AREA_M2, buildings: GeoDataFrame | None = None) -> tuple[GeoDataFrame, PatchReport]

Merge contiguous seeds until each reaches min_area_m2, most similar neighbour first.

Repeatedly takes the smallest surviving unit and merges it into the contiguous neighbour closest to it in features, subject to the combined area not exceeding max_area_m2. Where no neighbour satisfies that ceiling it takes the smallest neighbour instead, so the loop always makes progress and cannot stall on a unit hemmed in by large ones.

Deterministic. Units are processed smallest-area-first with ties broken on unit_id, candidate neighbours are scanned in sorted order, and equal distances break on unit_id too. The same seeds in a different row order give the same patches.

A partition in, a partition out. Only contiguous units are ever merged and the result is their union, so the covered ground and the absence of overlap are both preserved exactly from the seed partition. EnclosureUnits supplies that with clip=True.

features may be None, in which case every neighbour is equally similar and the merge runs on size alone. That is a worse unit and it is offered because it needs no building layer.

Pass buildings instead of features whenever seeds may exceed max_area_m2: oversized seeds are split first, and features are then computed on the pieces. A precomputed features frame cannot describe pieces that did not exist when it was built, so one that lacks a row for any seed the merge runs on is refused rather than read as all-missing, which _distance would score as identical to every neighbour.

Source code in src/lczkit/units/patches.py

def merge_to_patches(
    seeds: gpd.GeoDataFrame,
    features: pd.DataFrame | None = None,
    *,
    min_area_m2: float = DEFAULT_MIN_AREA_M2,
    max_area_m2: float | None = DEFAULT_MAX_AREA_M2,
    buildings: gpd.GeoDataFrame | None = None,
) -> tuple[gpd.GeoDataFrame, PatchReport]:
    """Merge contiguous `seeds` until each reaches `min_area_m2`, most similar neighbour first.

    Repeatedly takes the **smallest** surviving unit and merges it into the contiguous neighbour
    closest to it in `features`, subject to the combined area not exceeding `max_area_m2`. Where no
    neighbour satisfies that ceiling it takes the smallest neighbour instead, so the loop always
    makes progress and cannot stall on a unit hemmed in by large ones.

    **Deterministic.** Units are processed smallest-area-first with ties broken on `unit_id`,
    candidate neighbours are scanned in sorted order, and equal distances break on `unit_id` too.
    The same seeds in a different row order give the same patches.

    **A partition in, a partition out.** Only contiguous units are ever merged and the result is
    their union, so the covered ground and the absence of overlap are both preserved exactly from
    the seed partition. `EnclosureUnits` supplies that with `clip=True`.

    `features` may be `None`, in which case every neighbour is equally similar and the merge runs on
    size alone. That is a worse unit and it is offered because it needs no building layer.

    Pass `buildings` instead of `features` whenever seeds may exceed `max_area_m2`: oversized seeds
    are split first, and features are then computed on the pieces. A precomputed `features` frame
    cannot describe pieces that did not exist when it was built, so one that lacks a row for any
    seed the merge runs on is refused rather than read as all-missing, which `_distance` would
    score as identical to every neighbour.
    """
    assert_projected_crs(seeds, "seeds")
    if seeds.index.name != "unit_id":
        raise ValueError("seeds must be indexed by unit_id")
    if not seeds.index.is_unique:
        raise ValueError("seeds index must be unique")
    if min_area_m2 <= 0:
        raise ValueError(f"min_area_m2 must be positive, got {min_area_m2}")
    if max_area_m2 is not None and max_area_m2 < min_area_m2:
        raise ValueError(
            f"max_area_m2 ({max_area_m2}) must be at least min_area_m2 ({min_area_m2}); "
            "a ceiling below the target would block every merge"
        )
    if features is not None and buildings is not None:
        raise ValueError("pass either features or buildings, not both")

    # Before anything else, and before the seed quantiles are taken, so those describe the seeds
    # the merge actually ran on rather than the faces the barrier set happened to produce.
    seeds, n_split = split_oversized(seeds, max_area_m2)
    if buildings is not None:
        features = seed_features(seeds, buildings)
    elif features is not None:
        missing = seeds.index.difference(features.index)
        if len(missing):
            raise ValueError(
                f"features has no row for {len(missing)} seed(s), e.g. {missing[0]!r}: they were "
                "split from an oversized seed after the features were computed. Pass buildings= "
                "so features are computed on the split seeds."
            )

    seed_area = seeds.geometry.area
    quantiles = {
        q: float(seed_area.quantile(v)) for q, v in (("p10", 0.1), ("p50", 0.5), ("p90", 0.9))
    }
    if seeds.empty:
        return seeds, PatchReport(
            n_seeds=0,
            n_patches=0,
            n_merges=0,
            n_isolates=0,
            n_below_minimum=0,
            n_blocked_by_max_area=0,
            n_seeds_split=n_split,
            n_above_maximum=0,
            area_above_maximum=0.0,
            min_area_m2=min_area_m2,
            max_area_m2=max_area_m2 or float("inf"),
        )

    graph = Graph.build_contiguity(seeds, rook=False)
    neighbours: dict[str, set[str]] = {
        str(focal): {str(other) for other in tuple(others)}
        for focal, others in graph.neighbors.items()
    }
    n_isolates = sum(1 for others in neighbours.values() if not others)

    standardised = (
        _standardise(features.reindex(seeds.index)[list(MERGE_COLUMNS)])
        if features is not None
        else np.zeros((len(seeds), 0))
    )
    position = {str(uid): i for i, uid in enumerate(seeds.index)}
    vectors = {str(uid): standardised[i].copy() for uid, i in position.items()}
    areas = {str(uid): float(value) for uid, value in seed_area.items()}
    members: dict[str, list[str]] = {str(uid): [str(uid)] for uid in seeds.index}

    alive = set(areas)
    queue = [(area, uid) for uid, area in areas.items() if area < min_area_m2]
    heapq.heapify(queue)
    n_merges = 0
    n_blocked = 0

    while queue:
        area, uid = heapq.heappop(queue)
        # Lazy deletion: an entry whose area no longer matches describes a unit that has since
        # grown, and one whose id is gone describes a unit that was absorbed.
        if uid not in alive or areas[uid] != area or area >= min_area_m2:
            continue
        candidates = sorted(neighbours[uid] & alive)
        if not candidates:
            continue

        fits = [
            c for c in candidates if max_area_m2 is None or areas[uid] + areas[c] <= max_area_m2
        ]
        if fits:
            best = min(fits, key=lambda c: (_distance(vectors[uid], vectors[c]), c))
        else:
            # Every neighbour breaches the ceiling. Take the smallest so the loop still advances —
            # refusing would leave a sliver permanently, which is the state this exists to remove —
            # and count it, because a high rate means the two thresholds are fighting each other.
            n_blocked += 1
            best = min(candidates, key=lambda c: (areas[c], c))
        _absorb(uid, best, areas, members, neighbours, vectors, alive)
        n_merges += 1
        if areas[best] < min_area_m2:
            heapq.heappush(queue, (areas[best], best))

    return _assemble(
        seeds,
        members,
        alive,
        quantiles,
        min_area_m2,
        max_area_m2,
        n_merges,
        n_isolates,
        n_blocked,
        n_split,
    )

Aggregation

lczkit.units.aggregate

aggregate(): move attribute columns between two unit systems via an area overlay.

aggregate

aggregate(from_units: GeoDataFrame, to_units: GeoDataFrame, method: AggregateMethod) -> GeoDataFrame

Move every non-geometry column of from_units onto to_units' polygons.

Both must be indexed by unit_id, in the same projected CRS — typically both are outputs of a SpatialUnitStrategy.generate() call, with from_units carrying attribute columns already joined onto it elsewhere in the pipeline. This function is a geometric overlay plus a column-wise reduction; it does not compute any parameter itself.

"majority": every column takes the value from whichever from_units polygon overlaps a given to_units polygon by the largest area. Works for any column dtype. "area_weighted": every numeric column becomes the area-weighted mean across overlapping from_units polygons; non-numeric columns are dropped (there is no defined mean of a string) rather than silently guessing a reducer for them.

to_units rows with no overlapping from_units polygon get null attribute values.

Every result carries an aggregate_coverage column: the share of the target polygon that overlapping source polygons actually cover. "area_weighted" normalises by the summed overlap area rather than by the target's own area, so a target one tenth covered reports the mean over that tenth and is otherwise indistinguishable from a fully measured one. Reporting coverage beside the value is what makes the two distinguishable; the normalisation is left alone because changing it would move every arm-B projection.

Source code in src/lczkit/units/aggregate.py

def aggregate(
    from_units: gpd.GeoDataFrame,
    to_units: gpd.GeoDataFrame,
    method: AggregateMethod,
) -> gpd.GeoDataFrame:
    """Move every non-geometry column of `from_units` onto `to_units`' polygons.

    Both must be indexed by `unit_id`, in the same projected CRS — typically both are outputs
    of a `SpatialUnitStrategy.generate()` call, with `from_units` carrying attribute columns
    already joined onto it elsewhere in the pipeline. This function is a geometric overlay plus
    a column-wise reduction; it does not compute any parameter itself.

    `"majority"`: every column takes the value from whichever `from_units` polygon overlaps a
    given `to_units` polygon by the largest area. Works for any column dtype.
    `"area_weighted"`: every *numeric* column becomes the area-weighted mean across overlapping
    `from_units` polygons; non-numeric columns are dropped (there is no defined mean of a
    string) rather than silently guessing a reducer for them.

    `to_units` rows with no overlapping `from_units` polygon get null attribute values.

    Every result carries an **`aggregate_coverage`** column: the share of the target polygon that
    overlapping source polygons actually cover. `"area_weighted"` normalises by the summed overlap
    area rather than by the target's own area, so a target one tenth covered reports the mean over
    that tenth and is otherwise indistinguishable from a fully measured one. Reporting coverage
    beside the value is what makes the two distinguishable; the normalisation is left alone because
    changing it would move every arm-B projection.
    """
    assert_projected_crs(from_units, "from_units")
    assert_projected_crs(to_units, "to_units")
    if from_units.crs != to_units.crs:
        raise ValueError(f"from_units.crs ({from_units.crs}) != to_units.crs ({to_units.crs})")
    if from_units.index.name != "unit_id" or to_units.index.name != "unit_id":
        raise ValueError("both from_units and to_units must be indexed by unit_id")

    value_cols = [c for c in from_units.columns if c != "geometry"]

    left = to_units[["geometry"]].reset_index().rename(columns={"unit_id": "to_id"})
    right = (
        from_units[[*value_cols, "geometry"]].reset_index().rename(columns={"unit_id": "from_id"})
    )
    pieces = gpd.overlay(left, right, how="intersection")
    pieces["overlap_area"] = pieces.geometry.area
    pieces = pieces[pieces["overlap_area"] > 0]
    if pieces.empty:
        raise ValueError("from_units and to_units do not overlap with positive area")

    agg: pd.DataFrame
    if method == "majority":
        idx = pieces.groupby("to_id")["overlap_area"].idxmax()
        agg = pd.DataFrame(pieces.loc[idx].set_index("to_id")[value_cols])
    elif method == "area_weighted":
        numeric_cols = [c for c in value_cols if pd.api.types.is_numeric_dtype(from_units[c])]
        if not numeric_cols:
            raise ValueError("area_weighted requires at least one numeric column in from_units")
        weighted = pieces[numeric_cols].multiply(pieces["overlap_area"], axis=0)
        weighted["to_id"] = pieces["to_id"]
        sums = weighted.groupby("to_id").sum()
        area_sums = pieces.groupby("to_id")["overlap_area"].sum()
        agg = sums[numeric_cols].div(area_sums, axis=0)
    else:
        raise ValueError(f"Unknown aggregate method: {method!r}")

    target_area = to_units.geometry.area
    coverage = (
        pieces.groupby("to_id")["overlap_area"]
        .sum()
        .reindex(to_units.index)
        .div(target_area.where(target_area > 0))
    )

    joined = to_units[["geometry"]].join(agg, how="left")
    joined["aggregate_coverage"] = coverage
    result = gpd.GeoDataFrame(joined, geometry="geometry", crs=to_units.crs)
    result.index.name = "unit_id"
    return result

Overlay

The one definition of "intersect a layer with the units and measure the pieces" — the operation every per-unit area statistic in the package is built from, and which five near-identical private helpers each carried a copy of.

lczkit.units.overlay

Overlaying a layer against the units: the one definition every per-unit area statistic uses.

Intersect once with unit_pieces, then select, sum and share over the pieces. A footprint straddling a unit boundary contributes its share to each side, the same rule the height provenance uses, so every fraction built here shares a denominator with building_surface_fraction.

covered_fraction(dissolve=True) clips first and dissolves per unit. A whole-layer union_all is superlinear and raises GEOSException: side location conflict on real Overture land use even after make_valid; the union of clipped pieces within a unit equals the clip of the global union, so the bounded form is exact.

PIECE_AREA module-attribute

PIECE_AREA = 'piece_area'

Column unit_pieces measures each intersection into. Named rather than recomputed downstream so two consumers of the same pieces cannot disagree about what area means.

unit_pieces

unit_pieces(units: GeoDataFrame, layer: GeoDataFrame | GeoSeries, *, columns: Sequence[str] = (), keep_geom_type: bool = True) -> GeoDataFrame

layer intersected with units, one row per (unit, feature) pair.

Carries unit_id, the piece geometry, piece_area, and whichever of columns the layer has. Absent columns are skipped rather than raising, so a caller can ask for height and subtype without first checking which of them a hand-assembled layer happens to carry.

The layer is reset positionally before the overlay. The building layers carry no uniqueness guarantee and .loc[an_index] over a duplicated one silently returns extra rows — a wrong number rather than an error — so nothing here selects by index.

Returns an empty frame with the right columns when either side is empty, so callers branch on .empty and never on None.

Source code in src/lczkit/units/overlay.py

def unit_pieces(
    units: gpd.GeoDataFrame,
    layer: gpd.GeoDataFrame | gpd.GeoSeries,
    *,
    columns: Sequence[str] = (),
    keep_geom_type: bool = True,
) -> gpd.GeoDataFrame:
    """`layer` intersected with `units`, one row per (unit, feature) pair.

    Carries `unit_id`, the piece geometry, `piece_area`, and whichever of `columns` the layer has.
    Absent columns are skipped rather than raising, so a caller can ask for `height` and `subtype`
    without first checking which of them a hand-assembled layer happens to carry.

    The layer is reset positionally before the overlay. The building layers carry no uniqueness
    guarantee and `.loc[an_index]` over a duplicated one silently returns extra rows — a wrong
    number rather than an error — so nothing here selects by index.

    Returns an empty frame with the right columns when either side is empty, so callers branch on
    `.empty` and never on `None`.
    """
    # Checked here as well as in the blocks that call it. Everything below groups by `unit_id`, so
    # a frame indexed under another name fails several lines later with a `KeyError` naming a
    # column the caller never mentioned; a geographic CRS fails not at all, and silently reports
    # areas in square degrees.
    check_units(units)
    if isinstance(layer, gpd.GeoSeries):
        wanted: list[str] = []
        covering = gpd.GeoDataFrame(geometry=layer.reset_index(drop=True), crs=units.crs)
    else:
        wanted = [column for column in columns if column in layer.columns]
        name = layer.geometry.name
        covering = gpd.GeoDataFrame(
            layer[[*wanted, name]].reset_index(drop=True), geometry=name, crs=layer.crs
        )
        # `gpd.overlay` joins on the *active* geometry but the result keeps the left frame's
        # column name, so a right-hand layer whose geometry column is called something else still
        # works — renaming unconditionally does not, because geopandas refuses to rename a column
        # to the name it already has.
        if name != "geometry":
            covering = covering.rename_geometry("geometry")

    if units.empty or covering.empty:
        return _empty_pieces(units, wanted)

    pieces = gpd.overlay(
        units[["geometry"]].reset_index(),
        covering,
        how="intersection",
        keep_geom_type=keep_geom_type,
    )
    if pieces.empty:
        return _empty_pieces(units, wanted)
    return pieces.assign(**{PIECE_AREA: pieces.geometry.area})

area_in_units

area_in_units(units: GeoDataFrame, pieces: GeoDataFrame) -> Series

Total piece area per unit, zero where nothing reached the unit.

Zero rather than null: "nothing of this layer is here" is a measurement, unlike a land-cover fraction over ground the raster never covered. Callers that need the undefined case — a share of a unit holding no buildings — mask it themselves, where the reason is visible.

Source code in src/lczkit/units/overlay.py

def area_in_units(units: gpd.GeoDataFrame, pieces: gpd.GeoDataFrame) -> pd.Series:
    """Total piece area per unit, zero where nothing reached the unit.

    Zero rather than null: "nothing of this layer is here" is a measurement, unlike a land-cover
    fraction over ground the raster never covered. Callers that need the undefined case — a share
    of a unit holding no buildings — mask it themselves, where the reason is visible.
    """
    zero = pd.Series(0.0, index=units.index, dtype="float64")
    if pieces.empty:
        return zero
    summed = pieces.groupby("unit_id")[PIECE_AREA].sum()
    return summed.reindex(units.index).fillna(0.0)

covered_fraction

covered_fraction(units: GeoDataFrame, pieces: GeoDataFrame, *, dissolve: bool) -> Series

Share of each unit's area covered by pieces.

dissolve unions the pieces within each unit first, so ground under two overlapping features counts once. Required wherever the source layer has no overlap resolution — cleaning.land_use applies make_valid and nothing else, and Milan's parcels sum to 106.6% of its bbox — and wrong to pay for where it does: trim_overlaps has already made buildings_area disjoint.

Explicit rather than defaulted, because both mistakes are silent. Without it a fraction can exceed 1.0; with it unnecessarily, a run pays for a union nobody needed.

Source code in src/lczkit/units/overlay.py

def covered_fraction(
    units: gpd.GeoDataFrame, pieces: gpd.GeoDataFrame, *, dissolve: bool
) -> pd.Series:
    """Share of each unit's area covered by `pieces`.

    `dissolve` unions the pieces within each unit first, so ground under two overlapping features
    counts once. Required wherever the source layer has no overlap resolution — `cleaning.land_use`
    applies `make_valid` and nothing else, and Milan's parcels sum to 106.6% of its bbox — and
    wrong to pay for where it does: `trim_overlaps` has already made `buildings_area` disjoint.

    Explicit rather than defaulted, because both mistakes are silent. Without it a fraction can
    exceed 1.0; with it unnecessarily, a run pays for a union nobody needed.
    """
    unit_area = units.geometry.area
    if pieces.empty:
        return pd.Series(0.0, index=units.index, dtype="float64")
    covered = (
        pieces.dissolve(by="unit_id").geometry.area
        if dissolve
        else pieces.groupby("unit_id")[PIECE_AREA].sum()
    )
    return covered.reindex(units.index).fillna(0.0).div(unit_area.where(unit_area > 0))

share_of

share_of(numerator: Series, denominator: Series) -> Series

numerator / denominator, null where the denominator is zero.

The shape every "of what is built here, how much is X" column takes. Null rather than zero: a share of nothing is undefined, and reporting 0.0 tells a downstream rule that a cell definitely is not X rather than that there was nothing to judge.

Source code in src/lczkit/units/overlay.py

def share_of(numerator: pd.Series, denominator: pd.Series) -> pd.Series:
    """`numerator / denominator`, null where the denominator is zero.

    The shape every "of what is built here, how much is X" column takes. Null rather than zero: a
    share of nothing is undefined, and reporting 0.0 tells a downstream rule that a cell definitely
    is not X rather than that there was nothing to judge.
    """
    return numerator.div(denominator.where(denominator > 0))

Tessellation

A fourth strategy, built for lczkit.morphometrics rather than for classification: one polygon per building (an "enclosed tessellation cell"), not a partition of the bounding box — cells with no parent building are dropped, following Majer & Fleischmann (2026). Deliberately not reachable via UnitsConfig.strategy; see the module docstring for why.

lczkit.units.tessellation

TessellationUnits: building-level enclosed tessellation cells (ETCs).

The units lczkit.morphometrics computes the 2D attributes of Majer & Fleischmann (2026) on, built with momepy.enclosed_tessellation inside the street-bounded enclosures. Not offered as a classification unit through UnitsConfig.strategy.

Tessellation does not partition bbox, unlike GridUnits/EnclosureUnits. momepy.enclosed_tessellation assigns a negative index to any tessellation cell with no parent building (an enclosure interior, or a sliver of one, that never reached the input geometry) — those rows are dropped here, following Majer & Fleischmann (2026) §Supplementary A, which excludes cells without a parent building from the morphometric computation entirely. So a road median or an unbuilt block interior simply has no unit_id, which is correct for a strategy whose whole purpose is per-building 2D shape description, not spatial coverage.

DEFAULT_THRESHOLD module-attribute

DEFAULT_THRESHOLD = 0.05

momepy's own enclosed_tessellation defaults, restated as module constants rather than left implicit in a call site, so a run's manifest can record what actually ran.

TessellationReport dataclass

TessellationReport(n_enclosures: int, n_buildings_in: int, n_etc: int, n_excluded_no_parent_building: int, etc_area_quantiles: dict[str, float] = dict())

What tessellation produced, for the run manifest.

A unit strategy that drops rows (unlike GridUnits/EnclosureUnits, which partition the bbox exactly) has to say how many and why, or a city with an unusually sparse building layer is indistinguishable from a tessellation bug.

n_excluded_no_parent_building instance-attribute

n_excluded_no_parent_building: int

Tessellation cells momepy.enclosed_tessellation returned with a negative index — an enclosure interior (or a sliver of one) with no matching building. Excluded per Majer & Fleischmann (2026)'s own rule: ETCs without a parent building carry no morphometric meaning and are dropped before any attribute is computed.

etc_area_quantiles class-attribute instance-attribute

etc_area_quantiles: dict[str, float] = field(default_factory=dict)

p10/p50/p90 ETC area in m², after exclusion.

TessellationUnits

TessellationUnits(*, buildings: GeoDataFrame, shrink: float = DEFAULT_SHRINK, segment: float = DEFAULT_SEGMENT, threshold: float | None = DEFAULT_THRESHOLD, n_jobs: int = -1)

Building-level enclosed tessellation cells, satisfying SpatialUnitStrategy.

buildings is taken at construction rather than passed to generate, following the PatchUnits precedent (lczkit.units.patches) — enclosed_tessellation needs building footprints, and the protocol's generate(bbox, barriers) signature has no room for a second geometry layer without widening every strategy's interface for one caller.

The last TessellationReport is kept on the instance, the same way PatchUnits.report is, so a caller can put it in the manifest without the protocol growing a second return value.

Set the building layer tessellation is generated from, and momepy's own tuning knobs.

buildings should be buildings_area from lczkit.cleaning.pipeline.clean_vectors — the area-preserving layer, matching every other area statistic in this package. shrink/segment/threshold are passed straight through to momepy.enclosed_tessellation; the defaults are momepy's own.

n_jobs matters far more than its momepy default suggests, and it is the one knob a caller running its own process pool must set. momepy.enclosed_tessellation defaults to n_jobs=-1 and starts a joblib/loky pool over every core on the host — right for a single run that owns the machine, and catastrophic one level down. A multi-city driver with 8 worker processes measured 2 923 processes and ~70 cores against a 16-core budget, on a 256-core shared node whose load average reached 428. None of the OMP_NUM_THREADS/MKL/OpenBLAS pinning touches it: those cap threads inside a process, and this is a nested pool of processes. Pass n_jobs=1 when calling from inside a pool of your own; the default is left at momepy's so no existing run changes.

Source code in src/lczkit/units/tessellation.py

def __init__(
    self,
    *,
    buildings: gpd.GeoDataFrame,
    shrink: float = DEFAULT_SHRINK,
    segment: float = DEFAULT_SEGMENT,
    threshold: float | None = DEFAULT_THRESHOLD,
    n_jobs: int = -1,
) -> None:
    """Set the building layer tessellation is generated from, and momepy's own tuning knobs.

    `buildings` should be `buildings_area` from `lczkit.cleaning.pipeline.clean_vectors` —
    the area-preserving layer, matching every other area statistic in this package.
    `shrink`/`segment`/`threshold` are passed straight through to
    `momepy.enclosed_tessellation`; the defaults are momepy's own.

    **`n_jobs` matters far more than its momepy default suggests, and it is the one knob a
    caller running its own process pool must set.** `momepy.enclosed_tessellation` defaults
    to `n_jobs=-1` and starts a `joblib`/`loky` pool over *every core on the host* — right
    for a single run that owns the machine, and catastrophic one level down. A multi-city
    driver with 8 worker processes measured **2 923 processes and ~70 cores against a
    16-core budget**, on a 256-core shared node whose load average reached 428. None of the
    `OMP_NUM_THREADS`/`MKL`/`OpenBLAS` pinning touches it: those cap *threads* inside a
    process, and this is a nested pool of *processes*. Pass `n_jobs=1` when calling from
    inside a pool of your own; the default is left at momepy's so no existing run changes.
    """
    assert_projected_crs(buildings, "buildings")
    self.buildings = buildings
    self.shrink = shrink
    self.segment = segment
    self.n_jobs = n_jobs
    self.threshold = threshold
    self.report: TessellationReport | None = None

generate

generate(bbox: BBox, barriers: GeoDataFrame | None = None) -> GeoDataFrame

Enclosed tessellation over barriers, restricted to cells with a parent building.

Same contract as the other strategies: bbox lon/lat, the returned frame in the projected CRS, indexed by unit_id. barriers is required, as for EnclosureUnits — tessellation is generated within enclosures, so there is nothing to tessellate without them.

Extra columns enclosure_index and parent_building_id are carried alongside geometry, matching the precedent PatchUnits and EnclosureUnits set of returning more than a bare geometry frame where it is useful downstream — here, for joining ETC-level morphometrics back to the building that owns each cell.

Sets self.report as a side effect.

Source code in src/lczkit/units/tessellation.py

def generate(self, bbox: BBox, barriers: gpd.GeoDataFrame | None = None) -> gpd.GeoDataFrame:
    """Enclosed tessellation over `barriers`, restricted to cells with a parent building.

    Same contract as the other strategies: `bbox` lon/lat, the returned frame in the
    projected CRS, indexed by `unit_id`. `barriers` is required, as for `EnclosureUnits` —
    tessellation is generated *within* enclosures, so there is nothing to tessellate without
    them.

    Extra columns `enclosure_index` and `parent_building_id` are carried alongside geometry,
    matching the precedent `PatchUnits` and `EnclosureUnits` set of returning more than a
    bare geometry frame where it is useful downstream — here, for joining ETC-level
    morphometrics back to the building that owns each cell.

    Sets `self.report` as a side effect.
    """
    enclosures = EnclosureUnits().generate(bbox, barriers)

    buildings = usable_buildings(self.buildings)
    source_ids = building_ids(buildings)
    # Positional reset: `enclosed_tessellation` requires a unique non-negative integer index,
    # and its result is indexed by that same integer — `source_ids` stays aligned to it by
    # construction, since both are built from `buildings` in the same row order.
    positioned = buildings[["geometry"]].reset_index(drop=True)

    # momepy names its enclosure-linking output column after `enclosures`'s own index name —
    # `EnclosureUnits` indexes by `unit_id`, so that column would otherwise be called
    # `unit_id` too and collide with the one this method sets on its own result below.
    # Renaming the axis first keeps the real enclosure `unit_id` values (e.g. "enclosure_12")
    # while giving the column an unambiguous name.
    tess = momepy.enclosed_tessellation(
        positioned,
        enclosures.geometry.rename_axis("enclosure_index"),
        shrink=self.shrink,
        segment=self.segment,
        threshold=self.threshold,
        n_jobs=self.n_jobs,
    )
    matched = tess.loc[tess.index >= 0]
    n_excluded = len(tess) - len(matched)

    if matched.index.duplicated().any():
        # momepy's own documented edge case: a cell can split into a MultiPolygon during
        # shrinking, which can surface as more than one row sharing an index. Dissolving
        # keeps `unit_id` unique without discarding area — the parts union back to one cell.
        matched = matched.dissolve(by=matched.index)

    parent_ids = source_ids.loc[matched.index]
    result = gpd.GeoDataFrame(
        {
            "unit_id": "etc_" + parent_ids.to_numpy(),
            "enclosure_index": matched["enclosure_index"].to_numpy(),
            "parent_building_id": parent_ids.to_numpy(),
        },
        geometry=matched.geometry.to_numpy(),
        crs=enclosures.crs,
    ).set_index("unit_id")

    area = result.geometry.area
    self.report = TessellationReport(
        n_enclosures=int(len(enclosures)),
        n_buildings_in=int(len(buildings)),
        n_etc=int(len(result)),
        n_excluded_no_parent_building=int(n_excluded),
        etc_area_quantiles={
            "p10": float(area.quantile(0.1)) if len(area) else 0.0,
            "p50": float(area.quantile(0.5)) if len(area) else 0.0,
            "p90": float(area.quantile(0.9)) if len(area) else 0.0,
        },
    )
    return result

building_ids

building_ids(buildings: GeoDataFrame) -> Series

One id string per row of buildings, aligned to its row order, guaranteed unique.

Prefers building_id (lczkit.cleaning.buildings.BUILDING_ID), which is documented as unique on buildings_area — the layer this strategy is meant to be given — then Overture's own id, then lczkit's feature_id. The latter two are not trustworthy on their own: feature_id is stamped before cleaning explodes multipolygons specifically so height statistics can collapse split parts back together, which means the parts share it, and a raw Overture id behaves the same way after any explode. A duplicate unit_id would silently drop a cell at set_index, so whichever column is used is disambiguated with a per-group counter rather than trusted at face value; a building with none of these columns falls back to its row position, which is always unique.

This is the same id scheme TessellationUnits.generate builds unit_id from — exported so lczkit.morphometrics can reindex the raw building layer onto the ETC unit_id axis without re-deriving the scheme, which is what makes a building-level momepy result (e.g. momepy.circular_compactness(buildings)) directly joinable to an ETC-level one.

Source code in src/lczkit/units/tessellation.py

def building_ids(buildings: gpd.GeoDataFrame) -> pd.Series:
    """One id string per row of `buildings`, aligned to its row order, guaranteed unique.

    Prefers `building_id` (`lczkit.cleaning.buildings.BUILDING_ID`), which is documented as
    unique on `buildings_area` — the layer this strategy is meant to be given — then Overture's
    own `id`, then lczkit's `feature_id`. The latter two are **not** trustworthy on their own:
    `feature_id` is stamped before cleaning explodes multipolygons specifically so height
    statistics can collapse split parts back together, which means the parts share it, and a raw
    Overture `id` behaves the same way after any explode. A duplicate `unit_id` would silently
    drop a cell at `set_index`, so whichever column is used is disambiguated with a per-group
    counter rather than trusted at face value; a building with none of these columns falls back
    to its row position, which is always unique.

    This is the same id scheme `TessellationUnits.generate` builds `unit_id` from — exported so
    `lczkit.morphometrics` can reindex the raw building layer onto the ETC `unit_id` axis without
    re-deriving the scheme, which is what makes a building-level momepy result (e.g.
    `momepy.circular_compactness(buildings)`) directly joinable to an ETC-level one.
    """
    for column in ("building_id", "id", "feature_id"):
        if column in buildings.columns:
            base = buildings[column].astype(str)
            break
    else:
        base = buildings.index.to_series().astype(str)
    base = pd.Series(base.to_numpy())
    order = base.groupby(base).cumcount()
    return base.where(order == 0, base + "_" + order.astype(str))

usable_buildings

usable_buildings(buildings: GeoDataFrame) -> GeoDataFrame

buildings restricted to the rows the morphometrics stage can actually compute on.

One definition, used by both TessellationUnits.generate and buildings_for_etc, because the two must agree exactly. generate decides which buildings get a tessellation cell and buildings_for_etc reindexes buildings onto those cells, so a row either passes both filters or neither. Filtering in only one of them raises KeyError: [...] not in index from the reindex — met on real Overture data over Salvador, when the areal filter below was added to buildings_for_etc alone.

Areal geometries only, beyond the ordinary null/empty check. libpysal.graph.Graph.build_contiguity refuses anything outside {LineString, Polygon, MultiLineString, MultiPolygon} outright, so one stray point or GeometryCollection surviving cleaning takes the whole stage down with ValueError: This Graph type is only well-defined for geom_types: .... Dropped here rather than guarded at each graph call site, since every building-level metric in lczkit.morphometrics reads this frame and none of them is defined on a point either.

Source code in src/lczkit/units/tessellation.py

def usable_buildings(buildings: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """`buildings` restricted to the rows the morphometrics stage can actually compute on.

    **One definition, used by both `TessellationUnits.generate` and `buildings_for_etc`, because
    the two must agree exactly.** `generate` decides which buildings get a tessellation cell and
    `buildings_for_etc` reindexes buildings onto those cells, so a row either passes both filters
    or neither. Filtering in only one of them raises
    `KeyError: [...] not in index` from the reindex — met on real Overture data over Salvador,
    when the areal filter below was added to `buildings_for_etc` alone.

    Areal geometries only, beyond the ordinary null/empty check.
    `libpysal.graph.Graph.build_contiguity` refuses anything outside
    {LineString, Polygon, MultiLineString, MultiPolygon} outright, so one stray point or
    GeometryCollection surviving cleaning takes the whole stage down with `ValueError: This Graph
    type is only well-defined for geom_types: ...`. Dropped here rather than guarded at each graph
    call site, since every building-level metric in `lczkit.morphometrics` reads this frame and
    none of them is defined on a point either.
    """
    geometry = buildings.geometry
    areal = geometry.geom_type.isin(("Polygon", "MultiPolygon"))
    return buildings.loc[geometry.notna() & ~geometry.is_empty & areal]

buildings_for_etc

buildings_for_etc(buildings: GeoDataFrame, etc: GeoDataFrame) -> GeoDataFrame

buildings, filtered to those with a matching ETC and reindexed onto etc's unit_id.

The frame every building-level momepy call in lczkit.morphometrics should be given: its index lines up with etc.index row for row, so a Series computed from one is directly alignable with one computed from the other — no separate join step, and no risk of the two silently drifting apart if a caller filtered one but not the other.

Source code in src/lczkit/units/tessellation.py

def buildings_for_etc(buildings: gpd.GeoDataFrame, etc: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """`buildings`, filtered to those with a matching ETC and reindexed onto `etc`'s `unit_id`.

    The frame every building-level momepy call in `lczkit.morphometrics` should be given: its
    index lines up with `etc.index` row for row, so a `Series` computed from one is directly
    alignable with one computed from the other — no separate join step, and no risk of the two
    silently drifting apart if a caller filtered one but not the other.
    """
    assert_projected_crs(buildings, "buildings")
    assert_projected_crs(etc, "etc")
    valid = usable_buildings(buildings)
    keyed = valid.set_axis(pd.Index("etc_" + building_ids(valid).to_numpy(), name="unit_id"))
    # `etc.index` may hold repeats where `TessellationUnits` disambiguated a duplicate parent id
    # (`etc_<id>_1`, ...); `keyed`'s own duplicate suffixes were built the same way from the same
    # `building_ids`, so the two line up without a fuzzy join.
    return keyed.loc[etc.index]