Cleaning

Vector cleaning pipeline: building cleanup, street simplification, cross-layer topology.

Repairing and simplifying the raw Overture layers. Everything downstream depends on this, and the building surface fraction it produces carries the largest single weight in the classification.

Cleaning produces two building layers, not one. This is the most important structural rule in the package:

• buildings_topo — planar and non-overlapping, meaning no two footprints share any area. Feeds enclosure generation and anything needing a valid partition. Destructive operations are permitted.
• buildings_area — area-preserving. Feeds building surface fraction and every area statistic. Only validity fixes, multipolygon explosion, non-polygon removal, implausibly-large-footprint removal and genuine duplicate removal.

The split exists because building surface fraction carries roughly half the classification metric. Cleaning a single layer for topology measured 23.5% of Berlin's footprint area away before BSF was computed, which was worth 9.1 points of agreement. The cleaning report therefore records total footprint area in and out of every operation, not just feature counts.

lczkit.cleaning.pipeline

Top-level orchestration for vector cleaning.

Fetch raw vectors from a VectorSource and run the full cleaning pipeline against them.

CleanedVectors dataclass

CleanedVectors(buildings_area: GeoDataFrame, buildings_topo: GeoDataFrame, streets: GeoDataFrame, waterlines: GeoDataFrame, waterbodies: GeoDataFrame, land_use: GeoDataFrame, report: CleaningReport, crs: CRS)

The cleaned output of clean_vectors().

Holds live GeoDataFrames for in-process use. It is never itself serialised; lczkit.output writes the run's files.

There is no plain buildings attribute, deliberately. Which of the two layers a caller wants is never obvious from the name, and the single ambiguous layer that used to be here is what fed 23.5% of Berlin's footprint area into a statistic that needed all of it.

buildings_area instance-attribute

buildings_area: GeoDataFrame

Area-preserving. Building surface fraction, Hr, building count, mean building area and industrial_fraction all read this, and the height cascade runs on it.

buildings_topo instance-attribute

buildings_topo: GeoDataFrame

Planar and non-overlapping. The neatnet exclusion mask and momepy.street_profile read this; it has been through the road-buffer rule, so its facades stand outside the roadway.

tile_fingerprint

tile_fingerprint(config: CleaningConfig) -> str

Short hash of everything that changes a cached tile's contents.

Tile keys are aligned to the CRS origin and so are shared between runs over different bboxes — which is the point of the cache, and also why the key alone is not enough to identify a result. Anything that would make the same tile simplify differently has to move the fingerprint, or a later run silently reads an earlier run's answer.

The whole of CleaningConfig goes in, not only the tiling fields: the exclusion mask is buildings_topo, so every building threshold reaches a tile's result too. neatnet's version is included because its output is its algorithm. The pinned artifact threshold is not here — it depends on the full extent rather than on config, so simplify_streets_tiled folds it in once it has resolved it.

Source code in src/lczkit/cleaning/pipeline.py

def tile_fingerprint(config: CleaningConfig) -> str:
    """Short hash of everything that changes a cached tile's contents.

    Tile keys are aligned to the CRS origin and so are shared between runs over different
    bboxes — which is the point of the cache, and also why the key alone is not enough to
    identify a result. Anything that would make the same tile simplify differently has to move
    the fingerprint, or a later run silently reads an earlier run's answer.

    The whole of `CleaningConfig` goes in, not only the tiling fields: the exclusion mask is
    `buildings_topo`, so every building threshold reaches a tile's result too. `neatnet`'s
    version is included because its output is its algorithm. The pinned artifact threshold is
    *not* here — it depends on the full extent rather than on config, so `simplify_streets_tiled`
    folds it in once it has resolved it.
    """
    payload = {
        "cleaning": config.model_dump(mode="json"),
        "neatnet": version("neatnet"),
    }
    digest = hashlib.sha256(json.dumps(payload, sort_keys=True).encode()).hexdigest()
    return digest[:12]

reproject_to_local_utm

reproject_to_local_utm(bbox: BBox, **layers: GeoDataFrame) -> tuple[CRS, dict[str, GeoDataFrame], dict[str, int]]

Compute one UTM CRS from bbox, reproject every layer into it, and guarantee finiteness.

Computed once from the bbox itself, not from any individual layer — estimate_utm_crs() can differ between layers covering nearly the same area (e.g. near a UTM zone boundary), which would silently break cross-layer topology if each layer picked its own zone. Takes layers by keyword and returns them by name so that adding a layer cannot quietly reorder an unpacked tuple at a call site.

The third return value counts, per layer, the features that had to be clipped to the study extent to survive the projection at all — see _repair_unprojectable. Zero for most cities.

Source code in src/lczkit/cleaning/pipeline.py

def reproject_to_local_utm(
    bbox: BBox, **layers: gpd.GeoDataFrame
) -> tuple[CRS, dict[str, gpd.GeoDataFrame], dict[str, int]]:
    """Compute one UTM CRS from `bbox`, reproject every layer into it, and guarantee finiteness.

    Computed once from the bbox itself, not from any individual layer — `estimate_utm_crs()`
    can differ between layers covering nearly the same area (e.g. near a UTM zone boundary),
    which would silently break cross-layer topology if each layer picked its own zone. Takes
    layers by keyword and returns them by name so that adding a layer cannot quietly reorder
    an unpacked tuple at a call site.

    The third return value counts, per layer, the features that had to be clipped to the study
    extent to survive the projection at all — see `_repair_unprojectable`. Zero for most cities.
    """
    target = local_utm_crs(bbox)
    projected: dict[str, gpd.GeoDataFrame] = {}
    repaired: dict[str, int] = {}
    for name, layer in layers.items():
        frame, count = _repair_unprojectable(layer.to_crs(target), layer, bbox, target)
        projected[name] = frame
        repaired[name] = count
    return target, projected, repaired

clean_vectors

clean_vectors(source: VectorSource, bbox: BBox, config: CleaningConfig, *, cache_dir: Path | None = None) -> CleanedVectors

Fetch raw vector layers from source and run the full cleaning pipeline.

source is typed against the VectorSource protocol, not any concrete implementation, so this stays source-agnostic — usable with OvertureSource or any future VectorSource.

cache_dir memoises per-tile street simplification when tiling is configured; the caller resolves the path, since config owns every path in this package. None disables it.

Source code in src/lczkit/cleaning/pipeline.py

def clean_vectors(
    source: VectorSource,
    bbox: BBox,
    config: CleaningConfig,
    *,
    cache_dir: Path | None = None,
) -> CleanedVectors:
    """Fetch raw vector layers from `source` and run the full cleaning pipeline.

    `source` is typed against the `VectorSource` protocol, not any concrete implementation, so
    this stays source-agnostic — usable with `OvertureSource` or any future `VectorSource`.

    `cache_dir` memoises per-tile street simplification when tiling is configured; the caller
    resolves the path, since config owns every path in this package. `None` disables it.
    """
    max_area_m2 = _require(config.building_max_area_m2, "building_max_area_m2")
    min_area_m2 = _require(config.building_min_area_m2, "building_min_area_m2")
    merge_limit_m2 = _require(config.building_merge_limit_m2, "building_merge_limit_m2")
    overlap_limit = _require(config.building_overlap_limit, "building_overlap_limit")
    road_buffer_m = _require(config.building_road_buffer_m, "building_road_buffer_m")
    road_overlap_limit = _require(config.building_road_overlap_limit, "building_road_overlap_limit")

    raw_waterlines, raw_waterbodies = source.water(bbox)
    raw = {
        "buildings": source.buildings(bbox),
        "streets": source.streets(bbox),
        "waterlines": raw_waterlines,
        "waterbodies": raw_waterbodies,
        "land_use": source.land_use(bbox),
    }
    crs, layers, repaired = reproject_to_local_utm(bbox, **raw)
    buildings = layers["buildings"]
    streets = layers["streets"]
    waterlines = layers["waterlines"]
    waterbodies = layers["waterbodies"]

    steps: list[CleaningStep] = []
    if any(repaired.values()):
        # Recorded even though it is almost always zero: a feature clipped to the study extent
        # has been changed, and a change nobody can see in the report is the failure mode this
        # phase's own history keeps demonstrating.
        steps.append(
            CleaningStep(
                stage="ingestion",
                operation="clip_unprojectable_features",
                n_in=sum(len(frame) for frame in raw.values()),
                n_out=sum(len(frame) for frame in layers.values()),
                detail={name: count for name, count in repaired.items() if count},
            )
        )

    layers_out, building_steps = clean_buildings(
        buildings,
        max_area_m2=max_area_m2,
        min_area_m2=min_area_m2,
        merge_limit_m2=merge_limit_m2,
        overlap_limit=overlap_limit,
    )
    steps.extend(building_steps)

    # Land use takes no part in cross-layer topology — it is functional metadata, not a
    # physical surface that can conflict with a building or a street.
    land_use, land_use_step = clean_land_use(layers["land_use"])
    steps.append(land_use_step)

    # Simplification comes first: the road-buffer rule buffers the network, and buffering
    # unsimplified dual carriageways would cover the block between them.
    streets, street_step = _simplify(streets, layers_out.topo, config, cache_dir=cache_dir)
    steps.append(street_step)

    buildings_topo, streets, waterlines, waterbodies, topology_steps = apply_cross_layer_topology(
        layers_out.topo,
        streets,
        waterlines,
        waterbodies,
        road_buffer_m=road_buffer_m,
        road_overlap_limit=road_overlap_limit,
    )
    steps.extend(topology_steps)

    return CleanedVectors(
        buildings_area=layers_out.area,
        buildings_topo=buildings_topo,
        streets=streets,
        waterlines=waterlines,
        waterbodies=waterbodies,
        land_use=land_use,
        report=CleaningReport(steps=steps, footprints=layers_out.coverage),
        crs=crs,
    )

Buildings

lczkit.cleaning.buildings

Building-footprint cleaning, producing two layers rather than one.

Every function is a pure transform — (buildings, ...) -> (cleaned, CleaningStep) — with no shared mutable state, so each is testable in isolation. clean_buildings() composes them and is the only function callers outside this module need.

Why two layers. The original single-layer pipeline followed Majer & Fleischmann (arXiv:2603.00132) Supplementary D, where cleaning exists to produce a valid planar partition for tessellation and lost footprint area costs nothing. lczkit feeds the same layer to building surface fraction, which carries roughly 47% of the classification metric, so cleaning for topology silently destroyed the numerator — measured at 23.5% of Berlin's footprint area, worth 9.1 points of agreement. The answer is not weaker cleaning; it is two products with different contracts:

• buildings_area — shared prefix plus overlap trimming only. Feeds building surface fraction, Hr, building count, mean building area and industrial_fraction. Feature-preserving, so building_id is unique here and every area statistic has a complete population.
• buildings_topo — planar and non-overlapping, whatever that costs. Feeds the neatnet exclusion mask and momepy.street_profile. Destructive operations permitted.

Both derive from one shared base and carry building_id from it, so statistics stay joinable.

FEATURE_ID module-attribute

FEATURE_ID = 'feature_id'

Identifier of the source feature a footprint came from, stamped before the multipolygon split.

Several rows share it wherever one Overture feature arrived as a MultiPolygon. That is what makes "one building, one vote" expressible downstream: lczkit.ucp.buildings groups on it so a multi-wing complex contributes one term to Hr and counts once, rather than once per wing.

BUILDING_ID module-attribute

BUILDING_ID = 'building_id'

Stable per-footprint identifier, assigned once on the shared base and carried by both layers.

On buildings_area it is unique. On buildings_topo a dissolved feature keeps one constituent's id — arbitrary, and deliberately not relied on: heights reach buildings_topo by largest-overlap inheritance (lczkit.heights.inherit), not by this join.

RESIDUAL_OVERLAP_EPS_M module-attribute

RESIDUAL_OVERLAP_EPS_M = 1e-06

Width of the buffer subtracted to separate a residual zero-area overlap, in metres.

One micrometre — six orders of magnitude below any survey precision, and eleven below a footprint. It is a numerical tolerance for a floating-point artefact, not a decision about buildings, which is why it lives here rather than in CleaningConfig: nothing about a city would make a different value right. Measured on the Berlin fixture it removes 4.5e-5 m² of 3.13 km², or 1.4e-9 %.

MAX_PLANARITY_PASSES module-attribute

MAX_PLANARITY_PASSES = 8

Bound on the fixed-point loop in enforce_planarity.

Separating one pair can expose another, so the loop is genuinely iterative — Berlin needs two passes for three pairs. The bound exists so a pathological layer terminates instead of hanging.

MAX_PLANARITY_EPS_M module-attribute

MAX_PLANARITY_EPS_M = 0.001

Growth factor per pass, and the ceiling it grows to.

A fixed epsilon cannot converge: a pair the buffer fails to separate at one width will fail at that width however many passes it is given, so the loop spins to its bound and raises. That is what a 9 km² fixture cannot show and 891 km² of Berlin did — three pairs on the fixture all clear at one micrometre, and one pair in the metropolitan extent cleared at none of the eight passes.

The ceiling is one millimetre: still three orders of magnitude below survey precision, and six below a footprint, so escalating this far cannot move a building.

BuildingLayers dataclass

BuildingLayers(area: GeoDataFrame, topo: GeoDataFrame, coverage: FootprintCoverage)

The two layers clean_buildings() forks into, before cross-layer topology runs on topo.

area instance-attribute

area: GeoDataFrame

Area-preserving. Every area statistic reads this.

topo instance-attribute

topo: GeoDataFrame

Planar and non-overlapping. Topology and the street profile read this.

coverage instance-attribute

coverage: FootprintCoverage

Union-based footprint accounting for the area layer, against the raw input.

Carried here rather than returned as a third tuple element so adding it does not change the signature every caller unpacks.

union_area

union_area(buildings: GeoDataFrame) -> float

Ground actually covered by buildings, counting overlapping footprints once.

The denominator the cleaning retention criterion is stated against, and the quantity building surface fraction is trying to measure — BSF sums overlay pieces, so a self-overlapping layer inflates it. See FootprintCoverage for why the sum cannot serve.

Component-wise, because the obvious implementation is superlinear. A single shapely.union_all over every footprint was measured at five Berlin extents before being rejected: exponent 1.26 -> 1.80 in feature count, reaching 711 s for 892k footprints at 891 km2 — roughly doubling a metropolitan clean_vectors run, which is 9.8 minutes in total. Measuring the scaling exponent at several extents is what caught it.

Footprints are very nearly disjoint — Berlin's raw set overlaps itself by 0.19% of summed area at metropolitan extent — so almost every footprint is its own component and the global union is doing enormous work to discover that. Unioning only within groups that genuinely overlap is exact, not an approximation: area is additive over disjoint sets, and two footprints in different components share no area by construction.

Sharing only a boundary is not sharing area, so pairs are filtered on positive intersection area rather than on intersects alone. Otherwise a terrace row becomes one component and the saving goes with it.

Source code in src/lczkit/cleaning/buildings.py

def union_area(buildings: gpd.GeoDataFrame) -> float:
    """Ground actually covered by `buildings`, counting overlapping footprints once.

    The denominator the cleaning retention criterion is stated against, and the quantity building
    surface fraction is trying to measure — BSF sums overlay pieces, so a self-overlapping layer
    inflates it. See `FootprintCoverage` for why the sum cannot serve.

    **Component-wise, because the obvious implementation is superlinear.** A single
    `shapely.union_all` over every footprint was measured at five Berlin extents before being
    rejected: exponent 1.26 -> 1.80 in feature count, reaching **711 s for 892k footprints at
    891 km2** — roughly doubling a metropolitan `clean_vectors` run, which is 9.8 minutes in total.
    Measuring the scaling exponent at several extents is what caught it.

    Footprints are very nearly disjoint — Berlin's raw set overlaps itself by 0.19% of summed area
    at metropolitan extent — so almost every footprint is its own component and the global union is
    doing enormous work to discover that. Unioning only within groups that genuinely overlap is
    exact, not an approximation: area is additive over disjoint sets, and two footprints in
    different components share no area by construction.

    Sharing only a boundary is not sharing area, so pairs are filtered on positive intersection
    area rather than on `intersects` alone. Otherwise a terrace row becomes one component and the
    saving goes with it.
    """
    if buildings.empty:
        return 0.0
    assert_projected_crs(buildings, "buildings")
    geometry = np.asarray(buildings.geometry.values)
    areas = shapely.area(geometry)

    candidates = np.asarray(buildings.sindex.query(buildings.geometry, predicate="intersects"))
    left, right = candidates[0], candidates[1]
    keep = left < right
    left, right = left[keep], right[keep]
    if left.size:
        shared = shapely.area(shapely.intersection(geometry[left], geometry[right]))
        overlapping = shared > 0.0
        left, right = left[overlapping], right[overlapping]

    if not left.size:
        return float(areas.sum())

    n = len(geometry)
    graph = csr_matrix(
        (np.ones(left.size, dtype=np.int8), (left, right)),
        shape=(n, n),
    )
    _, labels = connected_components(graph, directed=False)

    order = np.argsort(labels, kind="stable")
    grouped = labels[order]
    boundaries = np.flatnonzero(np.r_[True, grouped[1:] != grouped[:-1], True])

    total = 0.0
    for start, stop in zip(boundaries[:-1], boundaries[1:], strict=True):
        members = order[start:stop]
        if members.size == 1:
            total += float(areas[members[0]])
        else:
            total += float(shapely.union_all(geometry[members]).area)
    return total

fix_invalid_geometries

fix_invalid_geometries(buildings: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Repair invalid geometries in place via make_valid(). Never changes feature count.

Source code in src/lczkit/cleaning/buildings.py

def fix_invalid_geometries(buildings: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Repair invalid geometries in place via `make_valid()`. Never changes feature count."""
    assert_projected_crs(buildings, "buildings")
    n_invalid = int((~buildings.geometry.is_valid).sum())
    fixed = buildings.copy()
    fixed["geometry"] = fixed.geometry.make_valid()
    return fixed, _step("fix_invalid_geometries", buildings, fixed, n_invalid_before=n_invalid)

explode_multipolygons

explode_multipolygons(buildings: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Split multi-part geometries into single-part rows.

MultiPolygons, and any GeometryCollections left over from make_valid().

Source code in src/lczkit/cleaning/buildings.py

def explode_multipolygons(buildings: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Split multi-part geometries into single-part rows.

    MultiPolygons, and any GeometryCollections left over from `make_valid()`.
    """
    assert_projected_crs(buildings, "buildings")
    exploded = buildings.explode(index_parts=False).reset_index(drop=True)
    return exploded, _step("explode_multipolygons", buildings, exploded)

drop_non_polygons

drop_non_polygons(buildings: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Drop features that are not (non-empty) Polygons.

E.g. stray Points or LineStrings left over from geometry repair, and any empty geometries.

Source code in src/lczkit/cleaning/buildings.py

def drop_non_polygons(buildings: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Drop features that are not (non-empty) Polygons.

    E.g. stray Points or LineStrings left over from geometry repair, and any empty geometries.
    """
    assert_projected_crs(buildings, "buildings")
    keep = (buildings.geometry.geom_type == "Polygon") & (~buildings.geometry.is_empty)
    filtered = buildings.loc[keep].reset_index(drop=True)
    return filtered, _step("drop_non_polygons", buildings, filtered)

drop_oversized

drop_oversized(buildings: GeoDataFrame, max_area_m2: float) -> tuple[GeoDataFrame, CleaningStep]

Drop footprints larger than max_area_m2 — implausible for a single building.

Source code in src/lczkit/cleaning/buildings.py

def drop_oversized(
    buildings: gpd.GeoDataFrame, max_area_m2: float
) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Drop footprints larger than `max_area_m2` — implausible for a single building."""
    assert_projected_crs(buildings, "buildings")
    filtered = buildings.loc[buildings.geometry.area <= max_area_m2].reset_index(drop=True)
    return filtered, _step("drop_oversized", buildings, filtered, max_area_m2=max_area_m2)

assign_feature_id

assign_feature_id(buildings: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Stamp FEATURE_ID onto the raw input, before the shared prefix explodes multipolygons.

This is the only point at which one source feature is still one row. After explode_multipolygons a courtyard block or a multi-wing complex is N rows, and every per-building statistic that treats a row as a building then counts it N times: building_count rises, mean_building_area_m2 falls, and - the one that reaches classification - Hr receives N equal terms in an unweighted geometric mean instead of one.

Not taken from Overture's GERS id, though that would usually work: it is nullable, this package must not assume a particular VectorSource supplies it, and a positional stamp is stable for a given input either way.

Source code in src/lczkit/cleaning/buildings.py

def assign_feature_id(buildings: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Stamp `FEATURE_ID` onto the raw input, *before* the shared prefix explodes multipolygons.

    This is the only point at which one source feature is still one row. After
    `explode_multipolygons` a courtyard block or a multi-wing complex is N rows, and every
    per-building statistic that treats a row as a building then counts it N times: `building_count`
    rises, `mean_building_area_m2` falls, and - the one that reaches classification - `Hr` receives
    N equal terms in an unweighted geometric mean instead of one.

    Not taken from Overture's GERS `id`, though that would usually work: it is nullable, this
    package must not assume a particular `VectorSource` supplies it, and a positional stamp is
    stable for a given input either way.
    """
    assert_projected_crs(buildings, "buildings")
    stamped = buildings.copy()
    stamped[FEATURE_ID] = [f"feat_{i}" for i in range(len(stamped))]
    return stamped, _step("assign_feature_id", buildings, stamped)

assign_building_id

assign_building_id(buildings: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Stamp BUILDING_ID onto the shared base, before the two layers diverge.

Assigned here rather than at ingestion because the shared prefix explodes multipolygons: an id taken from the source would be shared by every part of a multi-part footprint and would not identify a row. Positional, and stable for a given input and configuration.

Source code in src/lczkit/cleaning/buildings.py

def assign_building_id(buildings: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Stamp `BUILDING_ID` onto the shared base, before the two layers diverge.

    Assigned here rather than at ingestion because the shared prefix explodes multipolygons: an id
    taken from the source would be shared by every part of a multi-part footprint and would not
    identify a row. Positional, and stable for a given input and configuration.
    """
    assert_projected_crs(buildings, "buildings")
    stamped = buildings.copy()
    stamped[BUILDING_ID] = [f"bld_{i}" for i in range(len(stamped))]
    return stamped, _step("assign_building_id", buildings, stamped)

trim_overlaps

trim_overlaps(buildings: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Remove the shared part of every pair of overlapping footprints, keeping both features.

This is the only overlap operation buildings_area gets, and it is there for correctness rather than topology: lczkit.ucp.buildings sums overlay pieces per unit, so two footprints overlapping by 50 m² contribute that area twice and building surface fraction can exceed 1.0. Trimming removes exactly the double count. Merging, which would also dissolve the pair into one feature and corrupt building_count and mean_building_area_m2, is topology work and stays on buildings_topo.

Source code in src/lczkit/cleaning/buildings.py

def trim_overlaps(buildings: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Remove the shared part of every pair of overlapping footprints, keeping both features.

    This is the *only* overlap operation `buildings_area` gets, and it is there for correctness
    rather than topology: `lczkit.ucp.buildings` sums overlay pieces per unit, so two footprints
    overlapping by 50 m² contribute that area twice and building surface fraction can exceed 1.0.
    Trimming removes exactly the double count. Merging, which would also dissolve the pair into one
    feature and corrupt `building_count` and `mean_building_area_m2`, is topology work and stays on
    `buildings_topo`.
    """
    assert_projected_crs(buildings, "buildings")
    trimmed = geoplanar.trim_overlaps(buildings, strategy="largest")
    return trimmed, _step("trim_overlaps", buildings, trimmed, stage="buildings_area")

resolve_overlaps

resolve_overlaps(buildings: GeoDataFrame, merge_limit: float, overlap_limit: float) -> tuple[GeoDataFrame, CleaningStep]

Merge overlapping footprints, then trim whatever overlap remains. buildings_topo only.

Merges those below merge_limit, or above it if the shared overlap exceeds overlap_limit. merge_limit and overlap_limit map directly onto geoplanar.merge_overlaps' identically-named parameters.

Source code in src/lczkit/cleaning/buildings.py

def resolve_overlaps(
    buildings: gpd.GeoDataFrame, merge_limit: float, overlap_limit: float
) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Merge overlapping footprints, then trim whatever overlap remains. `buildings_topo` only.

    Merges those below `merge_limit`, or above it if the shared overlap exceeds `overlap_limit`.
    `merge_limit` and `overlap_limit` map directly onto `geoplanar.merge_overlaps`'
    identically-named parameters.
    """
    assert_projected_crs(buildings, "buildings")
    merged = geoplanar.merge_overlaps(buildings, merge_limit, overlap_limit)
    trimmed = geoplanar.trim_overlaps(merged, strategy="largest")
    return trimmed, _step(
        "resolve_overlaps",
        buildings,
        trimmed,
        stage="buildings_topo",
        merge_limit_m2=merge_limit,
        overlap_limit=overlap_limit,
    )

absorb_small_buildings

absorb_small_buildings(buildings: GeoDataFrame, min_area_m2: float) -> tuple[GeoDataFrame, CleaningStep]

Dissolve footprints smaller than min_area_m2 into a touching larger neighbour.

Keeping those that touch nothing. buildings_topo only.

geoplanar.merge_touching deletes any polygon in index that shares no boundary segment with a neighbour, and offers no way to turn that off. Deletion is wrong here: a free-standing garage is small, not spurious, and this operation must dissolve rather than delete. So the small set is partitioned on the same predicate merge_touching uses internally, only the touching part is passed to it, and the isolates are concatenated back untouched.

Measured on the Berlin fixture: 1186 footprints under 20 m², of which 1043 are isolated. The deletion was worth 0.12% of footprint area — a real bug, and not the one that mattered.

Source code in src/lczkit/cleaning/buildings.py

def absorb_small_buildings(
    buildings: gpd.GeoDataFrame, min_area_m2: float
) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Dissolve footprints smaller than `min_area_m2` into a touching larger neighbour.

    **Keeping those that touch nothing.** `buildings_topo` only.

    `geoplanar.merge_touching` deletes any polygon in `index` that shares no boundary segment with
    a neighbour, and offers no way to turn that off. Deletion is wrong here: a free-standing garage
    is small, not spurious, and this operation must dissolve rather than delete. So the small set
    is partitioned on the same predicate `merge_touching` uses internally,
    only the touching part is passed to it, and the isolates are concatenated back untouched.

    Measured on the Berlin fixture: 1186 footprints under 20 m², of which 1043 are isolated. The
    deletion was worth 0.12% of footprint area — a real bug, and not the one that mattered.
    """
    assert_projected_crs(buildings, "buildings")
    small = buildings.index[buildings.geometry.area < min_area_m2]
    if small.empty:
        return buildings, _step(
            "absorb_small_buildings",
            buildings,
            buildings,
            stage="buildings_topo",
            min_area_m2=min_area_m2,
            n_small=0,
            n_dissolved=0,
            n_isolated_retained=0,
        )

    # `source` indexes positionally into the query geometries, i.e. into `small`. This is the same
    # predicate `merge_touching` applies internally, so the partition matches exactly what it would
    # have dissolved and what it would have deleted.
    source, _ = buildings.boundary.sindex.query(buildings.loc[small].boundary, predicate="overlaps")
    touching = small[sorted(set(source.tolist()))]
    isolated = small.difference(touching)

    if touching.empty:
        merged = buildings
    else:
        merged = geoplanar.merge_touching(
            buildings.drop(index=isolated), index=touching.tolist(), largest=True
        )
        if not isolated.empty:
            merged = gpd.GeoDataFrame(
                pd.concat([merged, buildings.loc[isolated]], ignore_index=True),
                geometry="geometry",
                crs=buildings.crs,
            )

    return merged, _step(
        "absorb_small_buildings",
        buildings,
        merged,
        stage="buildings_topo",
        min_area_m2=min_area_m2,
        n_small=len(small),
        n_dissolved=len(touching),
        n_isolated_retained=len(isolated),
    )

enforce_planarity

enforce_planarity(buildings: GeoDataFrame, *, eps_m: float = RESIDUAL_OVERLAP_EPS_M, max_passes: int = MAX_PLANARITY_PASSES) -> tuple[GeoDataFrame, CleaningStep]

Clear the overlaps geoplanar.trim_overlaps cannot. buildings_topo only.

Why anything is left to clear. trim_overlaps subtracts one footprint of an overlapping pair from the other with a plain difference. Where the pair's overlap has collapsed to zero area — two polygons sharing a boundary that is collinear but carries mismatched vertices — that subtraction is a no-op, and the pair stays flagged however many times it is run. Measured on the Berlin fixture: three such pairs survive resolve_overlaps, relating as 212111212 (a 2-D interior intersection) while intersection() returns a MultiLineString of area exactly 0.0. That single flag is what made buildings_topo report is_planar_enforced: False for two phases, while momepy.enclosures() — which needs a planar input — was reading it.

Neither shapely.set_precision nor shapely.difference(grid_size=...) fixes them reliably; tried across 1e-9 to 1e-2, each resolved at most one of the three and coarser grids introduced new overlaps elsewhere. Subtracting a eps_m-wide buffer of the smaller polygon from the larger does resolve all three, because it forces a genuine re-noding of the shared boundary.

The larger of each pair is trimmed, matching trim_overlaps(strategy="largest"), so which feature loses the sliver does not depend on which operation reached it first.

Source code in src/lczkit/cleaning/buildings.py

def enforce_planarity(
    buildings: gpd.GeoDataFrame,
    *,
    eps_m: float = RESIDUAL_OVERLAP_EPS_M,
    max_passes: int = MAX_PLANARITY_PASSES,
) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Clear the overlaps `geoplanar.trim_overlaps` cannot. `buildings_topo` only.

    **Why anything is left to clear.** `trim_overlaps` subtracts one footprint of an overlapping
    pair from the other with a plain `difference`. Where the pair's overlap has collapsed to zero
    area — two polygons sharing a boundary that is collinear but carries mismatched vertices — that
    subtraction is a no-op, and the pair stays flagged however many times it is run. Measured on the
    Berlin fixture: three such pairs survive `resolve_overlaps`, relating as `212111212` (a 2-D
    interior intersection) while `intersection()` returns a MultiLineString of area exactly 0.0.
    That single flag is what made `buildings_topo` report `is_planar_enforced: False` for two
    phases, while `momepy.enclosures()` — which needs a planar input — was reading it.

    Neither `shapely.set_precision` nor `shapely.difference(grid_size=...)` fixes them reliably;
    tried across 1e-9 to 1e-2, each resolved at most one of the three and coarser grids introduced
    new overlaps elsewhere. Subtracting a `eps_m`-wide buffer of the smaller polygon from the larger
    does resolve all three, because it forces a genuine re-noding of the shared boundary.

    The larger of each pair is trimmed, matching `trim_overlaps(strategy="largest")`, so which
    feature loses the sliver does not depend on which operation reached it first.
    """
    assert_projected_crs(buildings, "buildings")
    working = buildings.reset_index(drop=True)
    # Edited as a positional array rather than through `.loc`: the spatial index reports positions,
    # and a pass rewrites the same footprint more than once when it overlaps two neighbours.
    geometries = working.geometry.to_numpy().copy()

    passes = 0
    n_pairs = 0
    eps = eps_m
    eps_used = eps_m
    pairs = _overlapping_pairs(working)
    while pairs and passes < max_passes:
        passes += 1
        n_pairs = max(n_pairs, len(pairs))
        # Record the width this pass actually subtracts, before the escalation at the foot of the
        # loop. Deriving it afterwards by dividing back is wrong once `eps` saturates at the
        # ceiling, because the last two passes then share a width and the division reports one
        # that was never used.
        eps_used = eps

        touched: set[int] = set()
        for first, second in pairs:
            larger = first if geometries[first].area >= geometries[second].area else second
            smaller = second if larger == first else first
            geometries[larger] = geometries[larger].difference(geometries[smaller].buffer(eps))
            touched.add(larger)

        positions = sorted(touched)
        geometries[positions] = largest_part(
            gpd.GeoSeries(geometries[positions], crs=working.crs)
        ).to_numpy()
        working = working.set_geometry(
            gpd.GeoSeries(geometries, index=working.index, crs=working.crs)
        )
        pairs = _overlapping_pairs(working)
        # Widen for the next pass: the same width that just failed on a pair will keep failing.
        eps = min(eps * PLANARITY_EPS_ESCALATION, MAX_PLANARITY_EPS_M)

    # Last resort: drop the smaller footprint of whatever still overlaps. `buildings_topo` is the
    # destructive layer and the one that must be planar for `momepy.enclosures()`; `buildings_area`
    # carries every area statistic and is untouched by this. Recorded, never silent — an
    # unexplained gap in the topology layer is exactly what the cleaning report exists to make
    # visible, and it is preferable to ending a run over a whole city on one bad pair.
    unresolvable = sorted(
        {
            second if geometries[first].area >= geometries[second].area else first
            for first, second in pairs
        }
    )
    if unresolvable:
        working = working.drop(index=working.index[unresolvable])

    working = working.loc[working.geometry.notna() & ~working.geometry.is_empty].reset_index(
        drop=True
    )
    return working, _step(
        "enforce_planarity",
        buildings,
        working,
        stage="buildings_topo",
        eps_m=eps_m,
        eps_final_m=eps_used,
        n_passes=passes,
        n_residual_pairs=n_pairs,
        n_dropped_unresolvable=len(unresolvable),
        area_removed_m2=_area(buildings) - _area(working),
    )

clean_buildings

clean_buildings(buildings: GeoDataFrame, *, max_area_m2: float, min_area_m2: float, merge_limit_m2: float, overlap_limit: float) -> tuple[BuildingLayers, list[CleaningStep]]

Run the shared prefix, then fork into the area-preserving and topological layers.

buildings_topo returns without cross-layer topology applied — the road-buffer rule needs the simplified street network, which pipeline.clean_vectors() produces after this point.

Source code in src/lczkit/cleaning/buildings.py

def clean_buildings(
    buildings: gpd.GeoDataFrame,
    *,
    max_area_m2: float,
    min_area_m2: float,
    merge_limit_m2: float,
    overlap_limit: float,
) -> tuple[BuildingLayers, list[CleaningStep]]:
    """Run the shared prefix, then fork into the area-preserving and topological layers.

    `buildings_topo` returns without cross-layer topology applied — the road-buffer rule needs the
    *simplified* street network, which `pipeline.clean_vectors()` produces after this point.
    """
    steps: list[CleaningStep] = []
    base, step = assign_feature_id(buildings)
    steps.append(step)
    for operation in (
        fix_invalid_geometries,
        explode_multipolygons,
        drop_non_polygons,
    ):
        base, step = operation(base)
        steps.append(step)
    base, step = drop_oversized(base, max_area_m2)
    steps.append(step)
    base, step = assign_building_id(base)
    steps.append(step)

    area, step = trim_overlaps(base)
    steps.append(step)

    # Measured on `base`, i.e. after the validity fixes and before the fork, which is how the
    # criterion is stated: at least 99% of the union of input footprint area is retained after
    # validity fixes. Taking it on the raw input instead would fold make_valid's repairs into
    # the denominator and make the criterion depend on how broken the source geometry was.
    coverage = FootprintCoverage(
        raw_summed_area_m2=_area(base),
        raw_union_area_m2=union_area(base),
        area_summed_m2=_area(area),
        area_union_m2=union_area(area),
    )

    topo, step = resolve_overlaps(base, merge_limit_m2, overlap_limit)
    steps.append(step)
    topo, step = absorb_small_buildings(topo, min_area_m2)
    steps.append(step)

    # Last, not earlier: `absorb_small_buildings` dissolves, and a dissolve can reintroduce the
    # artefact this clears. The invariant has to be established after the last thing that can
    # break it, or `validate_planarity` below is measuring a layer that no longer exists.
    topo, step = enforce_planarity(topo)
    steps.append(step)

    # allow_gaps=True: buildings are not a plane tessellation — gaps between separate
    # buildings are normal and expected; only overlaps are a real violation here.
    planar = bool(geoplanar.is_planar_enforced(topo, allow_gaps=True))
    steps.append(
        CleaningStep(
            stage="buildings_topo",
            operation="validate_planarity",
            n_in=len(topo),
            n_out=len(topo),
            area_in_m2=_area(topo),
            area_out_m2=_area(topo),
            detail={"is_planar_enforced": planar},
        )
    )
    return BuildingLayers(area=area, topo=topo, coverage=coverage), steps

Streets

lczkit.cleaning.streets

Street-network simplification via neatnet, whole-extent and tiled.

neatnet owns its own tuning parameters (node-merge tolerance, continuity-stroke angle threshold, and a dozen others) as documented package defaults — lczkit does not proxy them into config; doing so would wrap an entire downstream library's parameter surface for no current need.

Two entry points, same result shape:

• simplify_streets runs neatify over the whole extent. Correct at any size, but superlinear — see lczkit.cleaning.tiles.
• simplify_streets_tiled splits the extent, simplifies each tile independently over a buffered window, and stitches the cores back together.

The tiled path pins one face-artifact threshold across every tile. neatnet derives that threshold from the distribution of face-artifact-index values across whatever network it is handed — a kernel-density valley — so a tile left to itself computes a different threshold from its neighbour. On a 2x2 tiling of 16 km2 of Berlin the whole extent found no valley and fell back to 7.0, while two of the four tiles found 8.10 and 7.58: a face with an index of 7.5 would have been an artifact in one tile and ordinary urban fabric in the run next to it. Pinning one value across all tiles removes that class of seam disagreement outright.

Simplification is sensitive to input row order. neatnet re-nodes and re-merges a network in the order it receives it: a shuffled input yields the same feature count and the same total length with the edges split at different points, which reaches momepy.street_profile and so aspect_ratio. Nothing here can fix that, and nothing here should paper over it — what it means is that the row order arriving from a VectorSource has to be canonical already, or two runs over the same city disagree. lczkit.sources.overture._canonical_order is where that is established; this module only depends on it.

Where that value comes from is the second scaling problem this module had. Deriving it from the whole network, as resolve_artifact_threshold does, means running neatnet.fix_topology over the whole network — quadratic in feature count, measured at exponent 2.0 and projecting to ~8.6 hours at metropolitan scale, against 7.5 minutes for the tiles themselves. Detection is not the expensive part; the preprocessing it needs is, at 8.3 s against 12 392 s at 484 km2. The index is a per-face quantity, so pooled_artifact_threshold assembles the same distribution from the per-tile windows at k * (n/k)**2, in parallel, and resolve_artifact_threshold remains as the reference the pooled one is measured against.

ARTIFACT_THRESHOLD_FALLBACK module-attribute

ARTIFACT_THRESHOLD_FALLBACK = 7.0

neatnet.neatify's own default for when no face-artifact-index valley is found.

Restated here because the tiled path resolves the threshold itself and must fall back exactly as neatify would, so that a one-tile run and a whole-extent run agree.

TILE_RESULT_VERSION module-attribute

TILE_RESULT_VERSION = 3

Bumped whenever _simplify_window changes what it writes for a given input.

Part of the per-tile cache key. Config and the neatnet version cover everything outside this module that changes a tile's contents; this covers what is inside it, which nothing else in the key would detect.

Version 2 adds _tile_key and _failure, so a cached tile replays the report as well as the geometry.

Version 3: subset now preserves the layer's row order rather than returning spatial-index order, which changes the linework neatnet produces for an unchanged input. The rest of the key does not notice. Measured at 64 and 144 km2 of Berlin, the pooled threshold is identical under both orderings, so _threshold_tag moves not at all while tile contents differ by ~1.2% of linework — exactly the case this field exists for, and one that would otherwise have served pre-fix tiles to a post-fix run with nothing reporting it.

SIMPLIFIED_COLUMN module-attribute

SIMPLIFIED_COLUMN = '_simplified'

Per-edge flag: False marks linework that passed through a tile neatnet could not process.

Carried out of cleaning rather than reduced to a count, so a downstream oddity can be traced to the tile that produced it instead of being attributed to the classifier.

TILE_KEY_COLUMN module-attribute

TILE_KEY_COLUMN = '_tile_key'

Per-edge tile provenance: which tile emitted this line.

SIMPLIFIED_COLUMN alone says an edge came from a tile that failed, not which tile, so it could not actually be traced back — the claim that it could was made and never exercised. With the key present, a suspect enclosure downstream resolves to one tile and one window.

FAILURE_COLUMN module-attribute

FAILURE_COLUMN = '_failure'

Internal, dropped before the stitch: the exception class that made a tile pass through.

Written into the per-tile cache so a cached run reconstructs the same cleaning report a cold run would have produced. A cache that reproduces the geometry but not the record of how it was obtained is still not transparent.

THREAD_LIMIT_VARS module-attribute

THREAD_LIMIT_VARS = ('OMP_NUM_THREADS', 'OPENBLAS_NUM_THREADS', 'MKL_NUM_THREADS', 'NUMEXPR_NUM_THREADS')

The native thread-pool controls every worker must see set to 1.

Each worker runs GEOS and BLAS through geopandas, and each would otherwise start a thread pool sized to the whole machine. Thirty-two workers times thirty-two threads is not parallelism, it is the same oversubscription that makes the pool appear hung.

TileFaceIndex

Bases: NamedTuple

One tile's contribution to the pooled face-artifact-index distribution.

PooledThreshold dataclass

PooledThreshold(value: float, n_faces: int, n_faces_dropped: int, dropped_area_fraction: float, n_tiles_indexed: int, fallback_used: bool)

A face-artifact threshold pooled from per-tile distributions, with its error term.

as_detail

as_detail() -> dict[str, object]

The cleaning-report form. Every field travels; the error term is not a debug aid.

Source code in src/lczkit/cleaning/streets.py

def as_detail(self) -> dict[str, object]:
    """The cleaning-report form. Every field travels; the error term is not a debug aid."""
    return {
        "artifact_threshold": self.value,
        "threshold_source": "pooled",
        "threshold_fallback_used": self.fallback_used,
        "threshold_n_faces": self.n_faces,
        "threshold_n_faces_dropped": self.n_faces_dropped,
        "threshold_dropped_area_fraction": self.dropped_area_fraction,
        "threshold_n_tiles_indexed": self.n_tiles_indexed,
    }

simplify_streets

simplify_streets(streets: GeoDataFrame, buildings: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Simplify streets with neatnet.neatify, using buildings as the exclusion mask.

Required, not optional: unsimplified dual carriageways and roundabouts destroy enclosure generation downstream.

Source code in src/lczkit/cleaning/streets.py

def simplify_streets(
    streets: gpd.GeoDataFrame, buildings: gpd.GeoDataFrame
) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Simplify `streets` with `neatnet.neatify`, using `buildings` as the exclusion mask.

    Required, not optional: unsimplified dual carriageways and roundabouts destroy enclosure
    generation downstream.
    """
    assert_projected_crs(streets, "streets")
    assert_projected_crs(buildings, "buildings")
    simplified = neatnet.neatify(streets, exclusion_mask=buildings.geometry)
    step = CleaningStep(
        stage="streets",
        operation="simplify_streets",
        n_in=len(streets),
        n_out=len(simplified),
        detail={"tiled": False},
    )
    return simplified, step

resolve_artifact_threshold

resolve_artifact_threshold(streets: GeoDataFrame, *, fallback: float = ARTIFACT_THRESHOLD_FALLBACK) -> float

The face-artifact-index threshold neatify would derive from streets as a whole.

Computed once per extent and pushed into every tile — see this module's docstring for why that matters.

Measured after the same preprocessing neatify applies, not on the raw input. neatify runs fix_topology and then consolidate_nodes before it indexes any face, and those change which faces exist: an unnoded crossing forms no face at all, so a raw network can look artifact-free when the network actually simplified is not. Skipping this yields a threshold from a different network than the one it is pinned onto, which is the failure this function exists to prevent rather than one it may commit.

Returns fallback when no kernel-density valley is found, which is what neatify does internally with its artifact_threshold_fallback. Degenerate networks are the fallback's other job: one too sparse to polygonize has no faces to index, and a perfectly uniform one gives every face the same index, which makes the kernel density estimate singular and raises out of scipy. Neither is worth failing a run over, and at metropolitan scale neither can be ruled out.

Source code in src/lczkit/cleaning/streets.py

def resolve_artifact_threshold(
    streets: gpd.GeoDataFrame, *, fallback: float = ARTIFACT_THRESHOLD_FALLBACK
) -> float:
    """The face-artifact-index threshold `neatify` would derive from `streets` as a whole.

    Computed once per extent and pushed into every tile — see this module's docstring for why
    that matters.

    **Measured after the same preprocessing `neatify` applies**, not on the raw input.
    `neatify` runs `fix_topology` and then `consolidate_nodes` before it indexes any face, and
    those change which faces exist: an unnoded crossing forms no face at all, so a raw network
    can look artifact-free when the network actually simplified is not. Skipping this yields a
    threshold from a different network than the one it is pinned onto, which is the failure this
    function exists to prevent rather than one it may commit.

    Returns `fallback` when no kernel-density valley is found, which is what `neatify` does
    internally with its `artifact_threshold_fallback`. Degenerate networks are the fallback's
    other job: one too sparse to polygonize has no faces to index, and a perfectly uniform one
    gives every face the *same* index, which makes the kernel density estimate singular and
    raises out of scipy. Neither is worth failing a run over, and at metropolitan scale neither
    can be ruled out.
    """
    assert_projected_crs(streets, "streets")
    try:
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", message="No threshold found")
            warnings.filterwarnings("ignore", message="Input streets could not")
            threshold = neatnet.FaceArtifacts(_preprocess(streets)).threshold
    except (ValueError, np.linalg.LinAlgError):
        return fallback
    return fallback if threshold is None else float(threshold)

pooled_artifact_threshold

pooled_artifact_threshold(streets: GeoDataFrame, tiles: list[Tile], *, workers: int = 1, fallback: float = ARTIFACT_THRESHOLD_FALLBACK) -> PooledThreshold

Pin the artifact threshold from the tiles, instead of from the whole network.

resolve_artifact_threshold gets the same number by running neatnet.fix_topology over the entire extent, which is quadratic in feature count — measured at 481.7 s, 1916.9 s, 5004.6 s and 12392.6 s for 33.8k, 67.2k, 106.7k and 168.5k Berlin streets, an exponent of 2.0, and extrapolating to roughly 8.6 hours at the 267k streets of the metropolitan extent. That is the whole of the gap between the 7.5 minutes the tiles actually take and the 15 hours the first metropolitan run spent before it was killed.

Each tile already runs fix_topology over its own window inside neatify, so the same distribution can be assembled from k windows at k * (n/k)**2 — and, unlike the whole-network step, it parallelises.

This is an approximation with a measured error, not an identity. Faces spanning a window boundary are absent from the pool (see _tile_face_index), so the pooled distribution is missing its largest faces. dropped_area_fraction is that error. Measured against the whole-network threshold over six extents, the two estimators converge — the pooled value settles at 8.1876 against 8.1918, and the deviation shrinks as the extent grows.

Source code in src/lczkit/cleaning/streets.py

def pooled_artifact_threshold(
    streets: gpd.GeoDataFrame,
    tiles: list[Tile],
    *,
    workers: int = 1,
    fallback: float = ARTIFACT_THRESHOLD_FALLBACK,
) -> PooledThreshold:
    """Pin the artifact threshold from the tiles, instead of from the whole network.

    `resolve_artifact_threshold` gets the same number by running `neatnet.fix_topology` over the
    entire extent, which is **quadratic in feature count** — measured at 481.7 s, 1916.9 s,
    5004.6 s and 12392.6 s for 33.8k, 67.2k, 106.7k and 168.5k Berlin streets, an exponent of
    2.0, and extrapolating to roughly 8.6 hours at the 267k streets of the metropolitan extent.
    That is the whole of the gap between the 7.5 minutes the tiles actually take and the 15 hours
    the first metropolitan run spent before it was killed.

    Each tile already runs `fix_topology` over its own window inside `neatify`, so the same
    distribution can be assembled from k windows at k * (n/k)**2 — and, unlike the whole-network
    step, it parallelises.

    **This is an approximation with a measured error, not an identity.** Faces spanning a window
    boundary are absent from the pool (see `_tile_face_index`), so the pooled distribution is
    missing its largest faces. `dropped_area_fraction` is that error. Measured against the
    whole-network threshold over six extents, the two estimators converge — the pooled value
    settles at 8.1876 against 8.1918, and the deviation shrinks as the extent grows.
    """
    assert_projected_crs(streets, "streets")
    jobs = [(tile, subset(streets, tile.window)) for tile in tiles]
    n_workers = max(1, min(workers, len(jobs)))
    # The thread pinning wraps *both* branches. `n_workers` follows `os.sched_getaffinity`, so
    # whether this runs serially is a property of the machine, and the threshold it produces is
    # the tile cache key at full float precision. Pinning only the parallel branch would let the
    # same extent on a differently-sized node land on a different key and silently rebuild every
    # tile - a cache that misses for reasons the report cannot show.
    with _single_threaded_children():
        if n_workers == 1:
            indexed = [_tile_face_index(*job) for job in jobs]
        else:
            with _worker_pool(n_workers) as pool:
                indexed = list(pool.map(_tile_face_index, *zip(*jobs, strict=True)))

    populated = [tile_index for tile_index in indexed if tile_index.values.size]
    values = (
        np.concatenate([tile_index.values for tile_index in populated])
        if populated
        else np.empty(0, dtype=float)
    )
    n_dropped = sum(tile_index.n_dropped for tile_index in indexed)
    dropped_area = sum(tile_index.dropped_area_m2 for tile_index in indexed)
    kept_area = sum(tile_index.kept_area_m2 for tile_index in indexed)
    total_area = dropped_area + kept_area

    threshold = None
    if values.size > 1:
        try:
            threshold = _threshold_from_index(values)
        except (ValueError, np.linalg.LinAlgError):
            threshold = None
    return PooledThreshold(
        value=fallback if threshold is None else threshold,
        n_faces=int(values.size),
        n_faces_dropped=n_dropped,
        dropped_area_fraction=dropped_area / total_area if total_area else 0.0,
        n_tiles_indexed=len(populated),
        fallback_used=threshold is None,
    )

simplify_streets_tiled

simplify_streets_tiled(streets: GeoDataFrame, buildings: GeoDataFrame, *, tile_size_m: float, buffer_m: float, workers: int | None = None, cache_dir: Path | None = None, cache_fingerprint: str = 'default', artifact_threshold: float | None = None) -> tuple[GeoDataFrame, CleaningStep]

Tile the extent, simplify each tile over a buffered window, and stitch the cores.

workers defaults to every core the process is allowed to use. cache_dir, when given, memoises each tile under <cache_dir>/<cache_fingerprint>/<tile_key>.parquet; the fingerprint must cover everything that changes a tile's result — the Overture release, the tile geometry, and the pinned threshold — because tile keys alone are CRS-origin-aligned and therefore shared across runs.

artifact_threshold pins the face-artifact threshold explicitly; None pools it from the tiles via pooled_artifact_threshold.

Returns the same (GeoDataFrame, CleaningStep) shape as simplify_streets, so the caller does not branch on which path ran.

Source code in src/lczkit/cleaning/streets.py

def simplify_streets_tiled(
    streets: gpd.GeoDataFrame,
    buildings: gpd.GeoDataFrame,
    *,
    tile_size_m: float,
    buffer_m: float,
    workers: int | None = None,
    cache_dir: Path | None = None,
    cache_fingerprint: str = "default",
    artifact_threshold: float | None = None,
) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Tile the extent, simplify each tile over a buffered window, and stitch the cores.

    `workers` defaults to every core the process is allowed to use. `cache_dir`, when given,
    memoises each tile under `<cache_dir>/<cache_fingerprint>/<tile_key>.parquet`; the
    fingerprint must cover everything that changes a tile's result — the Overture release, the
    tile geometry, and the pinned threshold — because tile keys alone are CRS-origin-aligned and
    therefore shared across runs.

    `artifact_threshold` pins the face-artifact threshold explicitly; `None` pools it from the
    tiles via `pooled_artifact_threshold`.

    Returns the same `(GeoDataFrame, CleaningStep)` shape as `simplify_streets`, so the caller
    does not branch on which path ran.
    """
    assert_projected_crs(streets, "streets")
    assert_projected_crs(buildings, "buildings")

    extent = layer_extent(streets)
    tiles = build_tiles(extent, tile_size_m=tile_size_m, buffer_m=buffer_m, crs_hint="streets")
    n_workers = workers if workers is not None else len(os.sched_getaffinity(0))
    n_workers = max(1, min(n_workers, len(tiles)))

    if artifact_threshold is None:
        pooled = pooled_artifact_threshold(streets, tiles, workers=n_workers)
        threshold, threshold_detail = pooled.value, pooled.as_detail()
    else:
        threshold = artifact_threshold
        threshold_detail = {"artifact_threshold": threshold, "threshold_source": "configured"}

    # The threshold reaches the cache key because it is pinned from the whole extent: the same
    # tile genuinely simplifies differently under a different study area, and without this a run
    # over a larger city would read back tiles decided by the smaller one's threshold.
    # `TILE_RESULT_VERSION` covers the other half — a change to what `_simplify_window` writes,
    # which no amount of config hashing would notice.
    fingerprint = f"{cache_fingerprint}_v{TILE_RESULT_VERSION}_thr{_threshold_tag(threshold)}"
    jobs = [
        (
            tile,
            subset(streets, tile.window),
            subset(buildings, tile.window),
            threshold,
            None if cache_dir is None else _cache_path(cache_dir, tile, fingerprint),
            shared_edges(tile, tiles),
        )
        for tile in tiles
    ]

    # Counted before anything runs, because `_simplify_window` writes the file it would have
    # read. Asking afterwards reports every tile as a hit and tells you nothing.
    n_cached = sum(1 for job in jobs if job[4] is not None and job[4].exists())

    if n_workers == 1:
        parts = [_simplify_window(*job) for job in jobs]
    else:
        with _single_threaded_children(), _worker_pool(n_workers) as pool:
            parts = list(pool.map(_simplify_window, *zip(*jobs, strict=True)))

    # Aggregated before the stitch and sorted by tile key, so two runs that degraded on
    # different tiles are distinguishable from the report alone rather than by re-running.
    # `pd.notna`, not `is not None`: an all-null object column round-trips out of parquet as a
    # float NaN column, so a cached healthy tile would otherwise be reported as having failed
    # with the reason "nan" — a bug visible only on the cache path, which is the worst kind.
    passed_through = {
        str(part[TILE_KEY_COLUMN].iloc[0]): str(part[FAILURE_COLUMN].iloc[0])
        for part in parts
        if len(part) and pd.notna(part[FAILURE_COLUMN].iloc[0])
    }
    parts = [part.drop(columns=[FAILURE_COLUMN]) for part in parts]

    simplified = _stitch(parts, streets.crs)
    step = CleaningStep(
        stage="streets",
        operation="simplify_streets_tiled",
        n_in=len(streets),
        n_out=len(simplified),
        detail={
            "tiled": True,
            "n_tiles": len(tiles),
            "tile_size_m": tile_size_m,
            "buffer_m": buffer_m,
            "workers": n_workers,
            # Two separate facts. The old single `cached` flag recorded only that a directory
            # was configured, and was read as "this run reused tiles" — which is how a run that
            # had computed all 594 of its own tiles came to be reported as a cached run, and a
            # cold/cold difference came to be written up as a cache defect.
            "cache_dir_configured": cache_dir is not None,
            "n_tiles_reused": n_cached,
            "n_tiles_unsimplified": len(passed_through),
            "tiles_unsimplified": dict(sorted(passed_through.items())),
            **threshold_detail,
        },
    )
    return simplified, step

seam_disagreement

seam_disagreement(tiled: GeoDataFrame, untiled: GeoDataFrame, *, tolerance_m: float = 0.5) -> dict[str, float]

Length of linework present in one simplification and not the other, in kilometres.

The correctness instrument for the tiled path, and what its tests assert on. Both inputs are reduced to a single geometry and compared with a tolerance_m buffer, so a street counts as agreeing when it follows the same line, not when its vertices are equal — neatnet re-nodes and re-merges, so vertex equality is too strong a test.

Compare only over ground both runs cover: a whole-extent run keeps streets that merely intersect its window and so reaches beyond it, while the tiled run is clipped to tile cores. Comparing the two unclipped exaggerates disagreement by the perimeter overhang.

Source code in src/lczkit/cleaning/streets.py

def seam_disagreement(
    tiled: gpd.GeoDataFrame, untiled: gpd.GeoDataFrame, *, tolerance_m: float = 0.5
) -> dict[str, float]:
    """Length of linework present in one simplification and not the other, in kilometres.

    The correctness instrument for the tiled path, and what its tests assert on.
    Both inputs are reduced to a single geometry and compared with a `tolerance_m` buffer, so a
    street counts as agreeing when it follows the same line, not when its vertices are equal —
    neatnet re-nodes and re-merges, so vertex equality is too strong a test.

    Compare only over ground both runs cover: a whole-extent run keeps streets that merely
    *intersect* its window and so reaches beyond it, while the tiled run is clipped to tile
    cores. Comparing the two unclipped exaggerates disagreement by the perimeter overhang.
    """
    left = tiled.geometry.union_all()
    right = untiled.geometry.union_all()
    overlap = shapely.box(*left.bounds).intersection(shapely.box(*right.bounds))
    common = shapely.box(*overlap.bounds)
    left = left.intersection(common)
    right = right.intersection(common)
    missing = right.difference(left.buffer(tolerance_m)).length
    extra = left.difference(right.buffer(tolerance_m)).length
    return {
        "tiled_km": left.length / 1000.0,
        "untiled_km": right.length / 1000.0,
        "missing_km": missing / 1000.0,
        "extra_km": extra / 1000.0,
        "agreement": 1.0 - missing / right.length if right.length else 1.0,
    }

Tiling

Street simplification is superlinear in the size of the network, so a whole city does not finish in usable time. This is how an extent is cut into tiles, simplified in parallel and stitched back.

lczkit.cleaning.tiles

Square tiles with buffered working windows, for chunking work that is superlinear in extent.

This exists because neatnet.neatify is superlinear in the extent handed to it. Measured on Berlin, with everything else in clean_vectors held constant:

extent    streets   neatify    largest rook-connected artifact component
  1 km2       706      9.6 s       538
  9 km2      6428     93.4 s      3658
 36 km2     21168    580.3 s     13348
100 km2     51682   2983.3 s     32966

The exponent in area climbs from 0.95 (1 -> 4 km2) to 1.67 (36 -> 64 km2), because face artifacts percolate: 93 percent of them fall into a single rook-contiguous component, and neatify_clusters simplifies each component as one unit. Extent does not merely add work, it enlarges the single largest piece of work.

Superlinearity is exactly what makes tiling pay. Splitting an extent into k tiles cuts the largest component roughly k-fold, so the total shrinks by about k**(1 - p) for exponent p before any tile runs in parallel.

Tiles are aligned to the projected CRS's own coordinate origin rather than to the extent, the same convention lczkit.units.grid.GridUnits uses. Two overlapping extents therefore produce identical tile boundaries wherever they overlap, which is what lets a per-tile result be cached and reused across runs whose bboxes differ.

Tile dataclass

Tile(col: int, row: int, core: Polygon, window: Polygon)

One tile: an exclusive core and the larger window actually computed over.

core tiles partition the covered area exactly — no gaps, no overlaps — so concatenating per-tile results clipped to their cores neither duplicates nor drops ground. window is core expanded by the buffer, giving each tile enough surrounding context that a feature near the seam is decided with its neighbourhood present rather than truncated.

key property

key: str

Stable identifier: the tile's position on the CRS-origin-aligned grid.

Independent of which extent asked for it, so it is usable as a cache key.

build_tiles

build_tiles(extent: Polygon, *, tile_size_m: float, buffer_m: float, crs_hint: str = 'extent') -> list[Tile]

Tile extent into tile_size_m squares, each with a buffer_m working margin.

extent must be in a projected CRS — the caller enforces that; the sizes here are metres. Only tiles whose core actually meets extent are returned, so an L-shaped or diagonal study area does not pay for the empty corners of its bounding box.

Source code in src/lczkit/cleaning/tiles.py

def build_tiles(
    extent: Polygon, *, tile_size_m: float, buffer_m: float, crs_hint: str = "extent"
) -> list[Tile]:
    """Tile `extent` into `tile_size_m` squares, each with a `buffer_m` working margin.

    `extent` must be in a projected CRS — the caller enforces that; the sizes here are metres.
    Only tiles whose core actually meets `extent` are returned, so an L-shaped or diagonal
    study area does not pay for the empty corners of its bounding box.
    """
    if tile_size_m <= 0:
        raise ValueError(f"tile_size_m must be positive, got {tile_size_m}")
    if buffer_m < 0:
        raise ValueError(f"buffer_m must be non-negative, got {buffer_m}")
    if extent.is_empty:
        raise ValueError(f"{crs_hint} is empty; nothing to tile")

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

    tiles: list[Tile] = []
    for col in range(col_start, col_end + 1):
        x0 = col * tile_size_m
        for row in range(row_start, row_end + 1):
            y0 = row * tile_size_m
            core = box(x0, y0, x0 + tile_size_m, y0 + tile_size_m)
            # `touches` rather than `intersects`: an extent whose maximum falls exactly on a
            # tile boundary — which is every extent snapped to a round grid — otherwise spawns a
            # column of tiles meeting it along a line and holding no area at all. Those tiles
            # then fail artifact detection, since a network with no faces cannot be indexed.
            if not core.intersects(extent) or core.touches(extent):
                continue
            tiles.append(Tile(col=col, row=row, core=core, window=core.buffer(buffer_m)))
    if not tiles:
        raise ValueError(f"{crs_hint} produced no tiles; check tile_size_m against its extent")
    return tiles

shared_edges

shared_edges(tile: Tile, tiles: list[Tile]) -> BaseGeometry | None

The parts of tile's boundary that a neighbour in tiles will emit instead.

Cores partition by area, but they share their boundary lines, and a road running along a seam lies in both. Clipping each tile to its closed core therefore emits that road twice — measured at 30.9 km of output from 24.0 km of input on a grid aligned to the tiling, which is not a contrived case: tiles are axis-aligned to the CRS origin and so are many streets.

Ownership is resolved by giving every shared edge to the lower-indexed tile: each tile drops the linework lying along its right and top edges, because the neighbour on that side keeps it. Deterministic, and independent of the order tiles happen to be processed in.

Source code in src/lczkit/cleaning/tiles.py

def shared_edges(tile: Tile, tiles: list[Tile]) -> BaseGeometry | None:
    """The parts of `tile`'s boundary that a neighbour in `tiles` will emit instead.

    Cores partition by *area*, but they share their boundary *lines*, and a road running along a
    seam lies in both. Clipping each tile to its closed core therefore emits that road twice —
    measured at 30.9 km of output from 24.0 km of input on a grid aligned to the tiling, which
    is not a contrived case: tiles are axis-aligned to the CRS origin and so are many streets.

    Ownership is resolved by giving every shared edge to the lower-indexed tile: each tile drops
    the linework lying along its right and top edges, because the neighbour on that side keeps
    it. Deterministic, and independent of the order tiles happen to be processed in.
    """
    present = {(other.col, other.row) for other in tiles}
    minx, miny, maxx, maxy = tile.core.bounds
    edges: list[BaseGeometry] = []
    if (tile.col + 1, tile.row) in present:
        edges.append(shapely.LineString([(maxx, miny), (maxx, maxy)]))
    if (tile.col, tile.row + 1) in present:
        edges.append(shapely.LineString([(minx, maxy), (maxx, maxy)]))
    return shapely.union_all(edges) if edges else None

layer_extent

layer_extent(*layers: GeoDataFrame) -> Polygon

The bounding box covering every non-empty layer in layers, in their shared CRS.

Used to decide what to tile. Takes the union of bounds rather than of geometry: a tiling only needs to cover the data, and unioning a metropolitan street network to find that out would cost more than the tiling saves.

Source code in src/lczkit/cleaning/tiles.py

def layer_extent(*layers: gpd.GeoDataFrame) -> Polygon:
    """The bounding box covering every non-empty layer in `layers`, in their shared CRS.

    Used to decide what to tile. Takes the union of bounds rather than of geometry: a tiling
    only needs to cover the data, and unioning a metropolitan street network to find that out
    would cost more than the tiling saves.
    """
    boxes = [box(*layer.total_bounds) for layer in layers if len(layer)]
    if not boxes:
        raise ValueError("every layer is empty; nothing to tile")
    return box(*shapely.union_all(boxes).bounds)

subset

subset(layer: GeoDataFrame, window: Polygon) -> GeoDataFrame

Features of layer intersecting window, as a copy with a fresh index.

Uses the spatial index rather than .clip(): geometry must reach a tile whole, because neatnet decides a street's fate from its full shape and a truncated one would be simplified against a shape that does not exist.

The positions are sorted, and that is load-bearing. geopandas documents sindex.query as returning results in no guaranteed order ("often sorted, but there is no guarantee"), and neatnet simplifies a network as a function of the row order it receives — pinned by test_simplification_depends_on_input_row_order. Every layer arrives in one canonical order (overture._canonical_order, sorted by GERS id) so that two runs of a city agree; taking the query result unsorted here would silently undo that on the tiled path, leaving the order a property of the installed GEOS build rather than of the data. It also fed pooled_artifact_threshold, so a change in STRtree traversal would have moved the threshold and with it the tile cache key.

Source code in src/lczkit/cleaning/tiles.py

def subset(layer: gpd.GeoDataFrame, window: Polygon) -> gpd.GeoDataFrame:
    """Features of `layer` intersecting `window`, as a copy with a fresh index.

    Uses the spatial index rather than `.clip()`: geometry must reach a tile *whole*, because
    neatnet decides a street's fate from its full shape and a truncated one would be simplified
    against a shape that does not exist.

    **The positions are sorted, and that is load-bearing.** geopandas documents `sindex.query` as
    returning results in no guaranteed order ("often sorted, but there is no guarantee"), and
    `neatnet` simplifies a network as a function of the row order it receives — pinned by
    `test_simplification_depends_on_input_row_order`. Every layer arrives in one canonical order
    (`overture._canonical_order`, sorted by GERS id) so that two runs of a city agree; taking the
    query result unsorted here would silently undo that on the tiled path, leaving the
    order a property of the installed GEOS build rather than of the data. It also fed
    `pooled_artifact_threshold`, so a change in STRtree traversal would have moved the threshold and
    with it the tile cache key.
    """
    assert_projected_crs(layer, "layer")
    positions = np.sort(layer.sindex.query(window, predicate="intersects"))
    return layer.iloc[positions].reset_index(drop=True).copy()

Geometry, topology, land use, reporting

lczkit.cleaning.geometry

Geometry helpers shared by more than one cleaning operation.

Small by design. Anything here is used by both buildings and topology, and putting it in either would give the two modules a dependency direction neither wants.

largest_part

largest_part(geometry: GeoSeries) -> GeoSeries

Reduce each geometry to its single largest polygon part.

Subtracting something from a footprint can split it in two — a building a road crosses end to end, or one whose neighbour is subtracted through the middle. Keeping both parts would turn one building into two in building_count and give a building_id two rows; keeping the largest keeps the feature a feature. A geometry the subtraction emptied has no parts and comes back missing, for the caller to drop.

Source code in src/lczkit/cleaning/geometry.py

def largest_part(geometry: gpd.GeoSeries) -> gpd.GeoSeries:
    """Reduce each geometry to its single largest polygon part.

    Subtracting something from a footprint can split it in two — a building a road crosses end to
    end, or one whose neighbour is subtracted through the middle. Keeping both parts would turn one
    building into two in `building_count` and give a `building_id` two rows; keeping the largest
    keeps the feature a feature. A geometry the subtraction emptied has no parts and comes back
    missing, for the caller to drop.
    """
    parts, origin = shapely.get_parts(geometry.to_numpy(), return_index=True)
    if len(parts) == 0:
        empty = pd.Series(index=geometry.index, dtype="object")
        return gpd.GeoSeries(empty, crs=geometry.crs)
    winners = pd.DataFrame({"origin": origin, "area": shapely.area(parts)})
    picked = winners.groupby("origin")["area"].idxmax()
    largest = pd.Series(parts[picked.to_numpy()], index=geometry.index[picked.index])
    return gpd.GeoSeries(largest.reindex(geometry.index), crs=geometry.crs)

lczkit.cleaning.topology

Cross-layer topology cleanup, applied to buildings_topo only.

Strictly sequential: the waterbody check runs against buildings after the road-buffer rule, not against an independent snapshot of the layer.

buildings_area takes no part in any of this. Its contract is that footprint area survives, and every operation here removes some.

resolve_buildings_on_streets

resolve_buildings_on_streets(buildings: GeoDataFrame, streets: GeoDataFrame, *, buffer_m: float, overlap_limit: float) -> tuple[GeoDataFrame, CleaningStep]

Drop footprints mostly inside the roadway; trim those merely reaching into it.

The operation this replaces dropped every footprint whose geometry touched a street centreline. That is wrong for the ordinary European perimeter block, which fronts the street and routinely crosses a generalised centreline by a metre or two: measured on the Berlin fixture it removed 439 footprints carrying 22.5% of all building area, at a mean footprint of 1603 m² against 531 m² for those it kept — three times the size, which is the signature of a rule deleting blocks rather than artefacts.

The replacement measures how much of a footprint lies inside a road buffer of half-width buffer_m:

• above overlap_limit, the footprint is mostly roadway and is dropped
• at or below it, the roadway part is subtracted and the rest kept

Choosing the two values. On Berlin at buffer_m=4.0 the overlap fraction spans the whole [0, 1] range with no gap in the counts, but the median footprint falls monotonically across it, from 1652 m² in the lowest decile to 60 m² in the highest, against a fixture-wide median building of 230 m². That collapse is the separator, and overlap_limit=0.5 is where the features being dropped fall below the size of a typical building: it drops 54 footprints and recovers 95% of the lost area. Narrower buffers compress the distribution until nothing is droppable (p95 = 0.46 at 2 m); wider ones swallow the blocks (p90 = 0.98 at 8 m). Both values are configuration, and both were derived from two European fixtures.

Note what the rule does not do. Berlin's fixture is 6105 OpenStreetMap footprints against 88 Microsoft ML, so the high-overlap tail is overwhelmingly OSM kiosks, shelters and garages. This separates small structures standing in the roadway from blocks fronting it; it is not an ML-noise filter.

Source code in src/lczkit/cleaning/topology.py

def resolve_buildings_on_streets(
    buildings: gpd.GeoDataFrame,
    streets: gpd.GeoDataFrame,
    *,
    buffer_m: float,
    overlap_limit: float,
) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Drop footprints mostly inside the roadway; trim those merely reaching into it.

    The operation this replaces dropped every footprint whose geometry touched a street
    *centreline*. That is wrong for the ordinary European perimeter block, which fronts the street
    and routinely crosses a generalised centreline by a metre or two: measured on the Berlin
    fixture it removed 439 footprints carrying **22.5% of all building area**, at a mean footprint
    of 1603 m² against 531 m² for those it kept — three times the size, which is the signature of a
    rule deleting blocks rather than artefacts.

    The replacement measures how much of a footprint lies inside a road buffer of half-width
    `buffer_m`:

    - above `overlap_limit`, the footprint is mostly roadway and is dropped
    - at or below it, the roadway part is subtracted and the rest kept

    **Choosing the two values.** On Berlin at `buffer_m=4.0` the overlap fraction spans the whole
    [0, 1] range with no gap in the counts, but the median footprint falls monotonically across it,
    from 1652 m² in the lowest decile to 60 m² in the highest, against a fixture-wide median
    building of 230 m². That collapse is the separator, and `overlap_limit=0.5` is where the
    features being dropped fall below the size of a typical building: it drops 54 footprints and
    recovers 95% of the lost area. Narrower buffers compress the distribution until nothing is
    droppable (p95 = 0.46 at 2 m); wider ones swallow the blocks (p90 = 0.98 at 8 m). Both values
    are configuration, and both were derived from two European fixtures.

    Note what the rule does *not* do. Berlin's fixture is 6105 OpenStreetMap footprints against 88
    Microsoft ML, so the high-overlap tail is overwhelmingly OSM kiosks, shelters and garages. This
    separates small structures standing in the roadway from blocks fronting it; it is not an
    ML-noise filter.
    """
    assert_projected_crs(buildings, "buildings")
    if buildings.empty or streets.empty:
        return buildings, _passthrough(
            "resolve_buildings_on_streets",
            buildings,
            road_buffer_m=buffer_m,
            overlap_limit=overlap_limit,
        )
    assert_projected_crs(streets, "streets")

    # `road` is positional, one entry per input building, so every use of it below indexes by
    # position rather than by label.
    road = _road_near_each(buildings, streets, buffer_m=buffer_m)
    footprint_area = buildings.geometry.area
    inside = pd.Series(
        shapely.area(shapely.intersection(buildings.geometry.to_numpy(), road)),
        index=buildings.index,
    )
    fraction = inside.div(footprint_area.where(footprint_area > 0)).fillna(0.0)

    dropped = fraction > overlap_limit
    trimmed = (fraction > 0.0) & ~dropped

    kept = buildings.loc[~dropped].copy()
    to_trim = trimmed.loc[kept.index]
    if to_trim.any():
        positions = np.flatnonzero(trimmed.to_numpy())
        kept.loc[to_trim, "geometry"] = largest_part(
            gpd.GeoSeries(
                shapely.difference(buildings.geometry.to_numpy()[positions], road[positions]),
                index=kept.index[to_trim.to_numpy()],
                crs=buildings.crs,
            )
        ).to_numpy()
        kept = kept.loc[kept.geometry.notna() & ~kept.geometry.is_empty]
    kept = kept.reset_index(drop=True)

    step = CleaningStep(
        stage="buildings_topo",
        operation="resolve_buildings_on_streets",
        n_in=len(buildings),
        n_out=len(kept),
        area_in_m2=float(footprint_area.sum()),
        area_out_m2=float(kept.geometry.area.sum()),
        detail={
            "road_buffer_m": buffer_m,
            "overlap_limit": overlap_limit,
            "n_dropped": int(dropped.sum()),
            "n_trimmed": int(trimmed.sum()),
            "area_dropped_m2": float(footprint_area[dropped].sum()),
            "median_dropped_footprint_m2": float(np.nanmedian(footprint_area[dropped]))
            if dropped.any()
            else None,
        },
    )
    return kept, step

drop_buildings_on_waterbodies

drop_buildings_on_waterbodies(buildings: GeoDataFrame, waterbodies: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Drop buildings intersecting waterbodies.

Still a plain intersection test, unlike the street rule: a footprint reaching into a river is a conflation error rather than a building fronting it, and the measured cost is 6 features and 0.4% of area on Berlin rather than 439 and 22.5%.

Source code in src/lczkit/cleaning/topology.py

def drop_buildings_on_waterbodies(
    buildings: gpd.GeoDataFrame, waterbodies: gpd.GeoDataFrame
) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Drop buildings intersecting `waterbodies`.

    Still a plain intersection test, unlike the street rule: a footprint reaching into a river is a
    conflation error rather than a building fronting it, and the measured cost is 6 features and
    0.4% of area on Berlin rather than 439 and 22.5%.
    """
    return _drop_intersecting(
        buildings,
        waterbodies,
        stage="buildings_topo",
        operation="drop_buildings_on_waterbodies",
        areal=True,
    )

drop_waterlines_through_buildings

drop_waterlines_through_buildings(waterlines: GeoDataFrame, buildings: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Drop waterlines that pass through buildings.

Source code in src/lczkit/cleaning/topology.py

def drop_waterlines_through_buildings(
    waterlines: gpd.GeoDataFrame, buildings: gpd.GeoDataFrame
) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Drop waterlines that pass through `buildings`."""
    return _drop_intersecting(
        waterlines,
        buildings,
        stage="topology",
        operation="drop_waterlines_through_buildings",
        areal=False,
    )

apply_cross_layer_topology

apply_cross_layer_topology(buildings_topo: GeoDataFrame, streets: GeoDataFrame, waterlines: GeoDataFrame, waterbodies: GeoDataFrame, *, road_buffer_m: float, road_overlap_limit: float) -> CrossLayerResult

Run cross-layer topology cleanup on buildings_topo.

streets and waterbodies pass through unchanged — only the topological building layer and waterlines are filtered.

Source code in src/lczkit/cleaning/topology.py

def apply_cross_layer_topology(
    buildings_topo: gpd.GeoDataFrame,
    streets: gpd.GeoDataFrame,
    waterlines: gpd.GeoDataFrame,
    waterbodies: gpd.GeoDataFrame,
    *,
    road_buffer_m: float,
    road_overlap_limit: float,
) -> CrossLayerResult:
    """Run cross-layer topology cleanup on `buildings_topo`.

    `streets` and `waterbodies` pass through unchanged — only the topological building layer and
    `waterlines` are filtered.
    """
    cleaned, street_step = resolve_buildings_on_streets(
        buildings_topo, streets, buffer_m=road_buffer_m, overlap_limit=road_overlap_limit
    )
    cleaned, water_step = drop_buildings_on_waterbodies(cleaned, waterbodies)
    cleaned_waterlines, waterline_step = drop_waterlines_through_buildings(waterlines, cleaned)
    return (
        cleaned,
        streets,
        cleaned_waterlines,
        waterbodies,
        [street_step, water_step, waterline_step],
    )

lczkit.cleaning.land_use

Land-use cleaning — deliberately minimal.

Land use is functional metadata passed through to industrial_fraction, not a morphological layer. It gets no thresholded cleaning of any kind: no size filter, no overlap resolution, no absorption of small parcels. The one operation applied is geometry repair, so that an area overlay cannot die on an invalid upstream polygon.

MultiPolygons are kept rather than exploded — a land-use parcel is legitimately multipart, and area overlays handle multipart geometry without help.

clean_land_use

clean_land_use(land_use: GeoDataFrame) -> tuple[GeoDataFrame, CleaningStep]

Repair invalid land-use geometries via make_valid(). Never changes feature count.

Source code in src/lczkit/cleaning/land_use.py

def clean_land_use(land_use: gpd.GeoDataFrame) -> tuple[gpd.GeoDataFrame, CleaningStep]:
    """Repair invalid land-use geometries via `make_valid()`. Never changes feature count."""
    assert_projected_crs(land_use, "land_use")
    n_invalid = int((~land_use.geometry.is_valid).sum())
    fixed = land_use.copy()
    fixed["geometry"] = fixed.geometry.make_valid()
    step = CleaningStep(
        stage="land_use",
        operation="fix_invalid_geometries",
        n_in=len(land_use),
        n_out=len(fixed),
        detail={"n_invalid_before": n_invalid},
    )
    return fixed, step

lczkit.cleaning.report

Structured record of what every cleaning operation consumed and produced.

Every cleaning function returns the layer it produced alongside one or more CleaningStep fragments; the pipeline orchestrator in pipeline.py is the only place these fragments are assembled into a CleaningReport. This keeps every cleaning function a pure, independently testable transform with no shared mutable state.

Steps record area, not only feature counts. Counts alone are why a 23.5% loss of Berlin's building footprint area went unnoticed for a long time: the two operations responsible removed 1177 and 439 features respectively, and by count the second looks like the smaller of the two. By area the first costs 0.12% and the second 22.5%. Building surface fraction carries roughly 47% of the classification metric, so footprint area is the output here and a report that does not state it cannot be audited.

Stage module-attribute

Stage = Literal['ingestion', 'buildings', 'buildings_area', 'buildings_topo', 'streets', 'land_use', 'topology']

The pipeline stages a step can belong to.

"ingestion" covers repairs made to raw source data before any layer-specific cleaning starts — at present only the clipping of features that have no finite representation in the study area's UTM zone. It spans every layer at once, which is why it is not any one of the others.

CleaningStep

Bases: BaseModel

One recorded operation: how many features and how much area went in, and came out.

area_in_m2 / area_out_m2 are polygon area in the layer's projected CRS. They are 0.0 for stages whose geometry has no area — the street and waterline steps — which is a true statement about a linework layer rather than a missing measurement.

area_retained property

area_retained: float | None

Fraction of incoming area this step passed through, or None for an areal-free stage.

Can exceed 1.0: trimming a self-overlapping footprint set removes double-counted area, so a step that lowers the measured total can be the one making it correct.

FootprintCoverage

Bases: BaseModel

How much ground the raw footprints cover, against how much area they sum to.

The two differ because sources self-overlap. Overture conflates footprints from OSM, Esri, Google Open Buildings and Microsoft ML geometry-first and winner-takes-all, and where that leaves a podium and its tower as two features, or a duplicate that survived conflation, their areas add up while the ground beneath them is counted once. Measured on the committed fixtures: Kowloon's raw footprints double-count 7.52% of their summed area against Berlin's 0.61%.

This is why retention is measured against the union. Building surface fraction sums overlay pieces, so the union is what BSF is trying to measure, and a criterion stated against the sum is unmeetable wherever self-overlap exceeds the tolerance: trim_overlaps takes Kowloon to 98.40% of summed area without dropping a single feature, which reads as 1.6% attrition and is in fact 1.6% of double-counting removed. The sum-based criterion and "trim overlaps but do not merge" were jointly unsatisfiable there.

raw_self_overlap_fraction is reported in its own right: it is a real source-quality signal, and a city where it is high is a city where footprint provenance deserves a second look.

area_summed_m2 class-attribute instance-attribute

area_summed_m2: float = 0.0

Summed area of buildings_area, the layer every area statistic reads.

raw_self_overlap_fraction property

raw_self_overlap_fraction: float | None

Share of the raw summed area that is double-counted. 0.0 for a disjoint source.

residual_self_overlap_fraction property

residual_self_overlap_fraction: float | None

The same measure on buildings_area, after cleaning.

Non-zero means the BSF numerator still double-counts, because building surface fraction sums overlay pieces. trim_overlaps resolves pairwise overlaps and does not claim to resolve stacks, so this is expected to be small but not always zero, and it is reported rather than assumed away.

union_retention property

union_retention: float | None

buildings_area's summed area over the raw union — the retention criterion.

Above 1.0 means the layer holds more area than the ground it covers, i.e. residual double-counting; below 0.99 means ground was lost. The criterion is one-sided on the losing side, and the excess is reported through residual_self_overlap_fraction rather than being folded into a single number that cannot distinguish the two failures.

ground_retention property

ground_retention: float | None

Union out over union in: the share of covered ground the layer kept.

Immune to double-counting on both sides, so it is the honest "did cleaning lose any building?" figure, where union_retention is the one the BSF numerator cares about.

CleaningReport

Bases: BaseModel

The full sequence of cleaning steps applied by one clean_vectors() run.

footprints class-attribute instance-attribute

footprints: FootprintCoverage | None = None

Union-based footprint accounting. None where the run predates it or recorded no buildings.

stage_steps

stage_steps(stage: Stage) -> list[CleaningStep]

Every step recorded for stage, in the order it ran.

Source code in src/lczkit/cleaning/report.py

def stage_steps(self, stage: Stage) -> list[CleaningStep]:
    """Every step recorded for `stage`, in the order it ran."""
    return [step for step in self.steps if step.stage == stage]

area_retention

area_retention(stage: Stage) -> float | None

Area out of stage's last step over area into its first — end-to-end, not per-step.

Measured against the sum of incoming area, which is the wrong denominator for the buildings_area acceptance criterion and is retained because it is the per-stage figure every other stage wants. For buildings_area, read footprints.union_retention instead — see FootprintCoverage for why the sum cannot serve.

Source code in src/lczkit/cleaning/report.py

def area_retention(self, stage: Stage) -> float | None:
    """Area out of `stage`'s last step over area into its first — end-to-end, not per-step.

    Measured against the **sum** of incoming area, which is the wrong denominator for the
    `buildings_area` acceptance criterion and is retained because it is the per-stage figure
    every other stage wants. For `buildings_area`, read `footprints.union_retention` instead —
    see `FootprintCoverage` for why the sum cannot serve.
    """
    steps = [step for step in self.stage_steps(stage) if step.area_in_m2 > 0.0]
    if not steps:
        return None
    return steps[-1].area_out_m2 / steps[0].area_in_m2