Skip to content

Morphometrics

A descriptive output, not classifier input. Majer & Fleischmann (2026) computed 107 2D morphometric attributes — dimensional and shape, spatial distribution and intensity, street descriptors and connectivity — over enclosed tessellation cells (ETCs), and found that the correspondence between 2D morphometrics alone and LCZ types is "selective and inconsistent". This package ports the attribute computation, not the paper's RandomForest/CNN prediction schemes, and none of it feeds PrototypeClassifier: the classifier stays the existing prototype-distance metric over Stewart & Oke's urban canopy parameters.

An ETC is momepy's enclosed_tessellation — one cell per building, generated within street-bounded enclosures — restricted to cells with a parent building. It is not a partition of the bounding box, unlike GridUnits/EnclosureUnits: a cell with no building is dropped, since it carries no morphometric meaning. See Tessellation.

The full attribute-to-momepy-call mapping is transcribed in docs/references/tables/majer_2026_morphometrics_menu.md, the authoritative checklist this package's own test suite parses cell for cell against the registry below.

Output ships as its own run artefact, morphometrics.parquet — and a raster, one band per attribute, area-weighted, if a raster resolution is requested. Neither is joined into units.parquet or the classification table: ETCs are a different, finer-grained unit set than whatever the run's classification units are.

lczkit.morphometrics

2D urban morphometrics after Majer & Fleischmann (2026).

A descriptive output, not classifier input. compute_morphometrics returns 107 primary attributes (dimensional/shape, spatial distribution/intensity, street descriptors/connectivity), computed via momepy 1.0+ over enclosed tessellation cells (ETCs — see lczkit.units.tessellation), plus an opt-in contextual expansion. None of it feeds lczkit.classify.PrototypeClassifier: the paper's own finding is that 2D morphometrics alone predict LCZ "selectively and inconsistently", and this package's classifier stays the existing prototype-distance metric over Stewart & Oke's urban canopy parameters. Nor does it replicate the paper's RandomForest/CNN prediction schemes (S1-S4) — those, and satellite imagery, are out of scope for this feature.

Output ships as its own run artefact (morphometrics.parquet, and morphometrics.tif if a raster resolution is requested) — never joined into units.parquet or the classification table, since ETCs are a different, finer-grained unit set than the pipeline's classification units.

See docs/references/tables/majer_2026_morphometrics_menu.md for the full attribute list.

Deliberately no eager re-exports here. raster.py imports filename constants from lczkit.output.writer, which imports MorphometricsReport from report.py — importing any submodule always runs this file first, so re-exporting raster's names here would make importing lczkit.morphometrics.report alone drag in raster.py, and with it output.writer, while output.writer is itself mid-import. Import submodules directly: lczkit.morphometrics.compute.compute_morphometrics, lczkit.morphometrics.raster.*, lczkit.morphometrics.report.MorphometricsReport.

Computing the attributes

The one function lczkit.pipeline.run_pipeline calls for this stage.

lczkit.morphometrics.compute

compute_morphometrics: the one function lczkit.pipeline calls for this stage.

Builds the tessellation, builds every libpysal.graph.Graph the metric modules need exactly once, runs the three attribute blocks (dimensional, distribution, streets), validates the result against the registry, and optionally runs the opt-in contextual expansion.

BUILDING_DISTANCE_BANDS_M module-attribute 

BUILDING_DISTANCE_BANDS_M = (20.0, 100.0, 200.0)

Distance-band neighbourhoods for building-scale metrics: the paper's scales.

BUILDING_KNN module-attribute 

BUILDING_KNN = (10, 20, 30)

Nearest-neighbour counts for the mean-distance metrics: the paper's scales.

ETC_TOPOLOGICAL_STEPS module-attribute 

ETC_TOPOLOGICAL_STEPS = (1, 2, 3)

Contiguity depths for cell-scale metrics: the paper's scales.

Fixed rather than configurable because the registry names every column after them (..._w100m, ..._knn20, ..._w3steps); a different scale would be a different menu.

compute_morphometrics 

compute_morphometrics(bbox: BBox, buildings: GeoDataFrame, streets: GeoDataFrame, waterbodies: GeoDataFrame, *, config: MorphometricsConfig) -> tuple[GeoDataFrame, MorphometricsReport]

The 107 primary morphometric attributes, indexed by ETC unit_id, alongside geometry.

Plus, if config.contextual, their percentile expansion. buildings should be buildings_area from lczkit.cleaning.pipeline.clean_vectors — the area-preserving layer every other area statistic in this package reads. streets should be the cleaned, simplified network; waterbodies the cleaned waterbody polygons. All three must already share one projected CRS, as clean_vectors returns them.

Source code in src/lczkit/morphometrics/compute.py 
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def compute_morphometrics(
    bbox: BBox,
    buildings: gpd.GeoDataFrame,
    streets: gpd.GeoDataFrame,
    waterbodies: gpd.GeoDataFrame,
    *,
    config: MorphometricsConfig,
) -> tuple[gpd.GeoDataFrame, MorphometricsReport]:
    """The 107 primary morphometric attributes, indexed by ETC `unit_id`, alongside geometry.

    Plus, if `config.contextual`, their percentile expansion. `buildings` should be
    `buildings_area` from `lczkit.cleaning.pipeline.clean_vectors` — the
    area-preserving layer every other area statistic in this package reads. `streets` should be
    the cleaned, simplified network; `waterbodies` the cleaned waterbody polygons. All three must
    already share one projected CRS, as `clean_vectors` returns them.
    """
    barriers = assemble_barriers(streets, waterbodies)
    strategy = TessellationUnits(
        buildings=buildings,
        shrink=config.tessellation_shrink,
        segment=config.tessellation_segment,
        threshold=config.tessellation_threshold,
        n_jobs=config.tessellation_n_jobs,
    )
    etc = strategy.generate(bbox, barriers)
    assert strategy.report is not None  # generate() always sets it before returning

    if len(etc) > config.max_tessellation_cells:
        raise ValueError(
            f"tessellation produced {len(etc)} cells, over the configured ceiling of "
            f"{config.max_tessellation_cells} (MorphometricsConfig.max_tessellation_cells)"
        )

    matched_buildings = buildings_for_etc(buildings, etc)

    building_contiguity = graphs.building_contiguity(matched_buildings)
    building_tessellation_adjacency = graphs.tessellation_adjacency_by_building(
        etc, matched_buildings
    )
    building_distance_bands = {
        f"{int(distance)}m": graphs.building_distance_band(matched_buildings, distance)
        for distance in BUILDING_DISTANCE_BANDS_M
    }
    building_knn = {f"knn{k}": graphs.building_knn(matched_buildings, k) for k in BUILDING_KNN}
    etc_contiguity = graphs.etc_contiguity(etc)
    etc_higher_order = {
        steps: graphs.etc_higher_order(etc_contiguity, steps) for steps in ETC_TOPOLOGICAL_STEPS
    }

    dimensional = dimensional_metrics(
        matched_buildings,
        etc,
        building_contiguity=building_contiguity,
        building_w100m=building_distance_bands["100m"],
        building_w200m=building_distance_bands["200m"],
        etc_w3steps=etc_higher_order[3],
    )
    distribution = distribution_metrics(
        matched_buildings,
        etc,
        streets,
        building_contiguity=building_contiguity,
        building_tessellation_adjacency=building_tessellation_adjacency,
        building_adjacency_neighborhood=building_distance_bands["200m"],
        building_w100m=building_distance_bands["100m"],
        building_w200m=building_distance_bands["200m"],
        building_distance_bands=building_distance_bands,
        building_knn=building_knn,
        etc_higher_order=etc_higher_order,
    )
    street = street_metrics(
        etc,
        streets,
        buildings,
        profile_distance_m=config.street_profile_distance_m,
        profile_tick_length_m=config.street_profile_tick_length_m,
    )

    primary = pd.concat([dimensional, distribution, street], axis=1)
    missing = sorted(set(PARAMETER_COLUMNS) - set(primary.columns))
    extra = sorted(set(primary.columns) - set(PARAMETER_COLUMNS))
    if missing or extra:
        raise RuntimeError(
            "compute_morphometrics output does not match the registry — "
            f"missing: {missing}, unexpected: {extra}"
        )
    primary = primary[list(PARAMETER_COLUMNS)]

    n_contextual = 0
    if config.contextual:
        if len(etc) > config.max_contextual_cells:
            raise ValueError(
                f"tessellation produced {len(etc)} cells, over the configured ceiling of "
                f"{config.max_contextual_cells} for the contextual expansion "
                "(MorphometricsConfig.max_contextual_cells)"
            )
        contextual_graph = graphs.etc_higher_order(etc_contiguity, config.contextual_steps)
        contextual = contextual_expand(
            primary, contextual_graph, quantiles=config.contextual_quantiles
        )
        n_contextual = contextual.shape[1]
        primary = pd.concat([primary, contextual], axis=1)

    result = gpd.GeoDataFrame(primary.join(etc[["geometry"]]), geometry="geometry", crs=etc.crs)
    report = MorphometricsReport(
        tessellation=strategy.report,
        n_primary_attributes=len(PARAMETER_COLUMNS),
        contextual_enabled=config.contextual,
        n_contextual_attributes=n_contextual,
    )
    return result, report

Graphs

Every libpysal.graph.Graph the metric modules need, built once and shared. ETC contiguity is fuzzy, not the plain queen contiguity used elsewhere in this package — momepy's own documented requirement for enclosed_tessellation output, which does not form a precise polygonal coverage.

lczkit.morphometrics.graphs

Every libpysal.graph.Graph this package's morphometrics build, in one place.

Two modules building "the same" neighbourhood graph independently is how a weighted mean and a neighbour count end up describing different neighbourhoods without either call site knowing it. Every graph lczkit.morphometrics uses is constructed here and passed down, never rebuilt inline.

ETC contiguity is fuzzy, not the Graph.build_contiguity(gdf, rook=False) queen contiguity lczkit uses everywhere else (lczkit.units.patches, lczkit.classify.smoothing). This is momepy's own documented requirement, not a new pattern invented here: momepy.enclosed_tessellation's own docstring warns that its output "does not form a precise polygonal coverage" and that a contiguity graph over it should use Graph.build_fuzzy_contiguity(tessellation, buffer=...) instead — a plain queen graph misses edges that fall just short of touching because of the shrink/segment tolerances the tessellation algorithm applies.

FUZZY_CONTIGUITY_BUFFER_M module-attribute 

FUZZY_CONTIGUITY_BUFFER_M = 1e-06

momepy's own suggested tolerance for enclosed_tessellation output, applied here rather than left as a magic number at the one call site that needs it.

building_contiguity 

building_contiguity(buildings: GeoDataFrame) -> Graph

Queen contiguity between buildings.

What momepy.perimeter_wall groups joined structures by, and the contiguity_graph argument momepy.building_adjacency compares against.

Source code in src/lczkit/morphometrics/graphs.py 
27
28
29
30
31
32
33
def building_contiguity(buildings: gpd.GeoDataFrame) -> Graph:
    """Queen contiguity between buildings.

    What `momepy.perimeter_wall` groups joined structures by, and the `contiguity_graph`
    argument `momepy.building_adjacency` compares against.
    """
    return Graph.build_contiguity(buildings, rook=False)

tessellation_adjacency_by_building 

tessellation_adjacency_by_building(etc: GeoDataFrame, buildings: GeoDataFrame) -> Graph

ETC contiguity, relabelled onto the building index — what adjacency means between buildings.

momepy.mean_interbuilding_distance documents its adjacency_graph as "a contiguity graph derived from tessellation cells linked to buildings", and its own worked example passes a Delaunay triangulation. Building-footprint contiguity is not a substitute: footprints in the area-preserving layer mostly do not touch, so that graph is nearly edgeless and the metric has no adjacent pair to measure a distance between.

Measured on the Hong Kong fixture, 5 448 buildings:

building queen contiguity   4 242 edges   2 975 isolates (55%)   metric identically 0.0
ETC contiguity             23 128 edges      12 isolates         median 10.06 m, max 61.7

etc and buildings must be row-aligned, as buildings_for_etc returns them: the cells carry the adjacency and the buildings carry the index the metric's geometry is indexed by.

Source code in src/lczkit/morphometrics/graphs.py 
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def tessellation_adjacency_by_building(etc: gpd.GeoDataFrame, buildings: gpd.GeoDataFrame) -> Graph:
    """ETC contiguity, relabelled onto the building index — what adjacency means between buildings.

    `momepy.mean_interbuilding_distance` documents its `adjacency_graph` as "a contiguity graph
    derived from tessellation cells linked to buildings", and its own worked example passes a
    Delaunay triangulation. Building-footprint contiguity is not a substitute: footprints in the
    area-preserving layer mostly do not touch, so that graph is nearly edgeless and the metric
    has no adjacent pair to measure a distance between.

    Measured on the Hong Kong fixture, 5 448 buildings:

        building queen contiguity   4 242 edges   2 975 isolates (55%)   metric identically 0.0
        ETC contiguity             23 128 edges      12 isolates         median 10.06 m, max 61.7

    `etc` and `buildings` must be row-aligned, as `buildings_for_etc` returns them: the cells
    carry the adjacency and the buildings carry the index the metric's `geometry` is indexed by.
    """
    if len(etc) != len(buildings):
        raise ValueError(
            f"etc ({len(etc)}) and buildings ({len(buildings)}) must be row-aligned; pass the "
            "frame `buildings_for_etc` returned"
        )
    return Graph.build_fuzzy_contiguity(
        etc.set_axis(buildings.index), buffer=FUZZY_CONTIGUITY_BUFFER_M
    )

building_distance_band 

building_distance_band(buildings: GeoDataFrame, distance_m: float) -> Graph

Binary distance-band graph over building centroids at distance_m.

Binary (unweighted 0/1 edges) because every building-scale metric in this module reads counts or unweighted means over the neighbourhood, never a distance-decayed one.

Source code in src/lczkit/morphometrics/graphs.py 
63
64
65
66
67
68
69
def building_distance_band(buildings: gpd.GeoDataFrame, distance_m: float) -> Graph:
    """Binary distance-band graph over building centroids at `distance_m`.

    Binary (unweighted 0/1 edges) because every building-scale metric in this module reads counts
    or unweighted means over the neighbourhood, never a distance-decayed one.
    """
    return Graph.build_distance_band(buildings.geometry.centroid, threshold=distance_m, binary=True)

building_knn 

building_knn(buildings: GeoDataFrame, k: int) -> Graph

K-nearest-neighbour graph over building centroids.

The paper's "10/20/30 nearest neighbours" scale for mean distance to neighbours.

coplanar="clique", not libpysal's "raise" default. Two buildings can share a centroid exactly — concentric footprints, or a footprint and its own courtyard ring — and build_knn then refuses outright rather than degrading, taking the whole stage down with it. Measured on real Overture data: a 9 km² window over central Nairobi holds 8 014 buildings at 8 013 unique centroids, one coplanar pair, and that is enough to raise CoplanarError. The failure is data-dependent and silent until it fires, so it would surface unpredictably part-way through a multi-city extraction rather than up front.

"clique" over "jitter" — the other option that avoids the error — because jitter displaces the duplicates randomly, which would make mean_dist_neighbors_building_knn* differ between two runs over identical input. This package spends real effort on determinism (sorted Overture fetches, sorted tile subsets, pinned thread counts); a random tie-break in a shipped metric would undo that for the sake of a tie that "clique" resolves deterministically and defensibly: two buildings at one point genuinely are each other's nearest neighbour, at distance zero.

Source code in src/lczkit/morphometrics/graphs.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
def building_knn(buildings: gpd.GeoDataFrame, k: int) -> Graph:
    """K-nearest-neighbour graph over building centroids.

    The paper's "10/20/30 nearest neighbours" scale for mean distance to neighbours.

    **`coplanar="clique"`, not libpysal's `"raise"` default.** Two buildings can share a centroid
    exactly — concentric footprints, or a footprint and its own courtyard ring — and `build_knn`
    then refuses outright rather than degrading, taking the whole stage down with it. Measured on
    real Overture data: a 9 km² window over central Nairobi holds 8 014 buildings at 8 013 unique
    centroids, one coplanar pair, and that is enough to raise `CoplanarError`. The failure is
    data-dependent and silent until it fires, so it would surface unpredictably part-way through a
    multi-city extraction rather than up front.

    `"clique"` over `"jitter"` — the other option that avoids the error — because jitter displaces
    the duplicates *randomly*, which would make `mean_dist_neighbors_building_knn*` differ between
    two runs over identical input. This package spends real effort on determinism (sorted Overture
    fetches, sorted tile subsets, pinned thread counts); a random tie-break in a shipped metric
    would undo that for the sake of a tie that `"clique"` resolves deterministically and
    defensibly: two buildings at one point genuinely are each other's nearest neighbour, at
    distance zero.
    """
    return Graph.build_knn(buildings.geometry.centroid, k=k, coplanar="clique")

etc_contiguity 

etc_contiguity(etc: GeoDataFrame) -> Graph

Fuzzy queen contiguity between tessellation cells.

See the module docstring for why fuzzy rather than exact contiguity is required here.

Source code in src/lczkit/morphometrics/graphs.py 
 96
 97
 98
 99
100
101
def etc_contiguity(etc: gpd.GeoDataFrame) -> Graph:
    """Fuzzy queen contiguity between tessellation cells.

    See the module docstring for why fuzzy rather than exact contiguity is required here.
    """
    return Graph.build_fuzzy_contiguity(etc, buffer=FUZZY_CONTIGUITY_BUFFER_M)

etc_higher_order 

etc_higher_order(base: Graph, steps: int) -> Graph

base expanded to every cell reachable within steps topological hops.

base itself is included (lower_order=True) — "within N topological steps" in the paper's own wording, not "exactly N steps away".

Source code in src/lczkit/morphometrics/graphs.py 
104
105
106
107
108
109
110
111
112
def etc_higher_order(base: Graph, steps: int) -> Graph:
    """`base` expanded to every cell reachable within `steps` topological hops.

    `base` itself is included (`lower_order=True`) — "within N topological steps" in the
    paper's own wording, not "exactly N steps away".
    """
    if steps <= 1:
        return base
    return base.higher_order(steps, lower_order=True)

etc_granularity_graph 

etc_granularity_graph(base: Graph) -> Graph

base's 1-step neighbourhood with the focal cell included in its own neighbour set.

momepy.percentile/momepy.weighted_character/momepy.neighbors all read a graph's neighbours as other cells, so "granularity" — the paper's "sum of ETC area within 1 topological step" — needs the one graph in this module where a cell counts itself: without assign_self_weight, Graph.lag(area) would sum a cell's neighbours and silently exclude the cell's own area from a quantity meant to describe the immediate cluster it sits in.

Source code in src/lczkit/morphometrics/graphs.py 
115
116
117
118
119
120
121
122
123
124
def etc_granularity_graph(base: Graph) -> Graph:
    """`base`'s 1-step neighbourhood with the focal cell included in its own neighbour set.

    `momepy.percentile`/`momepy.weighted_character`/`momepy.neighbors` all read a graph's
    neighbours as *other* cells, so "granularity" — the paper's "sum of ETC area within 1
    topological step" — needs the one graph in this module where a cell counts itself: without
    `assign_self_weight`, `Graph.lag(area)` would sum a cell's neighbours and silently exclude the
    cell's own area from a quantity meant to describe the immediate cluster it sits in.
    """
    return base.assign_self_weight(1)

Dimensional & shape

lczkit.morphometrics.dimensional

Dimensional & Shape morphometrics — Majer & Fleischmann (2026), Supplementary A §1.

Every metric here is a direct 1:1 momepy 1.0+ call, computed once for buildings and once for ETCs (enclosed tessellation cells), plus area-weighted variants over the neighbourhoods lczkit.morphometrics.graphs builds. See docs/references/tables/majer_2026_morphometrics_menu.md for the full attribute-to-momepy-call mapping this module implements against.

buildings and etc must already share one index (unit_id) — buildings_for_etc (lczkit.units.tessellation) produces that alignment — since every weighted variant here is a plain Series operation across the two.

dimensional_metrics 

dimensional_metrics(buildings: GeoDataFrame, etc: GeoDataFrame, *, building_contiguity: Graph, building_w100m: Graph, building_w200m: Graph, etc_w3steps: Graph) -> DataFrame

The 62 Dimensional & Shape columns, indexed like buildings/etc (unit_id).

building_contiguity is queen contiguity over buildings — what momepy.perimeter_wall groups joined structures by. building_w100m/building_w200m are distance-band graphs over building centroids; etc_w3steps is the fuzzy-contiguity ETC graph expanded to 3 topological steps (lczkit.morphometrics.graphs.etc_higher_order).

Source code in src/lczkit/morphometrics/dimensional.py 
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
def dimensional_metrics(
    buildings: gpd.GeoDataFrame,
    etc: gpd.GeoDataFrame,
    *,
    building_contiguity: Graph,
    building_w100m: Graph,
    building_w200m: Graph,
    etc_w3steps: Graph,
) -> pd.DataFrame:
    """The 62 Dimensional & Shape columns, indexed like `buildings`/`etc` (`unit_id`).

    `building_contiguity` is queen contiguity over `buildings` — what `momepy.perimeter_wall`
    groups joined structures by. `building_w100m`/`building_w200m` are distance-band graphs over
    building centroids; `etc_w3steps` is the fuzzy-contiguity ETC graph expanded to 3 topological
    steps (`lczkit.morphometrics.graphs.etc_higher_order`).
    """
    area_building = buildings.geometry.area
    area_etc = etc.geometry.area
    columns: dict[str, pd.Series] = {
        "area_building": area_building,
        "area_etc": area_etc,
    }

    courtyard_area = momepy.courtyard_area(buildings)
    courtyard_index = momepy.courtyard_index(buildings, courtyard_area=courtyard_area).clip(
        lower=0.0, upper=1.0
    )
    # `momepy.courtyard_area` calls `shapely.get_exterior_ring`, which is defined only for a
    # single `Polygon` — on a `MultiPolygon` it returns `None`, silently making the "filled"
    # term of the courtyard-area formula 0 and the result `-area` instead of a real value.
    # `buildings_area` can legitimately hold MultiPolygons (`trim_overlaps` can cut a footprint in
    # two when it subtracts a neighbour crossing it), so this is nulled explicitly rather than
    # reported as a nonsensical negative — found on real Nairobi data, 3 of 7 214 buildings.
    multipart = buildings.geometry.geom_type == "MultiPolygon"
    # A hole cut by a footprint nested inside this one is ground another building covers, not a
    # courtyard; `cleaning.buildings.trim_overlaps` flags those footprints.
    if NESTED_FLAG in buildings.columns:
        multipart |= buildings[NESTED_FLAG].fillna(False).astype(bool)
    courtyard_area = courtyard_area.mask(multipart)
    courtyard_index = courtyard_index.mask(multipart)
    columns["courtyard_area_building"] = courtyard_area
    columns["courtyard_index_building"] = courtyard_index

    perimeter_wall = momepy.perimeter_wall(buildings, graph=building_contiguity)
    columns["perimeter_wall_building"] = perimeter_wall
    columns["perimeter_wall_building_w100m"] = momepy.weighted_character(
        perimeter_wall, area_building, building_w100m
    )
    columns["perimeter_wall_building_w200m"] = momepy.weighted_character(
        perimeter_wall, area_building, building_w200m
    )

    for name, fn in WEIGHTED_SHAPE_METRICS:
        building_value = fn(buildings)
        etc_value = fn(etc)
        if name in _BOUNDED_TO_UNIT_INTERVAL:
            building_value = building_value.clip(lower=0.0, upper=1.0)
            etc_value = etc_value.clip(lower=0.0, upper=1.0)
        columns[f"{name}_building"] = building_value
        columns[f"{name}_etc"] = etc_value
        columns[f"{name}_building_w100m"] = momepy.weighted_character(
            building_value, area_building, building_w100m
        )
        columns[f"{name}_building_w200m"] = momepy.weighted_character(
            building_value, area_building, building_w200m
        )
        columns[f"{name}_etc_w3steps"] = momepy.weighted_character(etc_value, area_etc, etc_w3steps)

    return pd.DataFrame(columns, index=buildings.index)

Spatial distribution & intensity

Coverage area ratio is not guaranteed to fall in [0, 1]. A building's footprint area over its own ETC's area is usually near 1, but momepy's tessellation shrink/segment step does not guarantee a cell fully contains the building that seeded it — measured on 6.4% of ETCs in the Hong Kong fixture. Reported as measured, not clipped.

lczkit.morphometrics.distribution

Spatial Distribution & Intensity morphometrics — Majer & Fleischmann (2026), Supplementary A §2.

Buildings and ETCs (enclosed tessellation cells) share one index throughout (unit_id) — see lczkit.units.tessellation.buildings_for_etc. Coverage Area Ratio in particular relies on this directly: it is a plain division of a building's area by its own ETC's area, with no join.

distribution_metrics 

distribution_metrics(buildings: GeoDataFrame, etc: GeoDataFrame, streets: GeoDataFrame, *, building_contiguity: Graph, building_tessellation_adjacency: Graph, building_adjacency_neighborhood: Graph, building_w100m: Graph, building_w200m: Graph, building_distance_bands: dict[str, Graph], building_knn: dict[str, Graph], etc_higher_order: dict[int, Graph]) -> DataFrame

The 23 Spatial Distribution & Intensity columns, indexed like buildings/etc (unit_id).

building_adjacency_neighborhood is the 200 m distance-band graph both building_adjacency and mean_interbuilding_distance read as their neighbourhood extent. The two differ in what they take as adjacency: building_contiguity for the first, as its docstring asks, and building_tessellation_adjacency for the second, which momepy documents as tessellation-cell contiguity and which is the difference between a real distance and a column of zeros. building_distance_bands and building_knn key on the paper's own labels ("20m", "100m", "200m", "knn10", "knn20", "knn30"); etc_higher_order keys on the step count (1, 2, 3), already expanded via lczkit.morphometrics.graphs.etc_higher_order.

Source code in src/lczkit/morphometrics/distribution.py 
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def distribution_metrics(
    buildings: gpd.GeoDataFrame,
    etc: gpd.GeoDataFrame,
    streets: gpd.GeoDataFrame,
    *,
    building_contiguity: Graph,
    building_tessellation_adjacency: Graph,
    building_adjacency_neighborhood: Graph,
    building_w100m: Graph,
    building_w200m: Graph,
    building_distance_bands: dict[str, Graph],
    building_knn: dict[str, Graph],
    etc_higher_order: dict[int, Graph],
) -> pd.DataFrame:
    """The 23 Spatial Distribution & Intensity columns, indexed like `buildings`/`etc` (`unit_id`).

    `building_adjacency_neighborhood` is the 200 m distance-band graph both `building_adjacency`
    and `mean_interbuilding_distance` read as their neighbourhood extent. The two differ in what
    they take as *adjacency*: `building_contiguity` for the first, as its docstring asks, and
    `building_tessellation_adjacency` for the second, which momepy documents as tessellation-cell
    contiguity and which is the difference between a real distance and a column of zeros.
    `building_distance_bands` and `building_knn` key on the paper's own labels (`"20m"`, `"100m"`,
    `"200m"`, `"knn10"`, `"knn20"`, `"knn30"`); `etc_higher_order` keys on the step count
    (`1`, `2`, `3`), already expanded via `lczkit.morphometrics.graphs.etc_higher_order`.
    """
    area_etc = etc.geometry.area
    area_building = buildings.geometry.area
    columns: dict[str, pd.Series] = {
        "building_adjacency_200m": momepy.building_adjacency(
            building_contiguity, building_adjacency_neighborhood
        ),
        # `building_tessellation_adjacency`, not `building_contiguity`. momepy documents this
        # argument as tessellation-cell contiguity linked to buildings, and the distinction is not
        # cosmetic: footprints in the area-preserving layer mostly do not touch, so passing their
        # queen contiguity leaves 55% of buildings isolated and the metric identically 0.0 on
        # every input. `building_adjacency` above keeps the building graph, which is what *its*
        # docstring asks for — one graph, two metrics, and only one of them wanted it.
        "mean_interbuilding_distance_200m": momepy.mean_interbuilding_distance(
            buildings, building_tessellation_adjacency, building_adjacency_neighborhood
        ),
        "shared_walls_building": momepy.shared_walls(buildings),
    }

    street_index = momepy.get_nearest_street(buildings, streets)
    building_orientation = momepy.orientation(buildings)
    street_orientation = momepy.orientation(streets)
    alignment = momepy.street_alignment(building_orientation, street_orientation, street_index)
    columns["street_alignment_building"] = alignment
    columns["street_alignment_building_w100m"] = momepy.weighted_character(
        alignment, area_building, building_w100m
    )
    columns["street_alignment_building_w200m"] = momepy.weighted_character(
        alignment, area_building, building_w200m
    )

    # ETC index is the parent building's, so this is a direct row-aligned division — the numerator
    # and denominator describe the same ground by construction, not a join that could drift.
    # Not clipped at 1.0: `momepy.enclosed_tessellation`'s shrink/segment step does not guarantee
    # a cell fully contains the building that seeded it, and on the Hong Kong fixture 7 of 110
    # ETCs (6.4%) are measurably smaller than their own building — a real property of the
    # tessellation on dense, irregular fabric, reported rather than hidden by a clip.
    coverage_area_ratio = (area_building / area_etc.where(area_etc > 0)).reindex(etc.index)
    columns["coverage_area_ratio_etc"] = coverage_area_ratio
    columns["coverage_area_ratio_etc_w3steps"] = momepy.weighted_character(
        coverage_area_ratio, area_etc, etc_higher_order[3]
    )

    columns["etc_granularity_1step"] = etc_granularity_graph(etc_higher_order[1]).lag(area_etc)

    for label, graph in building_distance_bands.items():
        columns[f"neighbors_building_{label}"] = momepy.neighbors(buildings, graph, weighted=False)
        columns[f"mean_dist_neighbors_building_{label}"] = momepy.neighbor_distance(
            buildings, graph
        )
    for label, graph in building_knn.items():
        columns[f"mean_dist_neighbors_building_{label}"] = momepy.neighbor_distance(
            buildings, graph
        )
    for steps, graph in etc_higher_order.items():
        columns[f"neighbors_etc_{steps}steps"] = momepy.neighbors(etc, graph, weighted=False)
        if steps > 1:
            columns[f"mean_dist_neighbors_etc_{steps}steps"] = momepy.neighbor_distance(etc, graph)

    return pd.DataFrame(columns, index=buildings.index)

Street descriptors & connectivity

5 m and 400 m are network-distance radii, not topological hops. Every connectivity function here reads its radius as a hop count unless an edge-length attribute is named — confirmed against the installed momepy rather than assumed, since 5 or 400 hops would be a nonsensical pairing on any real street network.

lczkit.morphometrics.streets

Street Descriptors & Connectivity morphometrics — Majer & Fleischmann (2026), Supplementary A §3.

Computed per street segment (street_length, street_linearity, street_width*) or per network node (everything else), then joined onto each ETC (enclosed tessellation cell) via momepy.get_nearest_street/momepy.get_nearest_node — the paper's own "linked to nearest street, node and edge" rule (Supplementary A).

5 m / 400 m are network-distance radii, not topological hops. Every connectivity function here (mean_node_degree, node_density, edge_node_ratio, cds_length, cyclomatic, gamma, meshedness) reads radius as a topological hop count unless distance names an edge attribute, in which case radius is measured along that attribute — confirmed against the installed momepy. distance="mm_len" (the length attribute momepy.gdf_to_nx writes by default) makes radius=5/radius=400 mean 5 m / 400 m of street, matching the only reading of the paper's own numbers that makes sense (5 or 400 topological hops would be a nonsensical pairing on any real network).

Two things about momepy's network functions this module works around, not just uses:

  • nx.get_node_attributes keys its result by momepy's own node identity, an (x, y) coordinate tuple. Handed straight to pd.Series(...), pandas reads a dict of tuple keys as a MultiIndex rather than as one column of coordinate pairs, which breaks a later .reindex() against nodeID-keyed data with an opaque error. Every extraction here goes through nodeID instead, via a coordinate-to-nodeID map built once from momepy.nx_to_gdf's stable node geometry.
  • Several functions have no name= parameter and always write the same attribute (node_density in particular: node_density/node_density_weighted, fixed). Calling one twice at two radii without reading the first result out in between would silently overwrite it — so this module processes one radius fully, in a fixed function order, reading every attribute out before moving to the next.

street_metrics 

street_metrics(etc: GeoDataFrame, streets: GeoDataFrame, buildings_for_profile: GeoDataFrame, *, profile_distance_m: float, profile_tick_length_m: float, node_radii_m: tuple[float, ...] = NODE_RADII_M) -> DataFrame

The 22 Street Descriptors & Connectivity columns, indexed like etc (unit_id).

streets is reset to a clean integer index first: momepy.gdf_to_nx/get_nearest_street both key off it, and a stray non-default index from upstream cleaning would silently change what "street index 8" means between the two calls. buildings_for_profile need not share etc's index — momepy.street_profile only needs building geometry, not identity — so the unfiltered building layer works and there is no reason to pass the ETC-reindexed one.

Source code in src/lczkit/morphometrics/streets.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
def street_metrics(
    etc: gpd.GeoDataFrame,
    streets: gpd.GeoDataFrame,
    buildings_for_profile: gpd.GeoDataFrame,
    *,
    profile_distance_m: float,
    profile_tick_length_m: float,
    node_radii_m: tuple[float, ...] = NODE_RADII_M,
) -> pd.DataFrame:
    """The 22 Street Descriptors & Connectivity columns, indexed like `etc` (`unit_id`).

    `streets` is reset to a clean integer index first: `momepy.gdf_to_nx`/`get_nearest_street`
    both key off it, and a stray non-default index from upstream cleaning would silently change
    what "street index 8" means between the two calls. `buildings_for_profile` need not share
    `etc`'s index — `momepy.street_profile` only needs building geometry, not identity — so the
    unfiltered building layer works and there is no reason to pass the ETC-reindexed one.
    """
    streets = streets.reset_index(drop=True)

    street_length = streets.geometry.length
    street_linearity = momepy.linearity(streets)
    profile = momepy.street_profile(
        streets,
        buildings_for_profile,
        distance=profile_distance_m,
        tick_length=profile_tick_length_m,
    )

    graph = momepy.gdf_to_nx(streets, length="mm_len")
    # Node identity and geometry never change after this point — only attributes are added — so
    # this map is built once and reused for every later extraction.
    coord_to_id = {
        (row.geometry.x, row.geometry.y): row.nodeID
        for row in momepy.nx_to_gdf(graph, points=True, lines=False).itertuples()
    }

    graph = momepy.node_degree(graph)
    graph = momepy.mean_node_dist(graph)
    graph = momepy.clustering(graph)
    node_degree = _by_node_id(graph, "degree", coord_to_id)
    mean_node_distance = _by_node_id(graph, "meanlen", coord_to_id)
    node_clustering = _by_node_id(graph, "cluster", coord_to_id)

    per_radius: dict[float, dict[str, pd.Series]] = {}
    for radius in node_radii_m:
        graph = momepy.mean_node_degree(graph, radius=radius, distance="mm_len", name="_mean_nd")
        graph = momepy.node_density(graph, radius=radius, distance="mm_len")
        density = _by_node_id(graph, "node_density", coord_to_id)
        graph = momepy.edge_node_ratio(graph, radius=radius, distance="mm_len", name="_enr")
        graph = momepy.cds_length(graph, radius=radius, distance="mm_len", name="_cds")
        graph = momepy.cyclomatic(graph, radius=radius, distance="mm_len", name="_cyc")
        graph = momepy.gamma(graph, radius=radius, distance="mm_len", name="_gamma")
        graph = momepy.meshedness(graph, radius=radius, distance="mm_len", name="_mesh")
        per_radius[radius] = {
            "mean_node_degree": _by_node_id(graph, "_mean_nd", coord_to_id),
            "node_density": density,
            "edge_node_ratio": _by_node_id(graph, "_enr", coord_to_id),
            "cds_length": _by_node_id(graph, "_cds", coord_to_id),
            "cyclomatic": _by_node_id(graph, "_cyc", coord_to_id),
            "gamma_index": _by_node_id(graph, "_gamma", coord_to_id),
            "meshedness": _by_node_id(graph, "_mesh", coord_to_id),
        }

    nodes, edges = momepy.nx_to_gdf(graph, points=True, lines=True)

    nearest_street = momepy.get_nearest_street(etc, streets)
    nearest_node = momepy.get_nearest_node(etc, nodes, edges, nearest_street)

    columns: dict[str, pd.Series] = {
        "street_length": _nearest_join(street_length, nearest_street),
        "street_linearity": _nearest_join(street_linearity, nearest_street),
        "street_width": _nearest_join(profile["width"], nearest_street),
        "street_width_deviation": _nearest_join(profile["width_deviation"], nearest_street),
        "street_openness": _nearest_join(profile["openness"], nearest_street),
        "node_degree": _nearest_join(node_degree, nearest_node),
        "mean_node_distance": _nearest_join(mean_node_distance, nearest_node),
        "node_clustering": _nearest_join(node_clustering, nearest_node),
    }
    for radius, values in per_radius.items():
        suffix = f"r{int(radius)}m"
        for label, series in values.items():
            columns[f"{label}_{suffix}"] = _nearest_join(series, nearest_node)

    return pd.DataFrame(columns, index=etc.index)

Contextual expansion

Opt-in: the 25th/50th/75th percentile of each primary attribute across neighbouring ETCs, via momepy.percentile. Off by default (MorphometricsConfig.contextual). Unlike the paper, primary attributes are kept once the expansion is computed rather than dropped.

lczkit.morphometrics.contextual

The opt-in contextual expansion — Majer & Fleischmann (2026) §3.1.

Per primary attribute, the 25th/50th/75th percentile across neighbouring ETCs within config.contextual_steps topological steps, via momepy.percentile (linearly weighted, "hazen" interpolation — momepy's own default and the paper's own choice). Unlike the paper, which drops the primary attributes once the expansion is computed, lczkit keeps both: MorphometricsConfig docstring records this as a deliberate divergence, not an oversight.

contextual_expand 

contextual_expand(primary: DataFrame, graph: Graph, *, quantiles: Iterable[int] = (25, 50, 75)) -> DataFrame

{column}_p{q} for every column of primary and every q in quantiles.

graph should be the ETC contiguity graph expanded to the configured number of topological steps (lczkit.morphometrics.graphs.etc_higher_order) — not self-weighted: the paper's own wording is percentiles across neighbouring cells, so the focal cell's own value must not be counted as one of its own neighbours.

Source code in src/lczkit/morphometrics/contextual.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def contextual_expand(
    primary: pd.DataFrame, graph: Graph, *, quantiles: Iterable[int] = (25, 50, 75)
) -> pd.DataFrame:
    """`{column}_p{q}` for every column of `primary` and every `q` in `quantiles`.

    `graph` should be the ETC contiguity graph expanded to the configured number of topological
    steps (`lczkit.morphometrics.graphs.etc_higher_order`) — *not* self-weighted: the paper's own
    wording is percentiles across *neighbouring* cells, so the focal cell's own value must not be
    counted as one of its own neighbours.
    """
    quantiles = list(quantiles)
    columns: dict[str, pd.Series] = {}
    for column in primary.columns:
        percentiles = momepy.percentile(primary[column], graph, q=quantiles)
        for q in quantiles:
            columns[f"{column}_p{q}"] = percentiles[q]
    return pd.DataFrame(columns, index=primary.index)

Registry

lczkit.morphometrics.registry

What every morphometric attribute column means, in what unit, from which source.

Reuses lczkit.ucp.registry.ParameterSpec rather than duplicating it — a column's shape (name, label, unit, description, reference) doesn't change because the module computing it does. Every one of the 107 primary attributes here has a reference of MOMEPY (the formula and, where the paper's own choice of which momepy call and which scale to use is what's being credited, MAJER_2026) or COMPUTED_HERE for the two lczkit-defined quantities (area_*, coverage_area_ratio_etc) that have no separate published definition beyond arithmetic.

Generated against docs/references/tables/majer_2026_morphometrics_menu.md — every name below must appear as a row name in that table, and tests/test_morphometrics_registry_matches_menu.py asserts it does, cell for cell, so the two cannot silently drift apart.

MAJER_2026 module-attribute 

MAJER_2026 = 'arXiv:2603.00132'

Majer & Fleischmann (2026), "Predicting Local Climate Zones using Urban Morphometrics and Satellite Imagery", arXiv. Credited where the paper's own choice of scale or neighbourhood (rather than momepy's formula itself) is what a column represents — every weighted variant, and the two lczkit-computed granularity/coverage quantities the paper specifies without naming a momepy function.

PARAMETER_COLUMNS module-attribute 

PARAMETER_COLUMNS: tuple[str, ...] = tuple(p.name for p in PARAMETERS)

Column names compute_morphometrics() must emit, in registry order.

CONTEXTUAL_QUANTILES module-attribute 

CONTEXTUAL_QUANTILES: tuple[int, ...] = (25, 50, 75)

The paper's own percentiles for the contextual expansion.

spec 

spec(name: str) -> ParameterSpec

The ParameterSpec for name, raising KeyError if it is not a registered column.

Source code in src/lczkit/morphometrics/registry.py 
1022
1023
1024
def spec(name: str) -> ParameterSpec:
    """The `ParameterSpec` for `name`, raising `KeyError` if it is not a registered column."""
    return _BY_NAME[name]

contextual_specs 

contextual_specs(quantiles: Iterable[int] = CONTEXTUAL_QUANTILES) -> tuple[ParameterSpec, ...]

ParameterSpecs for the opt-in contextual expansion: {name}_p{q} per primary attribute.

Generated from PARAMETERS the same way lczkit.ucp.registry.semantic_specs() generates its config-driven columns — a description template referencing the primary spec's own description, rather than 321 hand-written entries that would drift from the 107 they mirror.

Source code in src/lczkit/morphometrics/registry.py 
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
def contextual_specs(
    quantiles: Iterable[int] = CONTEXTUAL_QUANTILES,
) -> tuple[ParameterSpec, ...]:
    """`ParameterSpec`s for the opt-in contextual expansion: `{name}_p{q}` per primary attribute.

    Generated from `PARAMETERS` the same way `lczkit.ucp.registry.semantic_specs()` generates its
    config-driven columns — a description template referencing the primary spec's own
    description, rather than 321 hand-written entries that would drift from the 107 they mirror.
    """
    result: list[ParameterSpec] = []
    for primary in PARAMETERS:
        for q in quantiles:
            result.append(
                ParameterSpec(
                    name=f"{primary.name}_p{q}",
                    label=f"{primary.label} (p{q}, neighbouring ETCs)",
                    unit=primary.unit,
                    description=(
                        f"The {q}th percentile of `{primary.name}` across neighbouring ETCs "
                        '(`momepy.percentile`, linearly weighted, "hazen" interpolation). '
                        f"{primary.description}"
                    ),
                    reference=MAJER_2026,
                )
            )
    return tuple(result)

Rasterizing

Not part of the paper. Builds a fine grid over the ETC layer's own bounds at the requested resolution, reusing lczkit.units.overlay rather than a new vector-to-raster library — the overlay-and-measure primitive already in the package is exactly what "which ETCs does this pixel cover, and how much of each" is.

Three output formats. "gtiff" (default) and "cog" both go through rasterio's own GDAL drivers. "zarr" goes through xarray/rioxarray rather than GDAL's own Zarr driver — GDAL's driver writes the array correctly but stores the CRS only in a pam.aux.xml sidecar a plain zarr/xarray reader never sees, which was measured directly rather than assumed and defeats the interoperability a cloud-native format exists for. xr.open_zarr(path, decode_coords="all"), or any other CF-aware reader, round-trips the CRS this package writes.

Optional geographic tiling. MorphometricsConfig.raster_tile_deg splits the output into a geotessera-style tile grid (github.com/ucam-eo/geotessera) — tiles named grid_{lon:.2f}_{lat:.2f} after their centre — instead of one file for the whole extent. The expensive step (the area-weighted overlay) still runs exactly once; tiling only changes how the resulting array is windowed and written. Adjacent tiles' shared edges are reprojected and rounded to a pixel index exactly once, so they cannot leave a gap or overlap between two tiles — an earlier version rounded each tile's box independently and did exactly that, caught by a test that stitches tiled output back together and compares it to the untiled raster.

lczkit.morphometrics.raster

Rasterizing morphometric attributes to a multiband raster at a user-defined resolution.

Not part of the paper — Majer & Fleischmann (2026) never rasterize their 2D attributes, only the 20-attribute subset feeding their fusion schemes (out of scope here, see the module docstring of lczkit.morphometrics). This is a new lczkit capability: one raster, one band per attribute, each pixel the area-weighted mean of whichever ETCs overlap it.

Reuses lczkit.units.overlay.unit_pieces rather than a new vector-to-raster library — the overlay-and-measure primitive already in the package is exactly what "which ETCs does this pixel cover, and how much of each" is. A local, purpose-built grid builder is used rather than stretching GridUnits to a second contract: GridUnits takes a lon/lat bbox and estimates its own CRS, which is the wrong shape for a grid built directly over an already-projected ETC layer's own bounds at an arbitrary resolution.

Three output formats. "gtiff" (default, unchanged from the original single-file behaviour) and "cog" both go through rasterio's own GDAL drivers — verified to need no new dependency. "zarr" goes through xarray/rioxarray rather than GDAL's Zarr driver: measured directly before choosing this route (see lczkit.config.RasterFormat) — GDAL's own driver writes the array correctly but stores the CRS only in a GDAL-specific sidecar invisible to a plain zarr/xarray reader, which defeats the interoperability a cloud-native format is for.

Optional tiling. tile_deg splits the raster into a geotessera-style geographic tile grid (github.com/ucam-eo/geotessera) rather than one file covering the whole extent — see lczkit.config.MorphometricsConfig.raster_tile_deg. The expensive step (the area-weighted overlay) runs exactly once regardless of tiling; tiling only changes how the resulting array is windowed and written.

RasterExportReport dataclass 

RasterExportReport(resolution_m: float, n_rows: int, n_cols: int, band_names: tuple[str, ...], format: RasterFormat = 'gtiff', tile_deg: float | None = None, tiles: tuple[str, ...] = ())

What rasterize_attributes wrote, for the run manifest.

tiles class-attribute instance-attribute 

tiles: tuple[str, ...] = ()

Filenames written under out_path when tile_deg is set; empty for a single-file write, where out_path itself (known to the caller) is the whole answer.

as_manifest 

as_manifest() -> dict[str, Any]

The manifest.morphometrics_raster entry describing this raster.

Source code in src/lczkit/morphometrics/raster.py 
88
89
90
91
92
93
94
95
96
97
98
def as_manifest(self) -> dict[str, Any]:
    """The `manifest.morphometrics_raster` entry describing this raster."""
    return {
        "resolution_m": self.resolution_m,
        "n_rows": self.n_rows,
        "n_cols": self.n_cols,
        "band_names": list(self.band_names),
        "format": self.format,
        "tile_deg": self.tile_deg,
        "tiles": list(self.tiles),
    }

raster_filename 

raster_filename(format: RasterFormat) -> str

The single-file output name for format — a .zarr store is a directory, not a .tif.

Source code in src/lczkit/morphometrics/raster.py 
65
66
67
def raster_filename(format: RasterFormat) -> str:
    """The single-file output name for `format` — a `.zarr` store is a directory, not a `.tif`."""
    return RASTER_ZARR_FILE if format == "zarr" else RASTER_FILE

rasterize_attributes 

rasterize_attributes(gdf: GeoDataFrame, resolution_m: float, out_path: Path, *, columns: list[str] | None = None, max_cells: int = 50000000, format: RasterFormat = 'gtiff', tile_deg: float | None = None, max_bytes: int = 4000000000, origin: tuple[float, float] | None = None) -> RasterExportReport

Write out_path as a multiband raster, one band per attribute of gdf, area-weighted.

columns defaults to every non-geometry column of gdf. Raises ValueError before building the grid if it would exceed max_cells (pixel count) or max_bytes (the in-memory band stack — every format now goes through one array built once, see _compute_band_stack) — refusing is cheap; an unbounded allocation is not.

With tile_deg set, out_path is a directory and the output is one file per geographic tile (lczkit.morphometrics.raster._tile_grid) instead of one file for the whole extent. Without it, out_path is the single output file (or, for format="zarr", the single output store directory) — the original, still-default behaviour.

origin anchors the pixel lattice to an external grid rather than to gdf's own bounds — see _pixel_grid. Pass the top-left corner of the grid the output has to register against; the extent is grown outward to that lattice, never cropped to it.

Source code in src/lczkit/morphometrics/raster.py 
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
def rasterize_attributes(
    gdf: gpd.GeoDataFrame,
    resolution_m: float,
    out_path: Path,
    *,
    columns: list[str] | None = None,
    max_cells: int = 50_000_000,
    format: RasterFormat = "gtiff",
    tile_deg: float | None = None,
    max_bytes: int = 4_000_000_000,
    origin: tuple[float, float] | None = None,
) -> RasterExportReport:
    """Write `out_path` as a multiband raster, one band per attribute of `gdf`, area-weighted.

    `columns` defaults to every non-geometry column of `gdf`. Raises `ValueError` before building
    the grid if it would exceed `max_cells` (pixel count) or `max_bytes` (the in-memory band
    stack — every format now goes through one array built once, see `_compute_band_stack`) —
    refusing is cheap; an unbounded allocation is not.

    With `tile_deg` set, `out_path` is a directory and the output is one file per geographic tile
    (`lczkit.morphometrics.raster._tile_grid`) instead of one file for the whole extent. Without
    it, `out_path` is the single output file (or, for `format="zarr"`, the single output store
    directory) — the original, still-default behaviour.

    `origin` anchors the pixel lattice to an external grid rather than to `gdf`'s own bounds — see
    `_pixel_grid`. Pass the top-left corner of the grid the output has to register against; the
    extent is grown outward to that lattice, never cropped to it.
    """
    assert_projected_crs(gdf, "gdf")
    if resolution_m <= 0:
        raise ValueError(f"resolution_m must be positive, got {resolution_m}")
    if tile_deg is not None and tile_deg <= 0:
        raise ValueError(f"tile_deg must be positive, got {tile_deg}")
    assert gdf.crs is not None  # narrows for mypy; assert_projected_crs already guarantees this
    attribute_columns = (
        columns if columns is not None else [c for c in gdf.columns if c != "geometry"]
    )

    raw_bounds = gdf.total_bounds
    bounds = (
        _snap((raw_bounds[0], raw_bounds[1], raw_bounds[2], raw_bounds[3]), resolution_m, origin)
        if origin is not None
        else (raw_bounds[0], raw_bounds[1], raw_bounds[2], raw_bounds[3])
    )
    grid, n_rows, n_cols = _pixel_grid(bounds, resolution_m)
    if n_rows * n_cols > max_cells:
        raise ValueError(
            f"a {resolution_m} m grid over this extent would be {n_rows}x{n_cols} = "
            f"{n_rows * n_cols} cells, over the configured ceiling of {max_cells} "
            "(MorphometricsConfig.max_raster_cells)"
        )
    projected_bytes = len(attribute_columns) * n_rows * n_cols * 4
    if projected_bytes > max_bytes:
        raise ValueError(
            f"a {resolution_m} m grid over this extent with {len(attribute_columns)} attributes "
            f"would need ~{projected_bytes / 1e9:.2f} GB in memory (bands x rows x cols x 4 "
            f"bytes), over the configured ceiling of {max_bytes / 1e9:.2f} GB "
            "(MorphometricsConfig.max_raster_bytes) — rasterize fewer columns, a coarser "
            "resolution, or raise the ceiling"
        )
    grid = grid.set_geometry(grid.geometry, crs=gdf.crs)

    pieces = unit_pieces(grid, gdf, columns=attribute_columns)
    stack = _compute_band_stack(pieces, attribute_columns, grid, n_rows, n_cols)

    minx, _, _, maxy = bounds
    transform = from_origin(minx, maxy, resolution_m, resolution_m)
    crs = CRS.from_user_input(gdf.crs)

    if tile_deg is None:
        out_path.parent.mkdir(parents=True, exist_ok=True)
        _write_raster(stack, transform, crs, attribute_columns, out_path, format=format)
        return RasterExportReport(
            resolution_m=resolution_m,
            n_rows=n_rows,
            n_cols=n_cols,
            band_names=tuple(attribute_columns),
            format=format,
            tile_deg=None,
            tiles=(),
        )

    out_path.mkdir(parents=True, exist_ok=True)
    ext = _tile_ext(format)
    raster_bounds = (minx, maxy - n_rows * resolution_m, minx + n_cols * resolution_m, maxy)
    written: list[str] = []
    for name, window in _tile_grid(raster_bounds, crs, tile_deg, transform, n_rows, n_cols):
        row_off, col_off = int(window.row_off), int(window.col_off)
        tile_array = stack[
            :, row_off : row_off + int(window.height), col_off : col_off + int(window.width)
        ]
        if not np.isfinite(tile_array).any():
            continue
        tile_transform = window_transform(window, transform)
        filename = f"{name}{ext}"
        _write_raster(
            tile_array, tile_transform, crs, attribute_columns, out_path / filename, format=format
        )
        written.append(filename)

    return RasterExportReport(
        resolution_m=resolution_m,
        n_rows=n_rows,
        n_cols=n_cols,
        band_names=tuple(attribute_columns),
        format=format,
        tile_deg=tile_deg,
        tiles=tuple(sorted(written)),
    )

refresh_raster 

refresh_raster(run_dir: Path, resolution_m: float, *, columns: list[str] | None = None, max_cells: int = 50000000, format: RasterFormat = 'gtiff', tile_deg: float | None = None, max_bytes: int = 4000000000) -> RasterExportReport

Re-derive the morphometrics raster(s) from an already-written run's morphometrics.parquet.

The one function both run_pipeline (when --morphometrics-resolution is given) and lczkit morphometrics raster call — a run gets the same raster whether it was produced at run time or requested afterwards at a different resolution or format, because both paths go through this. Patches the manifest as JSON, the same way lczkit.output.gis's backfill does: this can run against a manifest written by an older version of RunManifest, and round-tripping it through today's model would rewrite fields the run never had.

A previous single-file or tiled output at a different format is never deleted — only appended to manifest.outputs — matching the project's standing convention of adding rather than overwriting a run's on-disk record (e.g. lczkit export's units.gpkg).

Source code in src/lczkit/morphometrics/raster.py 
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
def refresh_raster(
    run_dir: Path,
    resolution_m: float,
    *,
    columns: list[str] | None = None,
    max_cells: int = 50_000_000,
    format: RasterFormat = "gtiff",
    tile_deg: float | None = None,
    max_bytes: int = 4_000_000_000,
) -> RasterExportReport:
    """Re-derive the morphometrics raster(s) from an already-written run's `morphometrics.parquet`.

    The one function both `run_pipeline` (when `--morphometrics-resolution` is given) and
    `lczkit morphometrics raster` call — a run gets the same raster whether it was produced at
    run time or requested afterwards at a different resolution or format, because both paths go
    through this. Patches the manifest as JSON, the same way `lczkit.output.gis`'s backfill does:
    this can run against a manifest written by an older version of `RunManifest`, and
    round-tripping it through today's model would rewrite fields the run never had.

    A previous single-file or tiled output at a different format is never deleted — only appended
    to `manifest.outputs` — matching the project's standing convention of adding rather than
    overwriting a run's on-disk record (e.g. `lczkit export`'s `units.gpkg`).
    """
    morphometrics_path = run_dir / MORPHOMETRICS_FILE
    if not morphometrics_path.exists():
        raise FileNotFoundError(
            f"no {MORPHOMETRICS_FILE} in {run_dir} — this run has no morphometrics to rasterize"
        )
    gdf = gpd.read_parquet(morphometrics_path)
    target = (
        run_dir / RASTER_TILE_DIR if tile_deg is not None else run_dir / raster_filename(format)
    )
    report = rasterize_attributes(
        gdf,
        resolution_m,
        target,
        columns=columns,
        max_cells=max_cells,
        format=format,
        tile_deg=tile_deg,
        max_bytes=max_bytes,
    )

    manifest_path = run_dir / MANIFEST_FILE
    if manifest_path.exists():
        manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
        manifest["morphometrics_raster"] = report.as_manifest()
        outputs = list(manifest.get("outputs", []))
        output_entry = RASTER_TILE_DIR if tile_deg is not None else raster_filename(format)
        if output_entry not in outputs:
            outputs.append(output_entry)
        manifest["outputs"] = outputs
        manifest_path.write_text(json.dumps(manifest, indent=2) + "\n", encoding="utf-8")

    return report

Report

lczkit.morphometrics.report

What a morphometrics run produced, for the run manifest.

MorphometricsReport dataclass 

MorphometricsReport(tessellation: TessellationReport, n_primary_attributes: int, contextual_enabled: bool, n_contextual_attributes: int)

Counts and settings a reader of the manifest needs to interpret morphometrics.parquet.

n_contextual_attributes instance-attribute 

n_contextual_attributes: int

0 when contextual_enabled is False — the two fields are kept separate rather than one being inferred from the other, so a manifest reader never has to know the rule.