Pipeline

The end-to-end chain: clean vectors, generate units, fill heights, read land cover, compute parameters, classify, write outputs. run_pipeline is the only place the stages are wired together, and the command line calls it rather than restating any of it.

lczkit.pipeline

The whole chain, from a bbox to a run directory and optionally a map site.

run_pipeline is the only place the stages are wired together. The command line calls it rather than restating any of it, so there is one definition of what a run does.

What is deliberately absent. Validation. write_run takes a validation= report and the manifest has a slot for it, but the chain never populates them: agreement is measured separately, against reference datasets that are not always on disk. Wiring it in here would make every run depend on those datasets being present. Call lczkit.validation yourself when you have them.

StageObserver is how a caller watches a long run without this module choosing a rendering. The command line passes one backed by rich; any object with the same two methods will do.

STAGES module-attribute

STAGES = ('clean_vectors', 'heights', 'units', 'land_cover', 'provenance', 'parameters', 'classify', 'write_run', 'build_site')

Stage names, in order, so a caller can size a progress display before the run starts.

StageObserver

Bases: Protocol

Something that watches each stage begin and end.

A typing.Protocol rather than a base class, following the same decision the five data-source protocols were built under: the point is the seam, and a caller that already has a timer should not have to inherit anything to use it.

stage

stage(name: str) -> AbstractContextManager[None]

A context manager wrapping one stage's work.

Source code in src/lczkit/pipeline.py

def stage(self, name: str) -> AbstractContextManager[None]:
    """A context manager wrapping one stage's work."""
    ...

PipelineResult dataclass

PipelineResult(outputs: RunOutputs, site: SiteReport | None, site_skipped: str | None = None, stages: dict[str, float] = dict(), height_products: dict[str, str | None] = dict())

What a run produced, and how long each stage took.

site instance-attribute

site: SiteReport | None

None when the run was asked not to build one, or when tippecanoe is absent.

site_skipped class-attribute instance-attribute

site_skipped: str | None = None

Why no site was built, where one was asked for. None when one was built or not wanted.

The site is the last stage and everything else is already on disk by the time it runs, so a missing tippecanoe must not cost a caller the run. It used to: the error propagated out of run_pipeline, the command line turned it into an exit code, and the line naming the run directory was never printed — a ten-minute city reported as a failure with no mention that its output existed. lczkit site build <run_dir> completes it later.

stages class-attribute instance-attribute

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

Wall seconds per stage, in the order they ran.

height_products class-attribute instance-attribute

height_products: dict[str, str | None] = field(default_factory=dict)

Which areal height product file each enabled tier resolved to, by tier name.

None where the product has no coverage for this extent — Open Buildings stops at Europe — which is a different state from a tier that was disabled, and stays separable here.

run_dir property

run_dir: Path

Where everything was written.

seconds property

seconds: float

Total wall time across the stages that ran.

build_strategy

build_strategy(config: UnitsConfig, *, buildings: GeoDataFrame | None = None) -> SpatialUnitStrategy

The configured SpatialUnitStrategy.

buildings is only read by patch, and only when patch_merge_on_morphology is on. It is passed at construction rather than to generate because the protocol's signature is (bbox, barriers), and widening that for one strategy would put a building layer into an interface the other two have no use for.

Source code in src/lczkit/pipeline.py

def build_strategy(
    config: UnitsConfig, *, buildings: gpd.GeoDataFrame | None = None
) -> SpatialUnitStrategy:
    """The configured `SpatialUnitStrategy`.

    `buildings` is only read by `patch`, and only when `patch_merge_on_morphology` is on. It is
    passed at construction rather than to `generate` because the protocol's signature is
    `(bbox, barriers)`, and widening that for one strategy would put a building layer into an
    interface the other two have no use for.
    """
    if config.strategy == "grid":
        return GridUnits(cell_size_m=config.cell_size_m)
    if config.strategy == "enclosure":
        return EnclosureUnits()
    return PatchUnits(
        min_area_m2=config.patch_min_area_m2,
        max_area_m2=config.patch_max_area_m2,
        buildings=buildings if config.patch_merge_on_morphology else None,
    )

run_pipeline

run_pipeline(settings: Settings, bbox: BBox, *, build_site_after: bool = True, observer: StageObserver | None = None, extent: ExtentRecord | None = None) -> PipelineResult

Clean, fill heights, classify, write a run directory and optionally build its map site.

settings must already carry a runnable configuration — CleaningConfig and HeightConfig both have fields that default to None and raise at call time. lczkit.presets.apply_preset is what fills them.

Every path comes from settings: the run directory, the tile cache, and the input/ subdirectories the Overture and height-product fetchers own. Nothing existing under input/ is modified or removed.

extent records how bbox was chosen — a named place, a So2Sat window, or four numbers — and goes into the manifest. It defaults to the bbox alone, which is all a library caller who computed their own window can honestly claim.

Source code in src/lczkit/pipeline.py

def run_pipeline(
    settings: Settings,
    bbox: BBox,
    *,
    build_site_after: bool = True,
    observer: StageObserver | None = None,
    extent: ExtentRecord | None = None,
) -> PipelineResult:
    """Clean, fill heights, classify, write a run directory and optionally build its map site.

    `settings` must already carry a runnable configuration — `CleaningConfig` and `HeightConfig`
    both have fields that default to `None` and raise at call time. `lczkit.presets.apply_preset`
    is what fills them.

    Every path comes from `settings`: the run directory, the tile cache, and the `input/`
    subdirectories the Overture and height-product fetchers own. Nothing existing under `input/`
    is modified or removed.

    `extent` records **how** `bbox` was chosen — a named place, a So2Sat window, or four numbers —
    and goes into the manifest. It defaults to the bbox alone, which is all a library caller who
    computed their own window can honestly claim.
    """
    watch = observer if observer is not None else _NullObserver()
    covered = extent if extent is not None else ExtentRecord(kind="bbox", bbox=bbox)
    stages: dict[str, float] = {}

    @contextmanager
    def timed(name: str) -> Iterator[None]:
        """Run one stage under the observer, recording its wall time into `stages`."""
        started = time.perf_counter()
        with watch.stage(name):
            yield
        stages[name] = time.perf_counter() - started

    source = OvertureSource(settings)
    with timed("clean_vectors"):
        cleaned = clean_vectors(
            source,
            bbox,
            settings.cleaning,
            cache_dir=settings.tile_cache_dir,
        )

    with timed("heights"):
        # Places the products the configured cascade needs, and returns the config with each
        # tier's file resolved. Without this step `build_cascade` finds every areal tier's
        # `filename` unset and silently runs tier 1 alone, so the default cascade needs a step
        # that actually fetches.
        heights, placed = resolve_areal_tiers(settings, bbox)
        tiers = build_cascade(heights, settings.source_dir)
        buildings_area, height_fill = fill_heights(cleaned.buildings_area, tiers)
        buildings_topo = inherit_heights(cleaned.buildings_topo, buildings_area)
        availability = source_availability(cleaned.buildings_area)
        tags = tag_availability(cleaned.buildings_area, cleaned.land_use)

    with timed("units"):
        # The strategy is config, so the chain has to assemble barriers for the two that need
        # them rather than defaulting to the grid and being unable to reach anything else.
        strategy = build_strategy(settings.units, buildings=buildings_area)
        barriers = None
        measure_on_enclosures = settings.ucp.measure_on == "enclosures"
        if settings.units.strategy != "grid" or measure_on_enclosures:
            # `clean_vectors` does not carry rail — it is a barrier layer, not something the
            # cleaning pipeline has a rule for — so it comes straight off the source.
            streets = (
                filter_street_barriers(cleaned.streets)
                if settings.units.drop_pedestrian_barriers
                else cleaned.streets
            )
            barriers = assemble_barriers(
                streets, cleaned.waterbodies, rail=source.rail(bbox).to_crs(cleaned.crs)
            )
        units = strategy.generate(bbox, barriers)

        # A street canyon has to be measured against streets, and a grid cell is not bounded by
        # any. Off by default — see `UcpConfig.measure_on` — and where the target units *are* the
        # enclosures there is nothing to transfer, so the extra partition is skipped.
        measurement_units = units
        if measure_on_enclosures and settings.units.strategy != "enclosure":
            measurement_units = EnclosureUnits().generate(bbox, barriers)

    with timed("land_cover"):
        # `clip_worldcover` resolves the tiles the bbox actually spans and mosaics them. A single
        # hardcoded tile is correct for one city and a 0x0 window — `RasterioIOError` — for the
        # next one, or worse, a quarter of the map silently missing.
        worldcover = clip_worldcover(bbox, settings.run_dir / "worldcover.tif")
        raster = LocalRasterSource(
            settings.land_cover.dataset(settings.ucp.land_cover_dataset), worldcover
        )
        fractions = raster.fractions(units)
        # The surface fractions have to describe the units the parameters are measured on, or the
        # building share and the impervious share it is subtracted from would come from different
        # ground. A second zonal pass, and only when the two unit sets actually differ.
        measurement_fractions = (
            fractions if measurement_units is units else raster.fractions(measurement_units)
        )

    with timed("provenance"):
        # The column set comes from the configured cascade, not from which tiers happened to fire,
        # so a run with no areal product still reports its tier fractions as zeros rather than
        # omitting the columns and changing the output schema.
        provenance = height_metrics(buildings_area, units, cascade_height_sources(tiers))

    with timed("parameters"):
        parameters = compute_parameters(
            measurement_units,
            buildings_area,
            buildings_topo,
            cleaned.streets,
            cleaned.land_use,
            measurement_fractions,
            config=settings.ucp,
            land_cover_config=settings.land_cover,
        )
        if measurement_units is not units:
            parameters = transfer_parameters(parameters, measurement_units, units)

    with timed("classify"):
        classifier = PrototypeClassifier(config=settings.classification)
        classification = classifier.classify(parameters)
        # Off by default, so this is a no-op that still produces a report — a run has to be able to
        # say the filter did not fire as distinct from never having been configured.
        classification, smoothing = modal_filter(
            units,
            classification,
            enabled=settings.classification.modal_filter,
            min_like_neighbours=settings.classification.modal_filter_min_like_neighbours,
        )

    with timed("write_run"):
        outputs = write_run(
            settings,
            units,
            parameters,
            classification,
            classifier,
            extras=fractions.join(provenance),
            cleaning=cleaned.report,
            extent=covered,
            units_report=getattr(strategy, "report", None),
            height_fill=height_fill,
            height_source_availability=availability,
            tag_availability=tags,
            smoothing=smoothing,
            # The site draws its basemap and its extrusions from these, so that an archived run
            # directory rebuilds its own map with no access to `input/`.
            layers={
                "streets": cleaned.streets,
                "water": cleaned.waterbodies,
                "land_use": cleaned.land_use,
                "buildings": buildings_area,
            },
        )

    site: SiteReport | None = None
    skipped: str | None = None
    if build_site_after:
        with timed("build_site"):
            try:
                site = build_site(outputs.run_dir, config=settings.viz)
            except TippecanoeMissingError as error:
                # Caught rather than raised, and only this one: it is a statement about the
                # machine rather than about the run, and everything the run produced is already
                # written. Any other failure here is a defect and stays loud.
                skipped = str(error)

    return PipelineResult(
        outputs=outputs,
        site=site,
        site_skipped=skipped,
        stages=stages,
        height_products=dict(placed),
    )

Raster windows

Clipping a raster to a study window, and checking that what came back actually covers it. The check matters because neither half fails loudly on its own: read(window=…) returns a smaller array rather than raising when the window runs off the edge of a raster, and units the raster never reached come back as all-NaN fractions by design. Together that is a partly missing map with nothing raised.

lczkit.raster_window

Raster windowing shared by the places in the package that read or clip a raster.

The height cascade reads a mean per building footprint; the land-cover sources read class fractions per spatial unit. Those are different reductions over different libraries, but both begin by finding the one window of a raster that covers a set of geometries, and getting the edge padding wrong is a quiet off-by-one rather than a crash. It lives here, next to crs.py, because it belongs to neither phase.

clip_raster and coverage_shortfall join them because a run has to materialise a window before it can reduce over one — the global land-cover and reference products are read remotely and written into the run directory. See lczkit.sources.worldcover.

covering_window

covering_window(src: DatasetReader, bounds: ndarray | tuple[float, float, float, float]) -> Window | None

The raster window covering bounds, padded one cell and clipped to the raster.

bounds is (minx, miny, maxx, maxy) in the raster's own CRS. Returns None when it falls entirely outside the raster, which callers report as "this product cannot answer" rather than as an error.

The one-cell pad matters at the edges: a geometry whose boundary sits exactly on a cell boundary still touches the cell beyond it, and that cell has to be inside the window to be counted at all.

Source code in src/lczkit/raster_window.py

def covering_window(
    src: rasterio.DatasetReader, bounds: np.ndarray | tuple[float, float, float, float]
) -> windows.Window | None:
    """The raster window covering `bounds`, padded one cell and clipped to the raster.

    `bounds` is `(minx, miny, maxx, maxy)` in the raster's own CRS. Returns `None` when it falls
    entirely outside the raster, which callers report as "this product cannot answer" rather than
    as an error.

    The one-cell pad matters at the edges: a geometry whose boundary sits exactly on a cell
    boundary still touches the cell beyond it, and that cell has to be inside the window to be
    counted at all.
    """
    raw = windows.from_bounds(*bounds, transform=src.transform)
    col_off = max(0, math.floor(raw.col_off) - 1)
    row_off = max(0, math.floor(raw.row_off) - 1)
    col_end = min(src.width, math.ceil(raw.col_off + raw.width) + 1)
    row_end = min(src.height, math.ceil(raw.row_off + raw.height) + 1)
    if col_end <= col_off or row_end <= row_off:
        return None
    return windows.Window(col_off, row_off, col_end - col_off, row_end - row_off)

clip_raster

clip_raster(source: str, destination: Path, bbox: BBox) -> Path

Window source to bbox and write it into the run directory, preserving nodata and CRS.

Source code in src/lczkit/raster_window.py

def clip_raster(source: str, destination: Path, bbox: BBox) -> Path:
    """Window `source` to `bbox` and write it into the run directory, preserving nodata and CRS."""
    with rasterio.open(source) as src:
        window = from_bounds(*bbox, transform=src.transform)
        values = src.read(1, window=window)
        profile = src.profile | {
            "driver": "GTiff",
            "height": values.shape[0],
            "width": values.shape[1],
            "transform": src.window_transform(window),
            "compress": "deflate",
            "tiled": False,
            "count": 1,
        }
    with rasterio.open(destination, "w", **profile) as dst:
        dst.write(values, 1)
    return destination

coverage_shortfall

coverage_shortfall(bounds: tuple[float, float, float, float], bbox: BBox, res: float) -> dict[str, float]

How far a raster's bounds fall short of bbox on each side, in pixels.

Empty when the raster covers the window. A shortfall under one pixel is not reported: a clip lands on cell boundaries, so a fraction of a cell is rounding rather than missing ground.

Separated out because the failure this guards against is silent. clip_raster windows with from_bounds and then read(window=...), which returns a smaller array rather than raising when the window overruns the source, and LocalRasterSource.fractions turns units with no coverage into all-NaN rather than an error. Both behaviours are correct on their own; together they let a raster that covers a quarter of the requested window produce a map with a quarter of its land cover missing and nothing anywhere saying so.

Source code in src/lczkit/raster_window.py

def coverage_shortfall(
    bounds: tuple[float, float, float, float], bbox: BBox, res: float
) -> dict[str, float]:
    """How far a raster's `bounds` fall short of `bbox` on each side, in pixels.

    Empty when the raster covers the window. A shortfall under one pixel is not reported: a clip
    lands on cell boundaries, so a fraction of a cell is rounding rather than missing ground.

    Separated out because the failure this guards against is silent. `clip_raster` windows with
    `from_bounds` and then `read(window=...)`, which **returns a smaller array** rather than
    raising when the window overruns the source, and `LocalRasterSource.fractions` turns units
    with no coverage into all-`NaN` rather than an error. Both behaviours are correct on their
    own; together they let a raster that covers a quarter of the requested window produce a map
    with a quarter of its land cover missing and nothing anywhere saying so.
    """
    left, bottom, right, top = bounds
    minx, miny, maxx, maxy = bbox
    gaps = {
        "west": (left - minx) / res,
        "south": (bottom - miny) / res,
        "east": (maxx - right) / res,
        "north": (maxy - top) / res,
    }
    return {side: round(gap, 3) for side, gap in gaps.items() if gap > 1.0}