Height cascade

The building-height cascade: fill every footprint's height, and say how well.

Overture solves footprint coverage; it does not solve height. The answer here is a graded cascade — per-building height, height_source and height_confidence — plus per-unit height_completeness and the full tier distribution, so that "90% surveyed heights" and "90% coarse raster fallback" are distinguishable in the output. They produce the same LCZ label with very different trustworthiness.

This is what differentiates lczkit from an implementation running on OpenStreetMap alone. Overture Maps solves footprint coverage; it does not solve height. Overture merges several sources winner-takes-all per building, and only OpenStreetMap among them carries heights — so wherever a machine-learning footprint source won the geometry, there is no height at all. That is much of the Global South, and plenty of developed cities outside the centre: Cairo, Nairobi and Islamabad each carry a directly measured height for about 1% of their building area.

The answer is a graded cascade — a series of sources tried in order, each filling only what the last left empty — plus honest reporting of which one answered. The tiers, in order:

1. Overture height, else num_floors × storey_height
2. Google Open Buildings 2.5D — retired from the default cascade, measured harmful
3. WSF-3D, global ~90 m
4. GHS-BUILT-H, global 100 m

Tiers 2–4 are areal products: a raster giving one value per cell, so every building inside a cell is assigned the same neighbourhood average. That is a categorically weaker measurement than tier 1 and the output says so, per building via height_source and per unit via height_tier_fractions — because "90% real heights" and "90% coarse raster fallback" produce the same label with very different trustworthiness.

Per-building accuracy is the wrong acceptance test for a height product

Open Buildings 2.5D has the lowest per-building error of the three and the only within-unit skill, and still makes the map worse. Hr is a geometric mean, and dispersion depresses it: GOB's within-unit spread is 0.441 against reality's 0.195, so over half is noise. Evaluate any new tier on within-unit dispersion against reality, not on MAE.

lczkit.heights.cascade

Running the height tiers in order, and recording what each one managed.

The report is the point. A cascade that fills every building from surveyed heights and one that fell all the way through to a 100 m raster produce the same labels, so the output has to tell them apart.

UNRESOLVED module-attribute

UNRESOLVED = 'unresolved'

height_source tag for a building no tier could resolve. Its height stays null.

HeightTierResult

Bases: BaseModel

What one tier contributed to one cascade run.

n_candidates instance-attribute

n_candidates: int

Buildings still unresolved when this tier ran.

n_filled instance-attribute

n_filled: int

Buildings this tier resolved.

filled_by_source class-attribute instance-attribute

filled_by_source: dict[str, int] = Field(default_factory=dict)

Breakdown by height_source tag, for tiers that resolve rows by more than one route.

HeightFillReport

Bases: BaseModel

The full record of one fill_heights() run, for the output manifest.

height_sources class-attribute instance-attribute

height_sources: list[str] = Field(default_factory=list)

Every tag the configured cascade could emit, in tier order, ending with "unresolved". Downstream stages use this as the fixed column set for per-unit tier fractions, so the output schema depends on the configured cascade rather than on which tiers happened to fire.

cascade_height_sources

cascade_height_sources(tiers: Sequence[HeightSource]) -> list[str]

Every height_source tag tiers can emit, in order, plus "unresolved".

Source code in src/lczkit/heights/cascade.py

def cascade_height_sources(tiers: Sequence[HeightSource]) -> list[str]:
    """Every `height_source` tag `tiers` can emit, in order, plus `"unresolved"`."""
    return [source for tier in tiers for source in tier.height_sources] + [UNRESOLVED]

fill_heights

fill_heights(buildings: GeoDataFrame, tiers: Sequence[HeightSource]) -> tuple[GeoDataFrame, HeightFillReport]

Run tiers in order over buildings, returning the filled layer and a report.

Every returned row carries height_source and height_confidence. Rows no tier resolved are tagged "unresolved" with a null height rather than raising: with a tier's product simply absent from disk — the normal state of tiers 2-4 — refusing to return would make the package unusable where it should instead be honest. Callers that need completeness should read it off the report or off height_completeness per unit.

buildings is not mutated.

Source code in src/lczkit/heights/cascade.py

def fill_heights(
    buildings: gpd.GeoDataFrame, tiers: Sequence[HeightSource]
) -> tuple[gpd.GeoDataFrame, HeightFillReport]:
    """Run `tiers` in order over `buildings`, returning the filled layer and a report.

    Every returned row carries `height_source` and `height_confidence`. Rows no tier resolved
    are tagged `"unresolved"` with a null `height` rather than raising: with a tier's product
    simply absent from disk — the normal state of tiers 2-4 — refusing to return would make the
    package unusable where it should instead be honest. Callers that need completeness should
    read it off the report or off `height_completeness` per unit.

    `buildings` is not mutated.
    """
    out = prepare(buildings)
    results: list[HeightTierResult] = []

    for tier in tiers:
        pending = out["height_source"].isna()
        out = tier.fill(out)
        newly = pending & out["height_source"].notna()
        results.append(
            HeightTierResult(
                tier=tier.name,
                n_candidates=int(pending.sum()),
                n_filled=int(newly.sum()),
                filled_by_source={
                    str(source): int(count)
                    for source, count in out.loc[newly, "height_source"].value_counts().items()
                },
            )
        )

    unresolved = out["height_source"].isna()
    out.loc[unresolved, "height"] = np.nan
    out.loc[unresolved, "height_source"] = UNRESOLVED

    report = HeightFillReport(
        n_buildings=len(out),
        n_resolved=int(len(out) - unresolved.sum()),
        n_unresolved=int(unresolved.sum()),
        tiers=results,
        height_sources=cascade_height_sources(tiers),
    )
    return out, report

Tiers

lczkit.heights.tiers

The tiers of the height cascade.

Two HeightSource implementations cover all four tiers: OvertureAttributeTier is tier 1, and ArealRasterTier is tiers 2, 3 and 4 — Google Open Buildings 2.5D, WSF-3D and GHS-BUILT-H differ only in which file they read and how its values scale, which is why they are three configured instances of one class rather than three classes. Adding a fifth areal product is an entry in HeightConfig.areal_tiers, not a new implementation.

Tiers 2-4 are areal products: they assign a neighbourhood mean to individual buildings. That is a categorically weaker measurement than tier 1, and the per-building height_source and per-unit tier fractions exist so the output says so rather than presenting one number as if it were the other.

Every tier claims only rows no earlier tier tagged, so the cascade order in HeightConfig.areal_tiers is what makes the result correct. Run them through lczkit.heights.cascade.fill_heights rather than calling fill by hand.

OVERTURE_HEIGHT module-attribute

OVERTURE_HEIGHT = 'overture_height'

height_source tag for a height taken from Overture's own height attribute.

OVERTURE_NUM_FLOORS module-attribute

OVERTURE_NUM_FLOORS = 'overture_num_floors'

height_source tag for a height derived from num_floors x storey_height_m.

OvertureAttributeTier

OvertureAttributeTier(*, storey_height_m: float, height_confidence: float | None, num_floors_confidence: float | None)

Tier 1: Overture's height, else num_floors x storey_height_m.

The strongest tier available, and still not a guarantee of a surveyed measurement — a quarter of the tier-1 heights in the Berlin test fixture are Microsoft ML values conflated onto OSM footprints (see lczkit.heights.provenance). Where Overture attaches its own per-building confidence to a height, that real number is written to height_confidence in preference to the configured one, so the distinction survives into the output.

The num_floors path is a weaker measurement than the height path — storey height varies regionally and is a real error source — which is why it gets its own height_source tag and its own confidence rather than being folded into the first.

Configure the storey height and the two confidences this tier's routes carry.

Both confidences are required rather than defaulted: they are an ordinal ranking of measurement quality with no published value behind them, and a default would write a quality claim nobody chose into every manifest. _require_confidence says so if unset.

Source code in src/lczkit/heights/tiers.py

def __init__(
    self,
    *,
    storey_height_m: float,
    height_confidence: float | None,
    num_floors_confidence: float | None,
) -> None:
    """Configure the storey height and the two confidences this tier's routes carry.

    Both confidences are required rather than defaulted: they are an ordinal ranking of
    measurement quality with no published value behind them, and a default would write a
    quality claim nobody chose into every manifest. `_require_confidence` says so if unset.
    """
    if storey_height_m <= 0:
        raise ValueError(f"storey_height_m must be positive, got {storey_height_m}")
    self.storey_height_m = storey_height_m
    self.height_confidence = _require_confidence(
        height_confidence, "overture_height_confidence"
    )
    self.num_floors_confidence = _require_confidence(
        num_floors_confidence, "overture_num_floors_confidence"
    )

name property

name: str

The tier's name in the cascade. Its rows are tagged by height_sources, not by this.

height_sources property

height_sources: tuple[str, ...]

The height_source tags this tier can write.

Two tags, not one — the height and num_floors routes are different measurements and the output must keep them apart. See the class docstring.

fill

fill(buildings: GeoDataFrame) -> GeoDataFrame

Resolve heights from Overture's own attributes, height first and num_floors after.

A non-positive or non-finite height is cleared rather than kept: it is not a measurement, and leaving it would hand the next tier a resolved row and the output a zero-height building. Where Overture attaches its own per-building confidence to a height, that real number is written in preference to the configured one.

Rows this tier cannot resolve are returned unchanged for the next tier in the cascade. buildings is not mutated.

Source code in src/lczkit/heights/tiers.py

def fill(self, buildings: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Resolve heights from Overture's own attributes, `height` first and `num_floors` after.

    A non-positive or non-finite `height` is cleared rather than kept: it is not a
    measurement, and leaving it would hand the next tier a resolved row and the output a
    zero-height building. Where Overture attaches its own per-building confidence to a
    height, that real number is written in preference to the configured one.

    Rows this tier cannot resolve are returned unchanged for the next tier in the cascade.
    `buildings` is not mutated.
    """
    out = prepare(buildings)
    todo = out["height_source"].isna()

    # A non-positive or non-finite height is not a measurement; clear it so a later tier
    # sees an unresolved row rather than a zero-height building.
    from_height = todo & out["height"].gt(0) & np.isfinite(out["height"])
    out.loc[todo & ~from_height, "height"] = np.nan

    _, overture_confidence = height_attribution(out)
    out.loc[from_height, "height_source"] = OVERTURE_HEIGHT
    out.loc[from_height, "height_confidence"] = (
        overture_confidence.loc[from_height].fillna(self.height_confidence).astype("float64")
    )

    floors = numeric_column(out, "num_floors")
    from_floors = todo & ~from_height & floors.ge(1).fillna(False)
    out.loc[from_floors, "height"] = floors.loc[from_floors] * self.storey_height_m
    out.loc[from_floors, "height_source"] = OVERTURE_NUM_FLOORS
    out.loc[from_floors, "height_confidence"] = self.num_floors_confidence
    return out

ArealRasterTier

ArealRasterTier(*, name: str, path: Path, confidence: float | None, band: int = 1, scale: float = 1.0, nodata: float | None = None, min_height_m: float = 0.0)

Tiers 2-4: a neighbourhood mean read out of a local height raster.

One class for Google Open Buildings 2.5D, WSF-3D and GHS-BUILT-H, configured per product (lczkit.config.ArealTierConfig). Nothing product-specific is hardcoded — band, unit scale, nodata and minimum valid height are all config, because none of these products is present on the system this was written against and guessing at one is exactly the failure mode that produces a quietly wrong map.

Values at or below min_height_m are read as "no built volume in this cell" and left for the next tier, rather than written out as a zero-height building.

Configure one areal product: where its raster is, and how to read a height out of it.

band, scale, nodata and min_height_m are per-product and come from the product's own documentation via ArealTierConfig — nothing here is hardcoded, because guessing at a unit scale or a nodata value is what produces a quietly wrong map. confidence is required for the same reason the Overture tier's are.

Source code in src/lczkit/heights/tiers.py

def __init__(
    self,
    *,
    name: str,
    path: Path,
    confidence: float | None,
    band: int = 1,
    scale: float = 1.0,
    nodata: float | None = None,
    min_height_m: float = 0.0,
) -> None:
    """Configure one areal product: where its raster is, and how to read a height out of it.

    `band`, `scale`, `nodata` and `min_height_m` are per-product and come from the product's
    own documentation via `ArealTierConfig` — nothing here is hardcoded, because guessing at
    a unit scale or a nodata value is what produces a quietly wrong map. `confidence` is
    required for the same reason the Overture tier's are.
    """
    self.name = name
    self.path = path
    self.confidence = _require_confidence(confidence, f"areal_tiers[{name!r}].confidence")
    self.band = band
    self.scale = scale
    self.nodata = nodata
    self.min_height_m = min_height_m

height_sources property

height_sources: tuple[str, ...]

The single height_source tag this tier writes — its own product name.

from_config classmethod

from_config(config: ArealTierConfig, path: Path) -> ArealRasterTier

Build a tier from its serialised config and an already-resolved raster path.

The split is deliberate: config carries everything reproducible into the manifest, while path is placed per study area by lczkit.sources.height_products. A filename baked into the config would be a filename for one city.

Source code in src/lczkit/heights/tiers.py

@classmethod
def from_config(cls, config: ArealTierConfig, path: Path) -> ArealRasterTier:
    """Build a tier from its serialised config and an already-resolved raster path.

    The split is deliberate: `config` carries everything reproducible into the manifest,
    while `path` is placed per study area by `lczkit.sources.height_products`. A filename
    baked into the config would be a filename for one city.
    """
    return cls(
        name=config.name,
        path=path,
        confidence=config.confidence,
        band=config.band,
        scale=config.scale,
        nodata=config.nodata,
        min_height_m=config.min_height_m,
    )

fill

fill(buildings: GeoDataFrame) -> GeoDataFrame

Give each unresolved building the mean raster height under its own footprint.

A categorically weaker measurement than tier 1 and tagged as such: the product's cell is 90-100 m, so what a building receives is its neighbourhood's mean height, not its own. Values at or below min_height_m read as "no built volume in this cell" and are left for the next tier rather than written out as a zero-height building.

Rows this tier cannot resolve are returned unchanged. buildings is not mutated.

Source code in src/lczkit/heights/tiers.py

def fill(self, buildings: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Give each unresolved building the mean raster height under its own footprint.

    A categorically weaker measurement than tier 1 and tagged as such: the product's cell is
    90-100 m, so what a building receives is its neighbourhood's mean height, not its own.
    Values at or below `min_height_m` read as "no built volume in this cell" and are left for
    the next tier rather than written out as a zero-height building.

    Rows this tier cannot resolve are returned unchanged. `buildings` is not mutated.
    """
    out = prepare(buildings)
    assert_projected_crs(out, "buildings")
    todo = out["height_source"].isna()
    if not todo.any():
        return out

    sampled = zonal_mean(
        self.path,
        out.loc[todo].geometry,
        band=self.band,
        nodata=self.nodata,
    )
    heights = sampled * self.scale
    usable = np.isfinite(heights) & (heights > self.min_height_m)
    resolved = out.index[todo][usable]

    out.loc[resolved, "height"] = heights[usable]
    out.loc[resolved, "height_source"] = self.name
    out.loc[resolved, "height_confidence"] = self.confidence
    return out

numeric_column

numeric_column(buildings: GeoDataFrame, name: str) -> Series

buildings[name] coerced to float, or an all-null float Series if the column is absent.

Overture's height is float and num_floors a nullable integer, and either can be missing entirely from a non-Overture layer. Every read of both goes through here so that "column absent" and "column present but null" behave identically. Neither is an error at this layer: a missing height is what the rest of the cascade exists to answer.

Source code in src/lczkit/heights/tiers.py

def numeric_column(buildings: gpd.GeoDataFrame, name: str) -> pd.Series:
    """`buildings[name]` coerced to float, or an all-null float Series if the column is absent.

    Overture's `height` is float and `num_floors` a nullable integer, and either can be missing
    entirely from a non-Overture layer. Every read of both goes through here so that "column
    absent" and "column present but null" behave identically. Neither is an error at this layer:
    a missing height is what the rest of the cascade exists to answer.
    """
    if name not in buildings.columns:
        return pd.Series(np.nan, index=buildings.index, dtype="float64")
    return pd.to_numeric(buildings[name], errors="coerce").astype("float64")

prepare

prepare(buildings: GeoDataFrame) -> GeoDataFrame

Copy buildings with the cascade's three output columns present and correctly typed.

Idempotent, so a tier run standalone gets the same frame shape the cascade would give it. A row is "unresolved" exactly when height_source is null — not when height is null, because tier 1 reads a height that is already there and must still tag its provenance.

Source code in src/lczkit/heights/tiers.py

def prepare(buildings: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Copy `buildings` with the cascade's three output columns present and correctly typed.

    Idempotent, so a tier run standalone gets the same frame shape the cascade would give it.
    A row is "unresolved" exactly when `height_source` is null — not when `height` is null,
    because tier 1 reads a `height` that is already there and must still tag its provenance.
    """
    out = buildings.copy()
    out["height"] = numeric_column(out, "height")
    if "height_source" not in out.columns:
        out["height_source"] = pd.Series(pd.NA, index=out.index, dtype="object")
    if "height_confidence" not in out.columns:
        out["height_confidence"] = pd.Series(np.nan, index=out.index, dtype="float64")
    return out

build_cascade

build_cascade(config: HeightConfig, source_dir: Callable[[str], Path]) -> list[HeightSource]

Assemble the ordered tier list from config.

source_dir is Settings.source_dir in a real run; taking the callable rather than Settings keeps tests able to point tiers at a temporary directory without DATA_DIR.

An areal tier with no filename is skipped — the product is simply not available, and a shorter cascade with an honest height_completeness beats a failure. A tier that names a file which is not there raises, because that is a misconfiguration rather than an absence. A tier with enabled=False is skipped whether or not its file is there: the flag means off, and a flag whose meaning depended on which function read it would be worse than no flag.

Source code in src/lczkit/heights/tiers.py

def build_cascade(config: HeightConfig, source_dir: Callable[[str], Path]) -> list[HeightSource]:
    """Assemble the ordered tier list from config.

    `source_dir` is `Settings.source_dir` in a real run; taking the callable rather than
    `Settings` keeps tests able to point tiers at a temporary directory without `DATA_DIR`.

    An areal tier with no `filename` is skipped — the product is simply not available, and a
    shorter cascade with an honest `height_completeness` beats a failure. A tier that *names* a
    file which is not there raises, because that is a misconfiguration rather than an absence.
    A tier with `enabled=False` is skipped whether or not its file is there: the flag means off,
    and a flag whose meaning depended on which function read it would be worse than no flag.
    """
    tiers: list[HeightSource] = [
        OvertureAttributeTier(
            storey_height_m=config.storey_height_m,
            height_confidence=config.overture_height_confidence,
            num_floors_confidence=config.overture_num_floors_confidence,
        )
    ]
    seen = {OVERTURE_HEIGHT, OVERTURE_NUM_FLOORS}
    for tier_config in config.areal_tiers:
        if tier_config.name in seen:
            raise ValueError(f"duplicate height tier name: {tier_config.name!r}")
        seen.add(tier_config.name)
        if tier_config.filename is None or not tier_config.enabled:
            continue
        path = source_dir(tier_config.source_dir_name) / tier_config.filename
        if not path.is_file():
            raise FileNotFoundError(
                f"Height tier {tier_config.name!r} is configured to read {path}, which does not "
                "exist. Place the product there, or clear its `filename` to skip the tier."
            )
        tiers.append(ArealRasterTier.from_config(tier_config, path))
    return tiers

Raster reads

lczkit.heights.raster

A minimal local zonal read for the height cascade.

Tiers 2-4 need raster values, and the RasterSource protocol is about land-cover fractions per unit rather than a mean per footprint. This module is deliberately kept to one function, so a tier could be repointed at a RasterSource by changing one line.

It is not a general zonal-statistics implementation and should not grow into one — lczkit.landcover brings exactextract for that. What it does is read one window of one band and reduce it to a mean per footprint, which is all tiers 2-4 ask for.

zonal_mean

zonal_mean(path: Path, geoms: GeoSeries, *, band: int = 1, nodata: float | None = None) -> ndarray

Mean raster value under each geometry in geoms, as a float array aligned to geoms.

Returns nan for any geometry the raster cannot answer for: outside its extent, or covering only nodata cells. That is a real answer — "this product does not know" — and the caller passes those rows on to the next tier rather than inventing a value.

Cells are attributed with all_touched=True, so a footprint smaller than one cell still picks up the cell(s) it overlaps. For the coarse products this backs (~90-100 m), a building is typically far smaller than a cell and the result is that cell's neighbourhood mean, which is exactly what those products measure. Footprints that burn no cell at all fall back to the value at their representative point.

geoms is reprojected to the raster's CRS internally; it must have a CRS set, and so must the raster. Overlapping geometries are resolved last-writer-wins by rasterio.features; planar enforcement in cleaning means building footprints do not overlap, so this does not arise for the cascade's own use.

Source code in src/lczkit/heights/raster.py

def zonal_mean(
    path: Path,
    geoms: gpd.GeoSeries,
    *,
    band: int = 1,
    nodata: float | None = None,
) -> np.ndarray:
    """Mean raster value under each geometry in `geoms`, as a float array aligned to `geoms`.

    Returns `nan` for any geometry the raster cannot answer for: outside its extent, or covering
    only nodata cells. That is a real answer — "this product does not know" — and the caller
    passes those rows on to the next tier rather than inventing a value.

    Cells are attributed with `all_touched=True`, so a footprint smaller than one cell still
    picks up the cell(s) it overlaps. For the coarse products this backs (~90-100 m), a building
    is typically far smaller than a cell and the result is that cell's neighbourhood mean, which
    is exactly what those products measure. Footprints that burn no cell at all fall back to the
    value at their representative point.

    `geoms` is reprojected to the raster's CRS internally; it must have a CRS set, and so must
    the raster. Overlapping geometries are resolved last-writer-wins by `rasterio.features`;
    planar enforcement in cleaning means building footprints do not overlap, so this does not
    arise for the cascade's own use.
    """
    n = len(geoms)
    empty = np.full(n, np.nan, dtype="float64")
    if n == 0:
        return empty
    if geoms.crs is None:
        raise ValueError("geoms has no CRS set; cannot reproject to the raster's CRS.")

    with rasterio.open(path) as src:
        if src.crs is None:
            raise ValueError(f"{path} declares no CRS; cannot align it to the buildings layer.")
        projected = geoms.to_crs(CRS.from_user_input(src.crs.to_wkt()))
        window = covering_window(src, projected.total_bounds)
        if window is None:
            return empty

        # Read unmasked and build the validity mask here, rather than via `masked=True`: it
        # keeps the file's declared nodata and the caller's override going through one code
        # path, and sidesteps rasterio's masked-read, which trips a NumPy 2.5 deprecation.
        values = src.read(band, window=window).astype("float64")
        declared_nodata = src.nodatavals[band - 1]
        win_transform = src.window_transform(window)

    valid = np.isfinite(values)
    for sentinel in (declared_nodata, nodata):
        if sentinel is not None:
            valid &= values != sentinel

    # Label 0 is the "no footprint here" background, so labels are 1-based and the bincount
    # results are sliced from index 1.
    labels = features.rasterize(
        ((geom, index) for index, geom in enumerate(projected, start=1) if geom is not None),
        out_shape=values.shape,
        transform=win_transform,
        fill=0,
        all_touched=True,
        dtype="int32",
    )

    flat_labels = labels[valid]
    counts = np.bincount(flat_labels, minlength=n + 1)[1:]
    sums = np.bincount(flat_labels, weights=values[valid], minlength=n + 1)[1:]
    means = np.where(counts > 0, sums / np.maximum(counts, 1), np.nan)

    unburnt = np.flatnonzero(counts == 0)
    if unburnt.size:
        means[unburnt] = _sample_representative_points(
            gpd.GeoSeries(projected.iloc[unburnt]), values, valid, win_transform
        )
    return means

Completeness and provenance

height_completeness and height_tier_fractions are primary deliverables, not diagnostics.

lczkit.heights.completeness

Per-unit height provenance: height_completeness and the full tier distribution.

These are primary deliverables, not diagnostics. "90% surveyed heights" and "90% coarse raster fallback" produce the same LCZ label with very different trustworthiness, so the output carries the whole distribution across tiers rather than a single completeness number.

TIER1_SOURCES module-attribute

TIER1_SOURCES = (OVERTURE_HEIGHT, OVERTURE_NUM_FLOORS)

The height_source tags that count towards height_completeness.

Both tier-1 routes: a surveyed height, and a storey count multiplied by a storey height. The two are reported as separate fractions as well, so a stricter reading — completeness as surveyed heights only — stays computable from the same table without this module choosing it for everyone.

height_metrics

height_metrics(buildings: GeoDataFrame, units: GeoDataFrame, sources: Sequence[str]) -> DataFrame

Per-unit height provenance, indexed by unit_id to match units.

Returns height_completeness — the area fraction of building footprint resolved by tier 1 — plus one height_frac_<source> column for every tag in sources, which is HeightFillReport.height_sources in a real run. Passing the tag list explicitly is what fixes the output schema to the configured cascade rather than to whichever tiers happened to fire in a given city; Phases 6 and 7 need that stability.

Weighting is by the area of each footprint inside the unit, so a building straddling two grid cells contributes to both in proportion. For EnclosureUnits this is equivalent to assigning each building to one unit, since cross-layer topology cleaning already removes buildings that intersect the streets forming enclosure boundaries.

Units containing no building area are all-null, not zero: "no buildings here" and "0% tier-1 coverage" are different statements and collapsing them would misreport every park and water body as a height-data failure.

Source code in src/lczkit/heights/completeness.py

def height_metrics(
    buildings: gpd.GeoDataFrame,
    units: gpd.GeoDataFrame,
    sources: Sequence[str],
) -> pd.DataFrame:
    """Per-unit height provenance, indexed by `unit_id` to match `units`.

    Returns `height_completeness` — the area fraction of building footprint resolved by tier 1 —
    plus one `height_frac_<source>` column for every tag in `sources`, which is
    `HeightFillReport.height_sources` in a real run. Passing the tag list explicitly is what
    fixes the output schema to the configured cascade rather than to whichever tiers happened to
    fire in a given city; Phases 6 and 7 need that stability.

    Weighting is by the area of each footprint *inside* the unit, so a building straddling two
    grid cells contributes to both in proportion. For `EnclosureUnits` this is equivalent to
    assigning each building to one unit, since cross-layer topology cleaning already removes
    buildings that intersect the streets forming enclosure boundaries.

    Units containing no building area are all-null, not zero: "no buildings here" and "0% tier-1
    coverage" are different statements and collapsing them would misreport every park and water
    body as a height-data failure.
    """
    assert_projected_crs(units, "units")
    if units.index.name != "unit_id":
        raise ValueError("units must be indexed by unit_id")
    columns = [f"{FRACTION_PREFIX}{source}" for source in sources]

    if buildings.empty:
        return _all_null(units.index, columns)
    assert_projected_crs(buildings, "buildings")
    if buildings.crs != units.crs:
        raise ValueError(f"buildings.crs ({buildings.crs}) != units.crs ({units.crs})")
    if "height_source" not in buildings.columns:
        raise ValueError(
            "buildings has no height_source column; run lczkit.heights.cascade.fill_heights "
            "before computing per-unit height metrics."
        )

    pieces = gpd.overlay(
        units[["geometry"]].reset_index(),
        buildings[["height_source", "geometry"]].reset_index(drop=True),
        how="intersection",
    )
    if pieces.empty:
        return _all_null(units.index, columns)

    pieces["piece_area"] = pieces.geometry.area
    by_source = (
        pieces.groupby(["unit_id", "height_source"], observed=True)["piece_area"]
        .sum()
        .unstack(fill_value=0)
        .reindex(columns=list(sources), fill_value=0.0)
        .reindex(index=units.index)
    )
    totals = by_source.sum(axis=1)
    fractions = by_source.div(totals.where(totals > 0), axis=0)
    fractions.columns = pd.Index(columns)

    tier1 = [f"{FRACTION_PREFIX}{source}" for source in sources if source in TIER1_SOURCES]
    completeness = fractions[tier1].sum(axis=1) if tier1 else pd.Series(0.0, index=units.index)
    fractions.insert(0, "height_completeness", completeness.where(totals > 0))
    fractions.index.name = "unit_id"
    return fractions

Dispersion

Coverage is only half of what a substituted height does. Hr is a geometric mean, so it is depressed by spread and rises as spread collapses — and the tiers that shipped compress within-unit spread rather than inflating it, which is the opposite of the failure Open Buildings was rejected for. Median coefficient of variation across whole-city runs: 0.266 for real Overture heights in Berlin, 0.192 for WSF-3D in Nairobi, and 0.112 for GHS-BUILT-H in Bogotá, where 23.6% of units carry a single height throughout. Each run reports its own figures in the manifest.

lczkit.heights.dispersion

Within-unit height dispersion, per cascade tier — what an areal product costs Hr.

height_completeness says where a height came from. It says nothing about what the substitution did to the shape of the height distribution inside a unit, and that is the quantity Hr is sensitive to: it is a geometric mean, so it is depressed by spread and rises as spread collapses.

The sensitivity was established from one side. Google Open Buildings 2.5D had the lowest per-building error of any tier and the only within-unit skill, and it degraded the map, because its within-unit spread was 0.441 against reality's 0.195 — over half of it noise. A height tier is therefore accepted on within-unit dispersion and not on mean absolute error.

This module measures the other side, which nothing has: the tiers that were adopted compress dispersion rather than inflating it. Measured on the runs on disk, over units with buildings:

dominant source	city	median h_std	median CV	constant units
Overture height	Berlin	1.52 m	0.266	0.1%
WSF-3D	Nairobi	0.88 m	0.192	1.3%
WSF-3D	Bogota	1.05 m	0.207	1.1%
GHS-BUILT-H	Bogota	0.36 m	0.112	23.6%

A 90 m or 100 m product hands one height to every building it covers, so what survives inside a unit is variation between raster cells rather than between buildings. Same mechanism, opposite sign, and it biases Hr upward exactly where the cascade is doing the most work.

Reported per run because it is the target any future shrinkage work aims at: shrinking a fine-resolution product toward the unit mean is only worth doing against a measured statement of how much dispersion the incumbent has already lost.

TierDispersion

Bases: BaseModel

Within-unit height spread across the units one tier dominates.

source instance-attribute

source: str

The height_source tag, e.g. "wsf3d".

n_units instance-attribute

n_units: int

Units where this tier supplied more building area than any other.

median_h_std instance-attribute

median_h_std: float | None

Median of h_std — the area-weighted standard deviation of building height within a unit.

median_cv instance-attribute

median_cv: float | None

Median of h_std / h_mean_area_weighted. The scale-free form, and the one comparable against the 0.441-against-0.195 figures above.

constant_fraction instance-attribute

constant_fraction: float

Share of those units whose buildings all carry the same height to within a centimetre. A direct reading of how often the product resolves nothing inside a unit at all.

DispersionReport

Bases: BaseModel

Within-unit height dispersion for one run, per tier.

min_building_surface_fraction instance-attribute

min_building_surface_fraction: float

Units below this are excluded: a unit holding almost no building has a spread that is about its two buildings rather than about its fabric.

min_building_count instance-attribute

min_building_count: int

Units with fewer buildings are excluded, for the same reason. A spread over two buildings is not a description of a neighbourhood.

n_units instance-attribute

n_units: int

Units that passed both filters and carried a dispersion value.

dispersion_report

dispersion_report(parameters: DataFrame, *, min_building_surface_fraction: float = 0.05, min_building_count: int = 3) -> DispersionReport

Within-unit height dispersion per tier, from a finished parameter table.

parameters is what lczkit.ucp.compute_parameters() returns joined to the per-unit height fractions — the table a run assembles anyway — so this reads columns rather than recomputing anything, and it moves no measurement.

A unit is attributed to whichever tier supplied the largest share of its building area, which is a simplification and is stated as one: a unit split evenly between Overture and WSF-3D is counted wholly against the larger. The alternative, area-weighting every unit into every tier, would mix distributions and defeat the comparison the table exists to make.

Source code in src/lczkit/heights/dispersion.py

def dispersion_report(
    parameters: pd.DataFrame,
    *,
    min_building_surface_fraction: float = 0.05,
    min_building_count: int = 3,
) -> DispersionReport:
    """Within-unit height dispersion per tier, from a finished parameter table.

    `parameters` is what `lczkit.ucp.compute_parameters()` returns joined to the per-unit height
    fractions — the table a run assembles anyway — so this reads columns rather than
    recomputing anything, and it moves no measurement.

    A unit is attributed to whichever tier supplied the largest share of its building area, which
    is a simplification and is stated as one: a unit split evenly between Overture and WSF-3D is
    counted wholly against the larger. The alternative, area-weighting every unit into every tier,
    would mix distributions and defeat the comparison the table exists to make.
    """
    fractions = [column for column in parameters.columns if column.startswith(FRACTION_PREFIX)]
    empty = DispersionReport(
        min_building_surface_fraction=min_building_surface_fraction,
        min_building_count=min_building_count,
        n_units=0,
    )
    required = {"h_std", "h_mean_area_weighted", "building_surface_fraction", "building_count"}
    if not fractions or not required <= set(parameters.columns):
        return empty

    usable = parameters[
        (parameters["building_surface_fraction"] >= min_building_surface_fraction)
        & (parameters["building_count"] >= min_building_count)
        & parameters["h_std"].notna()
        & parameters[fractions].notna().any(axis=1)
    ]
    if usable.empty:
        return empty

    dominant = (
        usable[fractions]
        .fillna(0.0)
        .idxmax(axis=1)
        .astype("string")
        .str.removeprefix(FRACTION_PREFIX)
    )
    mean = usable["h_mean_area_weighted"]
    cv = usable["h_std"].div(mean.where(mean > 0))

    tiers: list[TierDispersion] = []
    for source, rows in usable.groupby(dominant.to_numpy(), sort=True):
        spread = rows["h_std"]
        tiers.append(
            TierDispersion(
                source=str(source),
                n_units=int(len(rows)),
                median_h_std=_finite(spread.median()),
                median_cv=_finite(cv.loc[rows.index].median()),
                constant_fraction=float((spread < 0.01).mean()),
            )
        )
    return DispersionReport(
        min_building_surface_fraction=min_building_surface_fraction,
        min_building_count=min_building_count,
        n_units=int(len(usable)),
        tiers=tiers,
    )

lczkit.heights.provenance

Reading per-attribute provenance out of Overture's sources column.

Overture records provenance per attribute, not just per feature. Each entry in sources carries a property naming what it applies to — '' for the footprint as a whole, or a JSON pointer such as /properties/height for one attribute — alongside the dataset it came from and, for machine-derived values, a confidence.

That distinction is load-bearing here, and it qualifies the usual description of Overture conflation as winner-takes-all with no attribute fusion. In release 2026-07-22.0, 394 of the 6195 footprints in the Berlin test fixture are OSM-won yet carry a height attributed to Microsoft ML Buildings — a quarter of every tier-1 height in that extent. Grouping provenance by the footprint's dataset alone would report those as surveyed OSM heights.

These two functions are the single parser for that column; the tier-1 height source and the source-availability diagnostic both read through them rather than each growing their own.

FOOTPRINT_PROPERTY module-attribute

FOOTPRINT_PROPERTY = ''

sources[].property value marking the entry that describes the feature as a whole.

HEIGHT_PROPERTY module-attribute

HEIGHT_PROPERTY = '/properties/height'

sources[].property value marking the entry that describes the height attribute.

footprint_datasets

footprint_datasets(buildings: GeoDataFrame) -> Series

The upstream dataset that won each footprint's geometry, as a string Series.

Null where the row carries no whole-feature provenance entry, and for every row when buildings has no sources column at all — a non-Overture VectorSource degrades to "no provenance known" rather than raising.

Source code in src/lczkit/heights/provenance.py

def footprint_datasets(buildings: gpd.GeoDataFrame) -> pd.Series:
    """The upstream dataset that won each footprint's geometry, as a string Series.

    Null where the row carries no whole-feature provenance entry, and for every row when
    `buildings` has no `sources` column at all — a non-Overture `VectorSource` degrades to "no
    provenance known" rather than raising.
    """
    if "sources" not in buildings.columns:
        return pd.Series(pd.NA, index=buildings.index, dtype="object")
    return buildings["sources"].map(
        lambda s: (_entry_for(s, FOOTPRINT_PROPERTY) or {}).get("dataset")
    )

height_attribution

height_attribution(buildings: GeoDataFrame) -> tuple[Series, Series]

(dataset, confidence) describing where each row's height value came from.

Where Overture records a /properties/height provenance entry, both come from it — this is the case that separates a conflated machine-learning height from a surveyed one, and the confidence is a real per-building number rather than anything this package assigns.

Where it does not, the dataset falls back to the footprint's own and the confidence is null: the height, if any, came in with the footprint, and Overture attaches no confidence to it.

Source code in src/lczkit/heights/provenance.py

def height_attribution(buildings: gpd.GeoDataFrame) -> tuple[pd.Series, pd.Series]:
    """`(dataset, confidence)` describing where each row's `height` value came from.

    Where Overture records a `/properties/height` provenance entry, both come from it — this is
    the case that separates a conflated machine-learning height from a surveyed one, and the
    `confidence` is a real per-building number rather than anything this package assigns.

    Where it does not, the dataset falls back to the footprint's own and the confidence is null:
    the height, if any, came in with the footprint, and Overture attaches no confidence to it.
    """
    footprint = footprint_datasets(buildings)
    if "sources" not in buildings.columns:
        return footprint, pd.Series(np.nan, index=buildings.index, dtype="float64")

    entries = buildings["sources"].map(lambda s: _entry_for(s, HEIGHT_PROPERTY))
    dataset = entries.map(lambda e: (e or {}).get("dataset")).fillna(footprint)
    confidence = pd.to_numeric(
        entries.map(lambda e: (e or {}).get("confidence")), errors="coerce"
    ).astype("float64")
    return dataset, confidence

lczkit.heights.diagnostic

Source-availability diagnostic: is this city viable before anyone waits for a full run?

It reports non-null height and num_floors counts grouped by Overture source dataset over the study area — twice, because the data forces it: once by the dataset that won the footprint, and once by the dataset that supplied the height.

Those two tables disagree, and the disagreement is the useful part. In the Berlin test fixture, grouping by footprint reports 1509 heights on OpenStreetMap footprints — but 394 of those heights are Microsoft ML Buildings values conflated onto an OSM footprint. Read the first table alone and a quarter of the city's tier-1 heights look surveyed when they are predicted.

UNKNOWN_DATASET module-attribute

UNKNOWN_DATASET = '(unknown)'

Stand-in for a row carrying no usable provenance, so it is counted rather than dropped.

DatasetAvailability

Bases: BaseModel

Height and floor-count availability for one upstream dataset.

HeightProvenance

Bases: BaseModel

How many heights one upstream dataset actually supplied.

SourceAvailability

Bases: BaseModel

The full diagnostic for one study area, for the output manifest.

by_footprint_dataset class-attribute instance-attribute

by_footprint_dataset: list[DatasetAvailability] = Field(default_factory=list)

Grouped by the dataset that won each footprint's geometry.

by_height_dataset class-attribute instance-attribute

by_height_dataset: list[HeightProvenance] = Field(default_factory=list)

Heights only, grouped by the dataset that supplied the height value rather than the one that won the footprint.

Deliberately narrower than the table above. Provenance in Overture is per attribute, and the only attribute conflated away from its footprint is height; reporting a building's num_floors under the dataset that supplied its height would credit a dataset with a value it never provided. Rows with no height contribute nothing here, which is why these counts do not sum to n_buildings.

source_availability

source_availability(buildings: GeoDataFrame) -> SourceAvailability

Count height and floor availability by upstream dataset over buildings.

Reads Overture's sources column. A layer without one still returns valid totals, with every row grouped under "(unknown)" — the diagnostic degrades rather than blocking a non-Overture VectorSource.

Source code in src/lczkit/heights/diagnostic.py

def source_availability(buildings: gpd.GeoDataFrame) -> SourceAvailability:
    """Count height and floor availability by upstream dataset over `buildings`.

    Reads Overture's `sources` column. A layer without one still returns valid totals, with
    every row grouped under `"(unknown)"` — the diagnostic degrades rather than blocking a
    non-Overture `VectorSource`.
    """
    height = numeric_column(buildings, "height")
    floors = numeric_column(buildings, "num_floors")
    has_height = height.notna() & height.gt(0)
    has_floors = floors.notna() & floors.ge(1)

    footprint = footprint_datasets(buildings)
    height_dataset, _ = height_attribution(buildings)

    return SourceAvailability(
        n_buildings=len(buildings),
        n_with_height=int(has_height.sum()),
        n_with_num_floors=int(has_floors.sum()),
        by_footprint_dataset=_by_footprint(footprint, has_height, has_floors),
        by_height_dataset=_by_height(height_dataset[has_height]),
    )

lczkit.heights.inherit

Carry resolved heights from one building layer onto another.

Cleaning produces two building layers from one source. The height cascade runs once, on buildings_area — that is the complete population, and so the honest denominator for height_completeness and for the source-availability diagnostic. buildings_topo still needs heights, because momepy.street_profile measures a canyon's height-to-width ratio against the buildings walling it, so the values have to reach it from somewhere.

Not by building_id. A dissolved buildings_topo feature keeps one arbitrary constituent's id, and on the Berlin fixture that constituent is as likely to be an absorbed sub-20 m² shed as the block that absorbed it — which would hand a perimeter block a garage's height. Largest overlap is the correct rule and costs one overlay.

INHERITED_COLUMNS module-attribute

INHERITED_COLUMNS = ('height', 'height_source', 'height_confidence')

What the cascade produces per building, and therefore what travels.

inherit_heights

inherit_heights(target: GeoDataFrame, source: GeoDataFrame) -> GeoDataFrame

Give every target footprint the height of the source footprint it overlaps most.

source must have been through lczkit.heights.cascade.fill_heights(). A target footprint overlapping no source footprint keeps a null height rather than being dropped or imputed: street_profile skips null-height buildings when averaging, so the aspect ratio's coverage degrades rather than its value. Neither input is mutated.

Source code in src/lczkit/heights/inherit.py

def inherit_heights(target: gpd.GeoDataFrame, source: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Give every `target` footprint the height of the `source` footprint it overlaps most.

    `source` must have been through `lczkit.heights.cascade.fill_heights()`. A `target` footprint
    overlapping no `source` footprint keeps a null height rather than being dropped or imputed:
    `street_profile` skips null-height buildings when averaging, so the aspect ratio's *coverage*
    degrades rather than its value. Neither input is mutated.
    """
    assert_projected_crs(target, "target")
    assert_projected_crs(source, "source")
    if target.crs != source.crs:
        raise ValueError(f"target.crs ({target.crs}) != source.crs ({source.crs})")
    missing = [column for column in INHERITED_COLUMNS if column not in source.columns]
    if missing:
        raise ValueError(
            f"source has no {', '.join(missing)}; run lczkit.heights.cascade.fill_heights on it "
            "before inheriting heights from it."
        )

    out = target.drop(columns=list(INHERITED_COLUMNS), errors="ignore").copy()
    if out.empty or source.empty:
        return out.assign(**dict.fromkeys(INHERITED_COLUMNS, None))

    pieces = gpd.overlay(
        out[["geometry"]].reset_index(names="_target"),
        source[[*INHERITED_COLUMNS, "geometry"]].reset_index(drop=True),
        how="intersection",
        keep_geom_type=True,
    )
    if pieces.empty:
        return out.assign(**dict.fromkeys(INHERITED_COLUMNS, None))

    # `idxmax` then `.loc`, rather than a groupby aggregation: `first`/`last` skip nulls per
    # column, so a target whose largest overlap is an unresolved building would silently inherit
    # a smaller neighbour's height instead of staying null. That is imputation, which the cascade
    # refuses — an unresolved height is information about the data, not a gap to be filled.
    pieces["_overlap"] = pieces.geometry.area
    winners = pieces.loc[pieces.groupby("_target")["_overlap"].idxmax()].set_index("_target")
    for column in INHERITED_COLUMNS:
        out[column] = pd.Series(winners[column]).reindex(out.index)
    return out