Validation

Agreement against a reference LCZ map, in the style of the lczexplore package.

Per-class agreement and a confusion matrix, never a single accuracy number; plus the breakdowns that turn the height caveat into a measured quantity - agreement by height_completeness band, and the two confusion axes reported apart: the height axis (1<->2<->3, 4<->5<->6) and the compactness axis (1<->4, 2<->5, 3<->6).

ranges adds the complementary instrument: agreement says whether a label was right, and parameter_ranges says whether the parameter behind it could ever have reached the published prototype it should have matched. uncertainty adds the third: how much of a difference between two agreement figures is signal, given that the units are a correlated sheet rather than independent draws.

Three references, not one, and they are three different kinds of object. reference_lcz reads the Demuzere global map, which is a model output carrying its own error; labelled_lcz reads hand-labelled So2Sat LCZ42 / DFC2017 patches, which are ground truth on a uniform grid; wudapt_lcz reads WUDAPT LCZ Generator training areas, which are hand-drawn too but irregular, overlapping, spread over four decades, and available for every city rather than fifty-one.

Where So2Sat exists it is primary and the agreement between it and lcz_v3 is the ceiling on any score against the map, and that distinction is easy to lose. WUDAPT does not give a second such ceiling: the LCZ Generator's training areas are the training data behind the Demuzere map, so the two are not independent and their agreement is inflated by construction.

All three return the same three columns, so agreement() consumes any of them without knowing which — which is exactly why every report records reference_file. A table that does not name the file it used is indistinguishable from one that used a different reference, and this project has found that mistake twice.

Agreement between a run's labels and a reference set, reported in the style of the lczexplore package — per-class figures and a sparse confusion matrix, never a single accuracy number.

Terms used throughout this page are defined in the glossary: ceiling, overall accuracy, built-class agreement and the two confusion axes.

Validate against labelled ground truth, not against another model

lcz_v3.tif — the Demuzere global LCZ map, derived from satellite imagery — is an estimate carrying its own error. Measuring against it compares two models and reports the disagreement as lczkit's error. Where hand-drawn polygons exist — So2Sat LCZ42, a labelled benchmark set over 51 cities, and WUDAPT, the community-contributed World Urban Database and Access Portal Tools — those are the reference and lcz_v3 is a secondary comparator. reference_file is recorded on every report: "the reference" names a role, and the three that can fill it disagree by up to 18 points, so a figure that does not say which file produced it cannot be read.

Three things any per-city figure needs stated beside it:

• Its ceiling. How well lcz_v3 itself agrees with the hand-drawn polygons on the same cells. That bounds what any map can score against lcz_v3, and it varies enormously: 22.8% in Mumbai, 83.2% in Rio.
• Built-class agreement, separately, with the natural-class share alongside. An overall figure dominated by trivially-classified water says nothing about the classifier.
• The label reproducibility. Two independent expert label sets over the same ground agree at a median 79.7% across twenty-eight cities, ranging 26.3% (Cairo — below what a constant predictor scores there, 52.1%) to 97.7% (Istanbul). Where two references disagree, no classifier can agree with both, so this is a floor under every residual reported here.

Never report '% of ceiling'

Vancouver scores 41.8% against a 36.7% ceiling — 114%. The comparator is another estimator, not an upper bound. Report raw agreement and ceiling side by side, or their difference.

lczkit.validation.agreement

Agreement against a reference LCZ map, reported the way the LCZ literature reports it.

Comparability matters more than a headline figure: per-class agreement and a confusion matrix, in the style of the lczexplore package, not a single accuracy number. A single number hides the thing a reader of an LCZ map actually needs to know, which is which classes are being confused for which - Bernard et al. (2024) Sect. 3.2 read their own results exactly this way, noting that their worst agreements sit on classes covering a negligible share of the area.

Four breakdowns beyond the standard ones are required here:

• built-class agreement, separately from overall, with the natural-class share stated beside it. An overall figure carried by water is not a statement about the classifier;
• agreement by height_completeness band, so a city where tier-1 heights are near-absent can be judged rather than disclaimed;
• the height axis, 1<->2<->3 and 4<->5<->6, holding compactness fixed. This is where error concentrates when heights come from an areal product, since such a product cannot resolve the <10 m / 10-25 m / >25 m bands within a heterogeneous unit. It is the axis that pairs with the stratification above; and
• the compactness axis, 1<->4, 2<->5, 3<->6, holding height fixed. This is the diagnostic for building footprint coverage and for whether the spatial unit is the right size to hold an LCZ patch at all.

The two are different instruments and are reported separately. An earlier revision of the spec named the second set the height axis, which inverted what a reader would conclude from it: a disagreement between LCZ 2 and LCZ 5 says nothing about heights, both being midrise.

Agreement is area-weighted throughout. On a regular grid that is the same as counting units; on enclosures it is not, and weighting by count there would let a thousand courtyards outvote a district. The axes are the exception, and deliberately carry both weightings: every axis figure published between Phases 6.7 and 11 is count-based, so share_of_disagreement keeps that definition and share_of_disagreement_area is reported beside it. Redefining the field in place would have silently moved every stored arm-B number while leaving arm A untouched, which is the hardest kind of change to notice.

Raw axis shares are not comparable across cities, which is easy to get wrong. Their denominator is all disagreement, so a city whose reference carries a large natural share dilutes both axes, while an all-built fixture dilutes neither; and which pairs can fire at all depends on which classes the reference happens to contain. axis_summary() below reports the raw share alongside two corrections for exactly this - see its docstring.

AXIS_ELIGIBLE_CODES module-attribute

AXIS_ELIGIBLE_CODES: frozenset[int] = frozenset(code for pair in HEIGHT_AXIS_PAIRS + COMPACTNESS_AXIS_PAIRS for code in pair)

LCZ 1-6. The only reference classes that can place a unit on either axis.

Both axis definitions are built from the compact (1-3) and open (4-6) types, so a unit whose reference is LCZ 7-10 or a natural class contributes to the denominator of a raw axis share while being unable to contribute to its numerator.

ClassAgreement

Bases: BaseModel

Agreement for one reference class.

n_reference instance-attribute

n_reference: int

Units the reference map assigns to this class.

agreement instance-attribute

agreement: float

area_agree_m2 / area_reference_m2, or 0.0 where the class is absent from the reference. This is producer's accuracy (recall): of the ground the reference calls this class, how much did lczkit agree on.

n_predicted class-attribute instance-attribute

n_predicted: int = 0

Units lczkit assigns to this class. The denominator of user_accuracy.

user_accuracy class-attribute instance-attribute

user_accuracy: float = 0.0

area_agree_m2 / area_predicted_m2, or 0.0 where the class is never predicted.

The complement of agreement, and the reason both are needed: producer's accuracy alone cannot see over-prediction. A classifier that labels every cell LCZ 5 scores 100% producer's accuracy on LCZ 5 while being useless, and nothing in this report would have said so. Demuzere et al. (2021) report the pair through F1 for exactly this reason.

f1 class-attribute instance-attribute

f1: float = 0.0

Harmonic mean of user_accuracy and agreement, 0.0 where either is zero.

The class-wise metric Demuzere et al. (2021), Sect. 2.4 use, following Verdonck et al. (2017): "the F1 metric, which is a harmonic mean of the user's and producer's accuracy". Reported so lczkit's per-class figures are directly comparable to published LCZ maps.

ConfusionCell

Bases: BaseModel

One (reference, predicted) pair with a non-zero count.

Stratum

Bases: BaseModel

Agreement within one band of a stratifying variable.

AxisConfusion

Bases: BaseModel

One pair of classes differing along a single axis, counted in both directions.

n_a_as_b instance-attribute

n_a_as_b: int

Units the reference calls a and the run called b.

share_of_disagreement instance-attribute

share_of_disagreement: float

n_total over all disagreeing units. If one axis is the limiting factor, its pairs should hold a disproportionate share of the errors - which is the whole reason for reporting the two axes apart rather than pooled.

share_of_disagreement_area class-attribute instance-attribute

share_of_disagreement_area: float = 0.0

The same share, area-weighted. Identical to share_of_disagreement on a regular grid, where every unit is one cell; different on enclosures, where units vary in size by orders of magnitude.

AxisSummary

Bases: BaseModel

One axis pooled over its pairs, with the denominators that make cities comparable.

Three shares of the same numerator, in increasing order of correction:

• share_of_disagreement is the raw figure Phases 6.7-11 published. Its denominator is all disagreement, so it is deflated by whatever built<->natural confusion a city happens to carry and cannot be compared across cities with different class composition;
• share_of_axis_eligible restricts the denominator to units whose reference is LCZ 1-6, the only classes that can land on either axis. Transparent, and enough to remove the natural-share dilution;
• lift divides by expected_share, what the axis would hold if misclassification ignored the axes entirely. This is the one that removes class composition: a reference carrying only LCZ 2 and 5 affords the compactness pair both its directions while giving each height pair only one, so the two axes are not on equal footing before it is applied.

lift is the figure to compare across cities. 1.0 means the axis holds exactly the share its reference affords it; above 1.0 means error concentrates there beyond what composition explains.

expected_share instance-attribute

expected_share: float

Share the axis would hold under a null that keeps each unit's reference class and the observed distribution of wrong labels, but breaks the association between them.

lift instance-attribute

lift: float

share_of_disagreement / expected_share, or 0.0 where the null affords the axis nothing.

AgreementReport

Bases: BaseModel

The full comparison against one reference map.

n_units instance-attribute

n_units: int

Units in the run.

n_compared instance-attribute

n_compared: int

Units that survived the coverage filter and carry both a prediction and a reference.

overall_agreement instance-attribute

overall_agreement: float

Area-weighted share of compared units whose label matches the reference.

built_agreement class-attribute instance-attribute

built_agreement: float = 0.0

Agreement over the units the reference calls built (LCZ 1-10).

Reported separately from overall_agreement always, because an overall figure can be dominated by trivially-classified natural cover and then says nothing about the classifier: Rotterdam's headline 42.5% was 266 water cells agreeing at 95.9% while LCZ 8 sat at 0.0% over 224. The denominator is the reference's family, not the run's - grouping by what the run predicted would let a classifier improve its built score by predicting water.

natural_agreement class-attribute instance-attribute

natural_agreement: float = 0.0

Agreement over the units the reference calls natural (LCZ A-G, codes 11-17).

natural_share class-attribute instance-attribute

natural_share: float = 0.0

Area share of the compared units the reference calls natural. Stated alongside the headline so a figure carried by water is recognisable as one.

weighted_agreement class-attribute instance-attribute

weighted_agreement: float = 0.0

OA_w: overall accuracy with partial credit for a near-miss, per Bechtel et al. (2020).

Every other agreement figure here treats calling a compact midrise cell "open midrise" as exactly as wrong as calling it "water". For a scheme whose classes lie on a near-continuum of built form that is plainly false, and it is the reason this module reports the two confusion axes apart in the first place. OA_w weights the confusion matrix by class similarity, so an adjacent-class error scores most of a correct one and a cross-family error scores little.

A generalisation of overall_agreement rather than a rival to it: plain OA is the same sum with a similarity matrix of ones on the diagonal and zeros off it, so the two coincide exactly when the matrix is the identity. Reported beside it, never instead of it — and the LCZ literature reports both, which is what makes lczkit's numbers comparable to a published map's.

Count-weighted, matching the confusion matrix it is computed from. overall_agreement is area-weighted; on a grid the two coincide, on enclosures they do not.

built_natural_agreement class-attribute instance-attribute

built_natural_agreement: float = 0.0

OA_bu: agreement on the built-versus-natural distinction alone, ignoring which built or which natural class.

Demuzere et al. (2021) Sect. 2.4 and Demuzere et al. (2022) Sect. 2.4 both report it beside OA and OA_u, and it separates two failures this report otherwise conflates: a city where lczkit finds the built fabric and misjudges its form, and one where it does not find the built fabric at all. overall_agreement charges both the same. It is also the one accuracy figure here that a height-blind pipeline can still score well on, which is the point when tier-1 coverage is 1%.

height_axis class-attribute instance-attribute

height_axis: list[AxisConfusion] = Field(default_factory=list)

1<->2<->3 and 4<->5<->6: compactness fixed, height band varies.

compactness_axis class-attribute instance-attribute

compactness_axis: list[AxisConfusion] = Field(default_factory=list)

1<->4, 2<->5, 3<->6: height fixed, building surface fraction varies.

compactness_axis_summary class-attribute instance-attribute

compactness_axis_summary: AxisSummary | None = None

The two axes pooled, with the normalised denominators. Compare cities on lift, never on share_of_disagreement, which is not comparable across cities whose references carry different classes.

n_disagree_axis_eligible class-attribute instance-attribute

n_disagree_axis_eligible: int = 0

Disagreeing units whose reference is LCZ 1-6. The denominator of share_of_axis_eligible.

agreement

agreement(predicted: Series, reference: Series, area_m2: Series, *, coverage: Series | None = None, height_completeness: Series | None = None, config: ValidationConfig | None = None, reference_file: str | None = None) -> AgreementReport

Compare predicted against reference, both indexed by unit_id.

Units missing either label, or whose reference coverage falls below the configured minimum, are excluded and counted separately rather than silently dropped - a run comparing a tenth of its units against the reference is a different claim from one comparing all of them.

Source code in src/lczkit/validation/agreement.py

def agreement(
    predicted: pd.Series,
    reference: pd.Series,
    area_m2: pd.Series,
    *,
    coverage: pd.Series | None = None,
    height_completeness: pd.Series | None = None,
    config: ValidationConfig | None = None,
    reference_file: str | None = None,
) -> AgreementReport:
    """Compare `predicted` against `reference`, both indexed by `unit_id`.

    Units missing either label, or whose reference coverage falls below the configured minimum,
    are excluded and counted separately rather than silently dropped - a run comparing a tenth of
    its units against the reference is a different claim from one comparing all of them.
    """
    settings = config or ValidationConfig()
    frame = pd.DataFrame(
        {
            "predicted": pd.to_numeric(predicted, errors="coerce"),
            "reference": pd.to_numeric(reference, errors="coerce"),
            "area": pd.to_numeric(area_m2, errors="coerce"),
        }
    )
    frame["coverage"] = 1.0 if coverage is None else pd.to_numeric(coverage, errors="coerce")

    no_prediction = frame["predicted"].isna()
    no_reference = frame["reference"].isna()
    low_coverage = (
        ~no_prediction & ~no_reference & (frame["coverage"] < settings.min_reference_coverage)
    ).fillna(True)

    compared = frame[~no_prediction & ~no_reference & ~low_coverage].copy()
    compared["agree"] = compared["predicted"] == compared["reference"]
    total_area = float(compared["area"].sum())

    report = AgreementReport(
        reference_citation=settings.reference_citation,
        reference_file=reference_file,
        n_units=int(len(frame)),
        n_compared=int(len(compared)),
        excluded_no_prediction=int(no_prediction.sum()),
        excluded_no_reference=int((~no_prediction & no_reference).sum()),
        excluded_low_coverage=int(low_coverage.sum()),
        min_reference_coverage=settings.min_reference_coverage,
        overall_agreement=_share(compared.loc[compared["agree"], "area"].sum(), total_area),
        area_compared_m2=total_area,
        n_disagree=int((~compared["agree"]).sum()),
    )
    if compared.empty:
        return report

    is_natural = compared["reference"].isin(NATURAL_CODES)
    built, natural = compared[~is_natural], compared[is_natural]
    report.n_built = int(len(built))
    report.n_natural = int(len(natural))
    report.built_agreement = _share(
        float(built.loc[built["agree"], "area"].sum()), float(built["area"].sum())
    )
    report.natural_agreement = _share(
        float(natural.loc[natural["agree"], "area"].sum()), float(natural["area"].sum())
    )
    report.natural_share = _share(float(natural["area"].sum()), total_area)
    # OA_bu: collapse both sides to the built/natural dichotomy and score that alone.
    same_family = is_natural == compared["predicted"].isin(NATURAL_CODES)
    report.built_natural_agreement = _share(
        float(compared.loc[same_family, "area"].sum()), total_area
    )

    report.per_class = _per_class(compared)
    report.confusion = _confusion(compared)
    report.weighted_agreement = weighted_agreement(report.confusion)
    report.height_axis = axis_pairs(report.confusion, HEIGHT_AXIS_PAIRS)
    report.compactness_axis = axis_pairs(report.confusion, COMPACTNESS_AXIS_PAIRS)
    report.height_axis_summary = axis_summary(report.confusion, HEIGHT_AXIS_PAIRS, axis="height")
    report.compactness_axis_summary = axis_summary(
        report.confusion, COMPACTNESS_AXIS_PAIRS, axis="compactness"
    )
    report.n_disagree_axis_eligible = report.height_axis_summary.n_axis_eligible
    if height_completeness is not None:
        report.by_height_completeness = _strata(
            compared,
            pd.to_numeric(height_completeness, errors="coerce").reindex(compared.index),
            settings.height_completeness_deciles,
        )
    return report

weighted_agreement

weighted_agreement(confusion: Sequence[ConfusionCell], weights: Mapping[tuple[int, int], float] | None = None) -> float

OA_w — Bechtel et al. (2020) Eq. 1: sum(w_ij * c_ij) / N.

w is the similarity matrix from docs/references/tables/lcz_class_similarity.md, one on the diagonal. Passing the complement would invert the measure without raising, which is why the default is not a parameter a caller has to get right and why similarity._check() refuses a matrix whose diagonal is not one.

Computed from the stored confusion list, so a manifest written before this existed can be re-scored without re-running the pipeline.

Source code in src/lczkit/validation/agreement.py

def weighted_agreement(
    confusion: Sequence[ConfusionCell],
    weights: Mapping[tuple[int, int], float] | None = None,
) -> float:
    """`OA_w` — Bechtel et al. (2020) Eq. 1: `sum(w_ij * c_ij) / N`.

    `w` is the **similarity** matrix from `docs/references/tables/lcz_class_similarity.md`, one on
    the diagonal. Passing the complement would invert the measure without raising, which is why the
    default is not a parameter a caller has to get right and why `similarity._check()` refuses a
    matrix whose diagonal is not one.

    Computed from the stored confusion list, so a manifest written before this existed can be
    re-scored without re-running the pipeline.
    """
    table = SIMILARITY if weights is None else weights
    total = sum(cell.n for cell in confusion)
    if total == 0:
        return 0.0
    scored = sum(table.get((cell.reference, cell.predicted), 0.0) * cell.n for cell in confusion)
    return float(scored / total)

axis_pairs

axis_pairs(confusion: Sequence[ConfusionCell], pairs: tuple[tuple[int, int], ...]) -> list[AxisConfusion]

One axis' pairs and their share of all disagreement, counted in both directions.

Derived from the confusion matrix rather than from the compared units, so that re-analysing a stored run cannot drift from what the run itself reported: the confusion list is what a run persists, and it holds exactly the counts and areas this needs.

Source code in src/lczkit/validation/agreement.py

def axis_pairs(
    confusion: Sequence[ConfusionCell], pairs: tuple[tuple[int, int], ...]
) -> list[AxisConfusion]:
    """One axis' pairs and their share of all disagreement, counted in both directions.

    Derived from the confusion matrix rather than from the compared units, so that re-analysing a
    stored run cannot drift from what the run itself reported: the confusion list is what a run
    persists, and it holds exactly the counts and areas this needs.
    """
    disagreeing = [cell for cell in confusion if cell.reference != cell.predicted]
    n_disagree = float(sum(cell.n for cell in disagreeing))
    area_disagree = float(sum(cell.area_m2 for cell in disagreeing))
    counts = {(cell.reference, cell.predicted): cell for cell in disagreeing}

    rows: list[AxisConfusion] = []
    for a, b in pairs:
        forward, backward = counts.get((a, b)), counts.get((b, a))
        n_forward = forward.n if forward else 0
        n_backward = backward.n if backward else 0
        area = (forward.area_m2 if forward else 0.0) + (backward.area_m2 if backward else 0.0)
        rows.append(
            AxisConfusion(
                a=a,
                b=b,
                n_a_as_b=n_forward,
                n_b_as_a=n_backward,
                n_total=n_forward + n_backward,
                share_of_disagreement=_share(float(n_forward + n_backward), n_disagree),
                area_m2=area,
                share_of_disagreement_area=_share(area, area_disagree),
            )
        )
    return rows

axis_summary

axis_summary(confusion: Sequence[ConfusionCell], pairs: tuple[tuple[int, int], ...], *, axis: str) -> AxisSummary

Pool one axis over its pairs and normalise it for the reference's class composition.

The null behind expected_share keeps each disagreeing unit's reference class and the observed distribution of wrong labels, but breaks the association between them: a unit whose reference is r takes a wrong label drawn from the run's own error distribution, conditioned on not being r. Formally, for axis X with partner sets P,

E[X] = sum_r (n_r / N) * (sum_{p in P(r)} q_p) / (1 - q_r)

with n_r the disagreeing units the reference calls r, N the disagreeing units in total, and q the distribution of predicted labels across them.

Keeping q from the run rather than assuming it uniform matters: a classifier that over-predicts LCZ 5 raises the chance of landing on 2<->5 for reasons that have nothing to do with compactness being the limiting axis, and the null has to be able to absorb that or lift would just rediscover the prediction histogram.

Source code in src/lczkit/validation/agreement.py

def axis_summary(
    confusion: Sequence[ConfusionCell],
    pairs: tuple[tuple[int, int], ...],
    *,
    axis: str,
) -> AxisSummary:
    """Pool one axis over its pairs and normalise it for the reference's class composition.

    The null behind `expected_share` keeps each disagreeing unit's *reference* class and the
    observed distribution of wrong labels, but breaks the association between them: a unit whose
    reference is `r` takes a wrong label drawn from the run's own error distribution, conditioned on
    not being `r`. Formally, for axis X with partner sets `P`,

        E[X] = sum_r (n_r / N) * (sum_{p in P(r)} q_p) / (1 - q_r)

    with `n_r` the disagreeing units the reference calls `r`, `N` the disagreeing units in total,
    and `q` the distribution of predicted labels across them.

    Keeping `q` from the run rather than assuming it uniform matters: a classifier that
    over-predicts LCZ 5 raises the chance of landing on 2<->5 for reasons that have nothing to do
    with compactness being the limiting axis, and the null has to be able to absorb that or `lift`
    would just rediscover the prediction histogram.
    """
    rows = axis_pairs(confusion, pairs)
    disagreeing = [cell for cell in confusion if cell.reference != cell.predicted]
    n_disagree = float(sum(cell.n for cell in disagreeing))
    n_eligible = sum(cell.n for cell in disagreeing if cell.reference in AXIS_ELIGIBLE_CODES)
    n_total = sum(row.n_total for row in rows)

    partners = _partners(pairs)
    predicted_marginal: dict[int, float] = {}
    reference_counts: dict[int, int] = {}
    for cell in disagreeing:
        predicted_marginal[cell.predicted] = predicted_marginal.get(cell.predicted, 0.0) + cell.n
        reference_counts[cell.reference] = reference_counts.get(cell.reference, 0) + cell.n
    if n_disagree > 0:
        predicted_marginal = {code: n / n_disagree for code, n in predicted_marginal.items()}

    expected = 0.0
    for code, count in reference_counts.items():
        remaining = 1.0 - predicted_marginal.get(code, 0.0)
        if remaining <= 0.0:
            continue
        reachable = sum(predicted_marginal.get(other, 0.0) for other in partners.get(code, ()))
        expected += (count / n_disagree) * (reachable / remaining)

    observed = _share(float(n_total), n_disagree)
    return AxisSummary(
        axis=axis,
        n_total=n_total,
        share_of_disagreement=observed,
        share_of_disagreement_area=sum(row.share_of_disagreement_area for row in rows),
        n_axis_eligible=int(n_eligible),
        share_of_axis_eligible=_share(float(n_total), float(n_eligible)),
        expected_share=expected,
        lift=_share(observed, expected),
    )

Confusion axes

Two different ways a label can be wrong, and they diagnose different things, so they must not be conflated:

• Height axis — confusions between 1↔2↔3 and between 4↔5↔6. Compactness is fixed and the height band varies, so this tracks the quality of the height data. Read it beside height_completeness.
• Compactness axis — confusions between 1↔4, 2↔5 and 3↔6. Height is fixed and building surface fraction varies, so this tracks how completely the footprints were captured and how the units were drawn.

Only pair-normalised lift against a composition-preserving null is reported. The raw share cannot compare the two: the height axis affords six pairs to compactness's three, so a null that never looks at the data awards height 3.9× more error on affordance alone.

Weighted accuracy

Plain overall accuracy (OA) counts a unit right only if its label matches the reference exactly. Weighted overall accuracy (OA_w) gives partial credit according to how similar the two classes are, using the class-similarity matrix of Bechtel, Demuzere & Stewart (2020) — so calling a compact midrise an open midrise scores most of a point, while calling it water scores near zero.

OA_w above OA therefore says the map is landing in neighbouring classes rather than at random. That is a real and useful thing to know, and it is not an accuracy figure. It is reported beside OA and never instead of it.

lczkit.validation.similarity

Bechtel et al. (2020) class-similarity weights, and the weighted accuracy they define.

Why this is not just another accuracy number. Overall accuracy treats every mistake as equally wrong: calling a compact midrise cell "open midrise" costs exactly what calling it "water" costs. For LCZ that is plainly false — the scheme's classes lie on a near-continuum of built form, and the whole point of reporting the two confusion axes separately is that adjacent-class error means something different from cross-family error. OA_w is the LCZ community's answer, and reporting it is what makes this package's per-class figures comparable to published LCZ maps.

Bechtel et al. frame it as a generalisation of overall accuracy rather than a different measure. Plain OA is already the Hadamard product of the confusion matrix with a weight matrix carrying ones on the diagonal and zeros off it — itself a similarity metric, in which every class is fully similar to itself and fully dissimilar to everything else. Relaxing those off-diagonal zeros to real similarities buys partial credit for a near-miss. So OA_w reduces to OA exactly when the matrix is the identity, which is the first thing test_validation_similarity.py asserts.

The weights live in docs/references/tables/lcz_class_similarity.md, not here. A published number is read from a committed transcription rather than reproduced from memory, and this matrix has 289 cells, of which a wrong one would be invisible. The table is parsed at import and the parsed values are the constants; there is no second copy to drift.

That file holds two matrices with identical headers — the similarity one and its complement. OA_w uses the similarity matrix. Substituting the other inverts the measure without raising: a perfect map scores 0.00, a map confusing LCZ 1 with LCZ G scores 1.00, and every cross-city comparison in a paper ranks backwards. Hence the parsing is by section heading and the direction is asserted by a test rather than trusted to a comment.

BECHTEL_2020 module-attribute

BECHTEL_2020 = '10.3390/rs12111769'

Bechtel, Demuzere & Stewart (2020), Remote Sensing 12(11), 1769. Defines OA_w and the class similarity matrix it weights the confusion matrix with.

SIMILARITY_TABLE module-attribute

SIMILARITY_TABLE = Path(__file__).resolve().parents[3] / 'docs' / 'references' / 'tables' / 'lcz_class_similarity.md'

The committed transcription. Repo content rather than data, so it needs no DATA_DIR — the same exception tests/fixtures/ relies on, and for the same reason: a clean checkout must be able to reproduce a classification.

SIMILARITY module-attribute

SIMILARITY: dict[tuple[int, int], float] = _parse(SIMILARITY_TABLE, SIMILARITY_HEADING)

(reference code, predicted code) -> similarity in [0, 1]. One on the diagonal.

References

lczkit.validation.labelled

Labelled LCZ ground truth (So2Sat LCZ42 / DFC2017), reduced to one label per spatial unit.

lcz_v3.tif is an estimate carrying its own error, so measuring against it compares two models and reports the disagreement as lczkit's. Where hand-labelled LCZ polygons exist they are the primary reference and lcz_v3 is a secondary comparator; the agreement between the two, on the same cells, is the ceiling on what lczkit can score against lcz_v3 at all.

Why the label is anchored on the patch centre rather than overlaid areally. So2Sat patches are 320 m squares sampled on a 100 m stride, so they overlap each other about sevenfold - measured on the Berlin fixture bbox: 473 patches, 48.4 km2 of patch area over a 7.1 km2 union, 16,560 overlapping pairs. An areal overlay would therefore count the same ground up to nine times, under labels that need not agree, and the resulting "majority" would be a property of the sampling density rather than of the city. Anchoring each label on the patch centre removes the double counting entirely and is exact at the scale lczkit validates on: on the Berlin fixture 438 patch centres fall into 438 distinct cells with no cell receiving two labels.

The 1:1 property holds, but not for the reason first recorded here. The patch centres are not aligned to the local UTM origin: measured on the Berlin fixture in EPSG:32633 they sit on an exact 100 m stride at a fixed phase offset of (40.0, 70.0) m from the GridUnits cell corners. That offset is what makes the property robust - being far from both 0 and 50 m, no centre can land on a cell boundary and no two centres can share a cell. An offset near either value would degrade, which is precisely what LabelMatch exists to expose, so it is reported per run rather than assumed.

What the centre rule does not control is support. The label describes a 320 x 320 m patch, 10.24 ha, and is attributed to one 1 ha cell whose centre is systematically ~22 m from the patch centre. A 100 m cell inside a compact-midrise patch can legitimately be a courtyard. That is an irreducible floor under any agreement figure measured this way, and it is the same patch-versus-cell mismatch the published parameter ranges show - the labels are patch-scale objects too.

The reduction is deliberately not raster-based. reference_lcz reads a categorical raster through LocalRasterSource because the Demuzere map is one; these patches are vector polygons whose value lies in their exact placement, and rasterising them to reuse that path would reintroduce the overlap problem in a less visible form.

Returns the same three columns as lczkit.validation.reference.reference_lcz, so agreement() consumes either without knowing which it was given.

LABEL_COLUMN module-attribute

LABEL_COLUMN = 'LCZ_class'

So2Sat LCZ42's own class column. Integer 1-17, matching Demuzere's coding for 1-10 and A-G, so no translation is needed - asserted by tests/test_validation_labelled.py.

SO2SAT_CITATION module-attribute

SO2SAT_CITATION = '10.1109/MGRS.2020.2964708'

Zhu et al. (2020), IEEE GRSM 8(3), 76-89. So2Sat LCZ42.

LabelMatch dataclass

LabelMatch(n_patches: int, n_centres_matched: int, n_centres_unmatched: int, n_centres_ambiguous: int, n_units_labelled: int, n_units_multi_label: int)

How the patch centres landed on the units, so a misaligned city is visible in the output.

The 1:1 mapping above holds because a So2Sat patch grid and a 100 m GridUnits grid share the UTM origin. Nothing guarantees that for every city or every unit strategy, and a silent degradation - centres falling on cell boundaries, or several centres per unit disagreeing - would look exactly like a well-measured run. These counts make it look like what it is.

n_centres_matched instance-attribute

n_centres_matched: int

Centres falling strictly inside exactly one unit.

n_centres_unmatched instance-attribute

n_centres_unmatched: int

Centres outside every unit, or on a shared boundary and so inside none of them.

n_centres_ambiguous instance-attribute

n_centres_ambiguous: int

Centres matching more than one unit. Zero for any valid partition; non-zero means the units overlap, and those centres are discarded rather than counted twice.

n_units_multi_label instance-attribute

n_units_multi_label: int

Units receiving centres of more than one class, and so decided by a majority rather than by a single label. Zero on a 100 m grid; non-zero on enclosures, which are larger than a patch.

labelled_lcz

labelled_lcz(units: GeoDataFrame, patches: GeoDataFrame, *, class_column: str = LABEL_COLUMN) -> tuple[DataFrame, LabelMatch]

The ground-truth class per unit_id, and a record of how the labels were matched.

Returns three columns, matching reference_lcz so the two are interchangeable downstream:

• reference_lcz - the majority class among the patch centres inside the unit, or null where no centre falls in it. Nullable Int8, never a sentinel.
• reference_coverage - 1.0 where the unit carries at least one centre, else 0.0. Labelled patches are a sample, not a map: a unit either has a label or has none, and reporting a fractional coverage would invite min_reference_coverage to filter on a number that means something different from the one it was written for.
• reference_majority_fraction - share of the unit's centres holding the winning class. 1.0 wherever a unit holds one centre, which on a 100 m grid is every labelled unit.

Neither input is mutated.

Source code in src/lczkit/validation/labelled.py

def labelled_lcz(
    units: gpd.GeoDataFrame,
    patches: gpd.GeoDataFrame,
    *,
    class_column: str = LABEL_COLUMN,
) -> tuple[pd.DataFrame, LabelMatch]:
    """The ground-truth class per `unit_id`, and a record of how the labels were matched.

    Returns three columns, matching `reference_lcz` so the two are interchangeable downstream:

    - `reference_lcz` - the majority class among the patch centres inside the unit, or null where
      no centre falls in it. Nullable `Int8`, never a sentinel.
    - `reference_coverage` - 1.0 where the unit carries at least one centre, else 0.0. Labelled
      patches are a *sample*, not a map: a unit either has a label or has none, and reporting a
      fractional coverage would invite `min_reference_coverage` to filter on a number that means
      something different from the one it was written for.
    - `reference_majority_fraction` - share of the unit's centres holding the winning class. 1.0
      wherever a unit holds one centre, which on a 100 m grid is every labelled unit.

    Neither input is mutated.
    """
    check_units(units)
    if class_column not in patches.columns:
        raise ValueError(
            f"patches must carry a {class_column!r} column; got {list(patches.columns)}"
        )
    if patches.crs is None:
        raise ValueError("patches must declare a CRS")
    target = units.crs
    assert target is not None  # noqa: S101 - check_units already raised if it were
    # Reprojected before the centroid is taken, never after: the centroid of a lat/lon polygon is
    # not the projection of the centroid, and this one decides which cell the label lands in.
    patches = patches.to_crs(target)

    centres = gpd.GeoDataFrame(
        {class_column: patches[class_column].to_numpy()},
        geometry=patches.geometry.centroid,
        crs=target,
    )
    # `within` rather than `intersects`: a centre on a shared cell boundary belongs to neither
    # cell, and counting it in both would let one patch label two units.
    matched = centres.sjoin(units[["geometry"]], predicate="within", how="left")

    ambiguous = matched.index.duplicated(keep=False) & matched["unit_id"].notna()
    unmatched = matched["unit_id"].isna()
    usable = matched.loc[~ambiguous & ~unmatched]

    votes = (
        usable.groupby(["unit_id", class_column]).size().rename("n").reset_index()
        if not usable.empty
        else pd.DataFrame(columns=["unit_id", class_column, "n"])
    )
    result = pd.DataFrame(
        {
            "reference_lcz": pd.Series(dtype="Int8"),
            "reference_coverage": 0.0,
            "reference_majority_fraction": pd.Series(dtype="float64"),
        },
        index=units.index,
    )
    n_multi = 0
    if not votes.empty:
        totals = votes.groupby("unit_id")["n"].sum()
        # `idxmax` on the group, not `.first()`/`.last()`: those skip nulls per column and would
        # silently return a class the unit's largest vote block does not hold.
        winners = votes.loc[votes.groupby("unit_id")["n"].idxmax()].set_index("unit_id")
        result["reference_lcz"] = winners[class_column].reindex(units.index).astype("Int8")
        result["reference_coverage"] = totals.reindex(units.index).notna().astype("float64")
        result["reference_majority_fraction"] = (
            winners["n"].div(totals).reindex(units.index).astype("float64")
        )
        n_multi = int((votes.groupby("unit_id").size() > 1).sum())

    match = LabelMatch(
        n_patches=int(len(patches)),
        n_centres_matched=int(len(usable)),
        n_centres_unmatched=int(unmatched.sum()),
        n_centres_ambiguous=int(matched.loc[ambiguous].index.nunique()),
        n_units_labelled=int(result["reference_lcz"].notna().sum()),
        n_units_multi_label=n_multi,
    )
    return result, match

lczkit.validation.reference

The reference LCZ map, reduced to one label per spatial unit.

The comparator map is the Demuzere global LCZ map, compared on the 100 m grid. The reduction that needs is areal class fractions per unit followed by the majority - which is exactly what LocalRasterSource already does for land cover, down to the CRS handling, the covering window, the exact cell weighting and the nodata policy. So the reference map is described as a LandCoverDatasetConfig and read through that source rather than growing a second exactextract path that could drift from the first.

The one adaptation is that nodata is assigned to its own class instead of being excluded, so the fractions still say how much of each unit the reference map reaches. Without that a unit lying half outside the map would report a confident majority computed from a corner of itself, and the agreement statistics would quietly include it.

"The 100 m grid" describes lczkit's units, not the reference. lcz_v3.tif is distributed in EPSG:4326 at 0.000898315 degrees, which is ~100 m north-south and 100 m east-west only at the equator: at Berlin's latitude the cells are 60.9 m wide, at Rotterdam 61.7 m, at Hong Kong 92.5 m. The lczkit unit is a 100 m square in a projected CRS, so the two grids share neither an origin nor an east-west cell size, and each unit takes the areal-majority label over roughly one-and-a-half reference cells across by one down.

Nothing here is approximated by that - units are reprojected into the raster's CRS rather than warping a categorical raster, and exactextract weights each cell by its exact covered fraction. But the majority is a reduction, and reference_majority_fraction reports how decisively it was won. A unit whose winner holds 0.34 of the observed area is close to a coin toss between three classes and is currently weighted the same as one holding 1.0.

NODATA_CLASS module-attribute

NODATA_CLASS = 'nodata'

Class name the reference dataset config assigns unmapped/nodata cells to.

reference_lcz

reference_lcz(units: GeoDataFrame, path: Path, config: LandCoverDatasetConfig, *, max_raster_cells: int = 200000000) -> DataFrame

The majority reference class per unit_id, with how well it is supported.

Returns three columns:

• reference_lcz - the integer code covering most of the unit, or null where the reference map covers none of it. Nullable Int8, never a sentinel: "the reference does not reach here" and "the reference says class 0" must not look the same.
• reference_coverage - fraction of the unit the reference map actually observed.
• reference_majority_fraction - fraction of the observed part the majority class holds. A unit split evenly between two classes agrees with whichever it is compared against about half the time, and that is worth knowing before reading a confusion matrix.

Neither input is mutated.

Source code in src/lczkit/validation/reference.py

def reference_lcz(
    units: gpd.GeoDataFrame,
    path: Path,
    config: LandCoverDatasetConfig,
    *,
    max_raster_cells: int = 200_000_000,
) -> pd.DataFrame:
    """The majority reference class per `unit_id`, with how well it is supported.

    Returns three columns:

    - `reference_lcz` - the integer code covering most of the unit, or null where the reference
      map covers none of it. Nullable `Int8`, never a sentinel: "the reference does not reach
      here" and "the reference says class 0" must not look the same.
    - `reference_coverage` - fraction of the unit the reference map actually observed.
    - `reference_majority_fraction` - fraction of the *observed* part the majority class holds.
      A unit split evenly between two classes agrees with whichever it is compared against about
      half the time, and that is worth knowing before reading a confusion matrix.

    Neither input is mutated.
    """
    check_units(units)
    if config.nodata_class != NODATA_CLASS:
        raise ValueError(
            f"reference dataset must assign nodata to a {NODATA_CLASS!r} class so unit coverage "
            f"is measurable; got nodata_class={config.nodata_class!r}"
        )

    fractions = LocalRasterSource(config, path, max_raster_cells=max_raster_cells).fractions(units)
    class_columns = [f"{config.column_prefix}lcz_{code}" for code in CODES]
    nodata_column = f"{config.column_prefix}{NODATA_CLASS}"

    observed = 1.0 - fractions[nodata_column]
    shares = fractions[class_columns]
    best = shares.max(axis=1)
    # `idxmax` would return a column name and a spurious winner for an all-zero row, so the
    # position is taken directly and masked where nothing was observed at all.
    winner = pd.Series(
        np.asarray(CODES)[shares.to_numpy(dtype="float64", na_value=0.0).argmax(axis=1)],
        index=fractions.index,
        dtype="float64",
    ).where(best > 0)

    return pd.DataFrame(
        {
            "reference_lcz": winner.astype("Int8"),
            "reference_coverage": observed,
            "reference_majority_fraction": best.div(observed.where(observed > 0)),
        },
        index=units.index,
    )

lczkit.validation.wudapt

WUDAPT LCZ training areas, reduced to one label per spatial unit.

The third reference, and the only one that reaches every city this package has been run on. It is the secondary reference, and the first where So2Sat has too few labels for a region of interest. Measured over sixteen study windows it carries 119-2374 polygons each, 11-18 classes, and 63-996 km2 of labelled ground against So2Sat's 7.1 km2 union on Berlin - one to two orders of magnitude more labelled area, at neighbourhood scale. (lczkit.cities now holds twenty; the four added afterwards are not in that range, and every figure quoted in this module is over the original sixteen.)

Why the reduction is areal here and not centroid-anchored. labelled_lcz anchors each So2Sat label on its patch centre, and that rule is justified by a property WUDAPT does not have: So2Sat patches are uniform 320 m squares on a 100 m stride, so a centroid is an exact, non-double-counting stand-in for the patch. WUDAPT polygons are hand-drawn and span 0 m2 to 18 680 km2, with a median of 4.8 ha across the whole file. A centroid rule would let a 4.8 ha polygon and a 1000 km2 polygon each label exactly one unit, which is not a sampling scheme but a discarding of almost all of the reference. So the reduction is an areal overlay, and reference_coverage is genuinely fractional here - unlike labelled_lcz, where it is deliberately binary.

The overlaps are real and are resolved before the overlay, not during it. Contributions come from different submissions in different years - representative_date spans 1983 to 2025 - and they overlap each other: in a 15 x 12 km Kowloon window, 803 polygons form 3330 overlapping pairs of which 560 carry different classes. Overlaying that directly would count the same ground several times under labels that contradict each other, and the "majority" would measure contributor enthusiasm rather than the city. resolve_overlaps gives each piece of ground to exactly one polygon by an explicit priority, and reports how much ground was contested and how much was merely duplicated - the contributors' own disagreement rate is a number worth having rather than hiding.

Two things about this file that will mislead a reader who trusts it.

• The stored area column is not area in any usable unit: it is km2 computed in Web Mercator, so it is inflated by 1/cos^2(latitude) - the median ratio to true area is 1 004 995 against Mollweide's 744 899. Nothing here reads it; test_validation_wudapt.py asserts that.
• class runs 1 to 19, not 1 to 17. 633 polygons globally carry codes 18 and 19, which are outside the Demuzere/So2Sat coding this package uses everywhere else. They are dropped and counted, never mapped onto a neighbouring class.

WUDAPT is not independent of lcz_v3. The LCZ Generator's training areas are the training data behind the Demuzere global map, so an agreement figure between the two is inflated by construction and is not a ceiling in the sense So2Sat gives one. Report it only with that stated beside it.

Returns the same three columns as reference_lcz and labelled_lcz, so agreement() consumes any of the three without knowing which it was given.

WUDAPT_SOURCE_DIR_NAME module-attribute

WUDAPT_SOURCE_DIR_NAME = 'WUDAPT'

input/<name>/ holding the LCZ Generator training-area export.

WUDAPT_CITATION module-attribute

WUDAPT_CITATION = '10.3390/ijgi4010199'

Bechtel et al. (2015), IJGI 4(1), 199-219 - the WUDAPT Level 0 protocol these areas follow. The LCZ Generator itself is Demuzere, Kittner & Bechtel (2021), 10.3389/fenvs.2021.637455.

CLASS_COLUMN module-attribute

CLASS_COLUMN = 'class'

WUDAPT's own class column. Integer, and 1-17 agrees with Demuzere's coding - but see UNSUPPORTED_CLASSES, because it does not stop at 17.

UNSUPPORTED_CLASSES module-attribute

UNSUPPORTED_CLASSES: tuple[int, ...] = (18, 19)

Codes present in the file and absent from lczkit.classify.labels.CODES.

633 polygons of 630 311 globally. They are dropped rather than folded into a neighbouring class: this package has no definition for them, and inventing one would put a label into the reference that no contributor ever drew.

QC_COLUMNS module-attribute

QC_COLUMNS: tuple[str, ...] = ('qc_step1', 'qc_step2', 'qc_step3')

The LCZ Generator's three submission quality-control flags, stored as strings.

PRIORITY_COLUMNS module-attribute

PRIORITY_COLUMNS: tuple[str, ...] = ('representative_date', 'oa', 'submission_date')

Columns resolve_overlaps reads to rank contested ground. Absent columns are skipped.

READ_COLUMNS module-attribute

READ_COLUMNS: tuple[str, ...] = (CLASS_COLUMN, *QC_COLUMNS, 'submission_id', 'submission_date', 'representative_date', 'city', 'oa', 'license')

Everything the loader needs. Reading 40 columns over 630 311 features to use nine of them is the difference between a fixture build that takes seconds and one that takes minutes.

OVERLAP_EPS_M2 module-attribute

OVERLAP_EPS_M2 = 1e-06

Intersection area below which two resolved polygons are treated as not overlapping.

shapely.difference leaves coordinate-noise slivers along the cut, so the topological overlaps predicate still fires on pairs whose shared area is around 1e-8 m² — a hundredth of a square micrometre, measured on the Berlin fixture. The resolution is exact in area (labelled + duplicate + conflict reproduces the raw sum to the full 13 595 047.0 m² on that fixture); it is the predicate that is not a statement about area. Tests assert on area against this tolerance.

A module constant and not config, for the same reason eps_m is one: a floating-point tolerance is not a domain threshold.

WudaptSelection dataclass

WudaptSelection(n_read: int, n_kept: int, n_dropped_invalid: int, n_dropped_unsupported_class: int, n_dropped_qc: int, n_dropped_accuracy: int, n_dropped_area: int, qc_pass_fraction: float, n_overlapping_pairs: int, n_conflicting_pairs: int, duplicate_area_m2: float, conflict_area_m2: float, labelled_area_m2: float, date_min: str | None, date_max: str | None, licences: tuple[str, ...] = (), cities: tuple[str, ...] = ())

What prepare_wudapt kept, what it dropped, and what the contributors disagreed about.

Every drop is counted rather than being visible only as a smaller frame, for the same reason the cleaning report counts what it removes: a reference that quietly loses half its polygons to a quality gate looks exactly like one that never had them.

n_dropped_invalid instance-attribute

n_dropped_invalid: int

Non-polygonal or unrepairable geometry. About 1.0% of the file self-intersects.

n_dropped_unsupported_class instance-attribute

n_dropped_unsupported_class: int

Codes in UNSUPPORTED_CLASSES, or null.

n_dropped_qc instance-attribute

n_dropped_qc: int

Failed the quality gate. Zero when WudaptConfig.require_qc is off, which is the default - the gate costs 51.8% of the file, so it is offered rather than imposed.

n_dropped_accuracy instance-attribute

n_dropped_accuracy: int

Below WudaptConfig.min_oa. oa is a property of the submission, not the polygon.

n_dropped_area instance-attribute

n_dropped_area: int

Outside WudaptConfig.min_area_m2 / max_area_m2.

qc_pass_fraction instance-attribute

qc_pass_fraction: float

Share of the read polygons passing all three QC flags, reported whether or not the gate is on, so the cost of turning it on is visible without a second run.

n_overlapping_pairs instance-attribute

n_overlapping_pairs: int

Pairs of kept polygons whose interiors overlap, in one direction each.

n_conflicting_pairs instance-attribute

n_conflicting_pairs: int

Of those, the pairs carrying different classes - contributors disagreeing about ground, as opposed to re-drawing it.

duplicate_area_m2 instance-attribute

duplicate_area_m2: float

Ground a polygon yielded to a higher-priority polygon of the same class. Redundancy.

conflict_area_m2 instance-attribute

conflict_area_m2: float

Ground a polygon yielded to a higher-priority polygon of a different class. This is the reference disagreeing with itself, and it belongs in the write-up beside any figure measured against it.

labelled_area_m2 instance-attribute

labelled_area_m2: float

Area of the resolved, non-overlapping polygons. The reference's actual support.

date_max instance-attribute

date_max: str | None

Range of representative_date among the kept polygons. A reference spanning four decades is not a snapshot, and a city that redeveloped inside that span will disagree with any single map for reasons that are not the map's error.

licences class-attribute instance-attribute

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

Distinct license values among the kept polygons, read from the data rather than assumed. Most of the file is CC BY-NC-SA 4.0 or CC BY-SA - non-commercial in the first case. That constrains the data, not this MIT package, and nothing here redistributes it; it is recorded so a run's manifest states the terms of what it was scored against.

cities class-attribute instance-attribute

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

Distinct city strings among the kept polygons. Free text and not normalised upstream - Wuhan and wuhan are different submissions - so this is provenance, never a join key.

WudaptMatch dataclass

WudaptMatch(n_polygons: int, n_units: int, n_units_labelled: int, n_units_multi_label: int, unit_area_m2: float, labelled_area_m2: float, mean_coverage: float, class_counts: dict[int, int] = dict())

How the resolved polygons landed on the units.

The analogue of LabelMatch, and it exists for the same reason: a reference that reaches a tenth of the study area produces exactly the same column set as one that reaches all of it.

n_units_labelled instance-attribute

n_units_labelled: int

Units receiving any labelled ground at all, before min_reference_coverage filters them.

n_units_multi_label instance-attribute

n_units_multi_label: int

Units whose labelled ground carries more than one class, and so are decided by an areal majority. Expected to be common on patch-scale units and rare on a 100 m grid.

labelled_area_m2 instance-attribute

labelled_area_m2: float

Labelled ground falling inside the units. Below WudaptSelection.labelled_area_m2 whenever the reference extends past the study window, which it usually does.

mean_coverage instance-attribute

mean_coverage: float

Area-weighted mean reference_coverage over the labelled units.

class_counts class-attribute instance-attribute

class_counts: dict[int, int] = field(default_factory=dict)

Units per assigned reference class, so a window carrying two classes is visible as such.

OverlapReport dataclass

OverlapReport(n_overlapping_pairs: int, n_conflicting_pairs: int, duplicate_area_m2: float, conflict_area_m2: float)

What resolve_overlaps had to arbitrate.

The three areas satisfy raw = labelled + duplicate + conflict exactly, which is the invariant worth checking: every square metre a contributor drew is either kept, redundant, or contested, and nothing is silently lost.

n_overlapping_pairs instance-attribute

n_overlapping_pairs: int

Ordered pairs, counted once per lower-priority polygon per higher-priority polygon it genuinely shares area with. Polygons that merely touch are excluded.

n_conflicting_pairs instance-attribute

n_conflicting_pairs: int

Of those, the pairs whose classes differ.

duplicate_area_m2 instance-attribute

duplicate_area_m2: float

Ground yielded to a higher-priority polygon of the same class, summed over polygons — so ground claimed by k agreeing polygons contributes k-1 times. It measures redundant drawing effort, not an area of the city.

conflict_area_m2 instance-attribute

conflict_area_m2: float

Ground yielded to a higher-priority polygon of a different class, summed the same way. Attributed to duplicate_area_m2 first where a polygon lost the same ground to both an agreeing and a disagreeing claimant, so this is the conservative reading of the disagreement.

read_wudapt

read_wudapt(path: Path, bbox: BBox, *, layer: str | None = None, columns: tuple[str, ...] = READ_COLUMNS) -> GeoDataFrame

The raw polygons intersecting bbox, in the file's own CRS (EPSG:4326).

Reads only columns and only the bbox: the file is 720 MB and 630 311 features, and every caller wants a city-sized window of it.

layer=None resolves to the first geometry-bearing layer rather than being passed through. The published export carries a second layer, layer_styles, which is a QGIS style table with no geometry; letting the driver pick would work today and read the style table the day the layer order changes.

Source code in src/lczkit/validation/wudapt.py

def read_wudapt(
    path: Path,
    bbox: BBox,
    *,
    layer: str | None = None,
    columns: tuple[str, ...] = READ_COLUMNS,
) -> gpd.GeoDataFrame:
    """The raw polygons intersecting `bbox`, in the file's own CRS (EPSG:4326).

    Reads only `columns` and only the bbox: the file is 720 MB and 630 311 features, and every
    caller wants a city-sized window of it.

    `layer=None` resolves to the first *geometry-bearing* layer rather than being passed through.
    The published export carries a second layer, `layer_styles`, which is a QGIS style table with
    no geometry; letting the driver pick would work today and read the style table the day the
    layer order changes.
    """
    if not path.exists():
        raise FileNotFoundError(
            f"WUDAPT training areas not found at {path}. Place the LCZ Generator export under "
            f"input/{WUDAPT_SOURCE_DIR_NAME}/ and set ValidationConfig.wudapt.filename."
        )
    frame = gpd.read_file(
        path,
        layer=layer if layer is not None else _first_geometry_layer(path),
        bbox=bbox,
        columns=list(columns),
    )
    return gpd.GeoDataFrame(frame)

prepare_wudapt

prepare_wudapt(polygons: GeoDataFrame, *, crs: CRS, config: WudaptConfig | None = None) -> tuple[GeoDataFrame, WudaptSelection]

Clean, filter and de-overlap polygons, returning them in crs with a record of the cost.

The returned frame carries class and a non-overlapping geometry: every piece of ground belongs to exactly one polygon, so an areal overlay onto units cannot double-count it.

polygons is not mutated.

Source code in src/lczkit/validation/wudapt.py

def prepare_wudapt(
    polygons: gpd.GeoDataFrame,
    *,
    crs: CRS,
    config: WudaptConfig | None = None,
) -> tuple[gpd.GeoDataFrame, WudaptSelection]:
    """Clean, filter and de-overlap `polygons`, returning them in `crs` with a record of the cost.

    The returned frame carries `class` and a non-overlapping geometry: every piece of ground
    belongs to exactly one polygon, so an areal overlay onto units cannot double-count it.

    `polygons` is not mutated.
    """
    settings = config or WudaptConfig()
    class_column = settings.class_column
    n_read = int(len(polygons))
    if class_column not in polygons.columns:
        raise ValueError(
            f"polygons must carry a {class_column!r} column; got {list(polygons.columns)}"
        )
    if polygons.crs is None:
        raise ValueError("polygons must declare a CRS")

    frame = polygons.to_crs(crs)
    qc_pass = _qc_pass(frame)
    qc_pass_fraction = float(qc_pass.fillna(False).mean()) if n_read else 0.0

    # Force 2D before anything geometric: every feature in the file is `Polygon Z` carrying a
    # dummy Z, and a Z coordinate survives `make_valid` into overlay operations that do not want it.
    geometry = shapely.force_2d(shapely.make_valid(frame.geometry.to_numpy()))
    frame = frame.set_geometry(gpd.GeoSeries(geometry, index=frame.index, crs=crs))
    valid = frame.geom_type.isin(("Polygon", "MultiPolygon")) & ~frame.geometry.is_empty
    n_dropped_invalid = int((~valid).sum())
    frame = frame.loc[valid]

    codes = pd.to_numeric(frame[class_column], errors="coerce")
    supported = codes.isin(CODES)
    n_dropped_unsupported = int((~supported).sum())
    frame = frame.loc[supported].assign(**{class_column: codes[supported].astype("int64")})

    kept_qc = qc_pass.reindex(frame.index)
    if settings.require_qc:
        gate = kept_qc.fillna(False).astype(bool)
        n_dropped_qc = int((~gate).sum())
        frame = frame.loc[gate]
    else:
        n_dropped_qc = 0

    if settings.min_oa is not None and "oa" in frame.columns:
        accurate = pd.to_numeric(frame["oa"], errors="coerce").ge(settings.min_oa).fillna(False)
        n_dropped_accuracy = int((~accurate).sum())
        frame = frame.loc[accurate.astype(bool)]
    else:
        n_dropped_accuracy = 0

    # Recomputed in the projected CRS, never read from the stored `area` column - see the module
    # docstring. `test_validation_wudapt.py` asserts the stored column is not consulted.
    areas = frame.geometry.area
    big_enough = areas.ge(settings.min_area_m2)
    if settings.max_area_m2 is not None:
        big_enough &= areas.le(settings.max_area_m2)
    n_dropped_area = int((~big_enough).sum())
    frame = frame.loc[big_enough]

    resolved, overlaps = resolve_overlaps(frame, class_column=class_column)

    dates = _dates(frame["representative_date"]) if "representative_date" in frame else None
    return resolved, WudaptSelection(
        n_read=n_read,
        n_kept=int(len(resolved)),
        n_dropped_invalid=n_dropped_invalid,
        n_dropped_unsupported_class=n_dropped_unsupported,
        n_dropped_qc=n_dropped_qc,
        n_dropped_accuracy=n_dropped_accuracy,
        n_dropped_area=n_dropped_area,
        qc_pass_fraction=qc_pass_fraction,
        n_overlapping_pairs=overlaps.n_overlapping_pairs,
        n_conflicting_pairs=overlaps.n_conflicting_pairs,
        duplicate_area_m2=overlaps.duplicate_area_m2,
        conflict_area_m2=overlaps.conflict_area_m2,
        labelled_area_m2=float(resolved.geometry.area.sum()),
        date_min=_stamp(dates.min()) if dates is not None and dates.notna().any() else None,
        date_max=_stamp(dates.max()) if dates is not None and dates.notna().any() else None,
        licences=_distinct(frame, "license"),
        cities=_distinct(frame, "city"),
    )

priority_order

priority_order(polygons: GeoDataFrame) -> Index

Index of polygons, best claim to contested ground first.

Most recent representative_date, then the higher submission accuracy, then the smaller polygon, then the index. The last two are what make it deterministic; the first two are the judgement:

• Recency wins because a WUDAPT polygon describes the city at representative_date, the dates span 1983-2025, and lczkit is classifying a current Overture release. Where two contributors disagree about the same ground, the newer one is describing the city the package is looking at.
• Accuracy breaks the tie because oa is the LCZ Generator's own cross-validated score for the submission the polygon came from. It is a property of the submission, not the polygon, so it is a weak signal - which is exactly why it ranks below recency rather than above it.
• Smaller wins last, on the reasoning that a contributor who drew a 2 ha polygon inside someone else's 200 ha one was being more specific about that ground.

A missing date or accuracy sorts last rather than dropping the polygon: it costs the polygon its claim on contested ground, not its label on ground nobody else drew.

Source code in src/lczkit/validation/wudapt.py

def priority_order(polygons: gpd.GeoDataFrame) -> pd.Index:
    """Index of `polygons`, best claim to contested ground first.

    Most recent `representative_date`, then the higher submission accuracy, then the *smaller*
    polygon, then the index. The last two are what make it deterministic; the first two are the
    judgement:

    - **Recency wins** because a WUDAPT polygon describes the city at `representative_date`, the
      dates span 1983-2025, and lczkit is classifying a current Overture release. Where two
      contributors disagree about the same ground, the newer one is describing the city the
      package is looking at.
    - **Accuracy breaks the tie** because `oa` is the LCZ Generator's own cross-validated score
      for the submission the polygon came from. It is a property of the submission, not the
      polygon, so it is a weak signal - which is exactly why it ranks below recency rather than
      above it.
    - **Smaller wins** last, on the reasoning that a contributor who drew a 2 ha polygon inside
      someone else's 200 ha one was being more specific about that ground.

    A missing date or accuracy sorts last rather than dropping the polygon: it costs the polygon
    its claim on contested ground, not its label on ground nobody else drew.
    """
    frame = pd.DataFrame(index=polygons.index)
    if "representative_date" in polygons.columns:
        frame["date"] = _dates(polygons["representative_date"])
    elif "submission_date" in polygons.columns:
        frame["date"] = pd.to_datetime(polygons["submission_date"], errors="coerce", utc=True)
    if "oa" in polygons.columns:
        frame["oa"] = pd.to_numeric(polygons["oa"], errors="coerce")
    frame["small"] = -polygons.geometry.area

    by = [*(column for column in ("date", "oa") if column in frame.columns), "small"]
    ordered = frame.sort_values(by=by, ascending=False, na_position="last", kind="stable")
    return ordered.index

resolve_overlaps

resolve_overlaps(polygons: GeoDataFrame, *, class_column: str = CLASS_COLUMN) -> tuple[GeoDataFrame, OverlapReport]

Give every piece of ground to exactly one polygon, and report what that cost.

Each polygon keeps itself minus every higher-priority polygon. Subtracting the higher-priority polygons' original geometries rather than their remainders is equivalent and cheaper: ground a higher-priority polygon itself lost went to something higher still, so it is claimed either way and never returns to a lower-priority claimant.

Polygons reduced to nothing are dropped. Their label is not lost - it agreed with, or was overruled by, a polygon that kept the same ground.

The result is exact in area but not in topology: the cuts leave coordinate-noise slivers, so overlaps still fires on some pairs at around 1e-8 m2. See OVERLAP_EPS_M2.

Source code in src/lczkit/validation/wudapt.py

def resolve_overlaps(
    polygons: gpd.GeoDataFrame, *, class_column: str = CLASS_COLUMN
) -> tuple[gpd.GeoDataFrame, OverlapReport]:
    """Give every piece of ground to exactly one polygon, and report what that cost.

    Each polygon keeps itself minus every higher-priority polygon. Subtracting the higher-priority
    polygons' *original* geometries rather than their remainders is equivalent and cheaper: ground
    a higher-priority polygon itself lost went to something higher still, so it is claimed either
    way and never returns to a lower-priority claimant.

    Polygons reduced to nothing are dropped. Their label is not lost - it agreed with, or was
    overruled by, a polygon that kept the same ground.

    The result is exact in area but not in topology: the cuts leave coordinate-noise slivers, so
    `overlaps` still fires on some pairs at around 1e-8 m2. See `OVERLAP_EPS_M2`.
    """
    if polygons.empty:
        return polygons, OverlapReport(0, 0, 0.0, 0.0)

    order = priority_order(polygons)
    ordered = polygons.loc[order]
    geometry = ordered.geometry.to_numpy()
    codes = ordered[class_column].to_numpy()
    tree = shapely.STRtree(geometry)

    remainders: list[shapely.Geometry] = []
    duplicate_area = 0.0
    conflict_area = 0.0
    overlapping = 0
    conflicting = 0

    for position, geom in enumerate(geometry):
        hits = tree.query(geom, predicate="intersects")
        higher = hits[hits < position]
        if higher.size == 0:
            remainders.append(geom)
            continue

        # `intersects` includes polygons that merely touch, which take no area and are not
        # overlaps. Counting them would inflate the disagreement rate with shared edges.
        claimed = shapely.intersection(geom, shapely.union_all(geometry[higher]))
        lost = float(shapely.area(claimed))
        if lost > 0.0:
            same = higher[codes[higher] == codes[position]]
            different = higher.size - same.size
            overlapping += int(higher.size)
            conflicting += int(different)
            duplicate = (
                float(shapely.area(shapely.intersection(geom, shapely.union_all(geometry[same]))))
                if same.size
                else 0.0
            )
            duplicate_area += duplicate
            conflict_area += max(lost - duplicate, 0.0)
        remainders.append(shapely.difference(geom, shapely.union_all(geometry[higher])))

    cut = np.asarray(remainders, dtype=object)
    resolved = ordered.set_geometry(gpd.GeoSeries(cut, index=ordered.index, crs=polygons.crs))
    # Tested on the shapely array rather than through `GeoSeries.notna`, whose treatment of empty
    # geometries changed in geopandas 1.x and now warns whenever any are present — which is exactly
    # the case here, since a fully absorbed polygon is how this function reports "kept nothing".
    kept = ~(shapely.is_empty(cut) | shapely.is_missing(cut))
    return gpd.GeoDataFrame(resolved.loc[kept]), OverlapReport(
        n_overlapping_pairs=overlapping,
        n_conflicting_pairs=conflicting,
        duplicate_area_m2=duplicate_area,
        conflict_area_m2=conflict_area,
    )

wudapt_lcz

wudapt_lcz(units: GeoDataFrame, polygons: GeoDataFrame, *, class_column: str = CLASS_COLUMN) -> tuple[DataFrame, WudaptMatch]

The reference class per unit_id by areal majority, and how the labels landed.

Returns the three columns reference_lcz and labelled_lcz return, with the same meanings:

• reference_lcz - the class holding most of the unit's labelled area, or null where no polygon reaches it. Nullable Int8, never a sentinel.
• reference_coverage - fraction of the unit any polygon covers. Genuinely fractional, unlike labelled_lcz's binary flag, because WUDAPT polygons tile ground rather than sampling it.
• reference_majority_fraction - the winner's share of the labelled part.

polygons must already be non-overlapping - pass what prepare_wudapt returns. Overlapping input is not detected here (the check costs as much as the fix); it would inflate every coverage above 1.0 and let one contributor outvote a city.

Neither input is mutated.

Source code in src/lczkit/validation/wudapt.py

def wudapt_lcz(
    units: gpd.GeoDataFrame,
    polygons: gpd.GeoDataFrame,
    *,
    class_column: str = CLASS_COLUMN,
) -> tuple[pd.DataFrame, WudaptMatch]:
    """The reference class per `unit_id` by areal majority, and how the labels landed.

    Returns the three columns `reference_lcz` and `labelled_lcz` return, with the same meanings:

    - `reference_lcz` - the class holding most of the unit's *labelled* area, or null where no
      polygon reaches it. Nullable `Int8`, never a sentinel.
    - `reference_coverage` - fraction of the unit any polygon covers. Genuinely fractional, unlike
      `labelled_lcz`'s binary flag, because WUDAPT polygons tile ground rather than sampling it.
    - `reference_majority_fraction` - the winner's share of the labelled part.

    `polygons` must already be non-overlapping - pass what `prepare_wudapt` returns. Overlapping
    input is not detected here (the check costs as much as the fix); it would inflate every
    coverage above 1.0 and let one contributor outvote a city.

    Neither input is mutated.
    """
    check_units(units)
    if class_column not in polygons.columns:
        raise ValueError(
            f"polygons must carry a {class_column!r} column; got {list(polygons.columns)}"
        )
    if polygons.crs is None:
        raise ValueError("polygons must declare a CRS")

    unit_area = units.geometry.area
    result = pd.DataFrame(
        {
            "reference_lcz": pd.Series(pd.NA, index=units.index, dtype="Int8"),
            "reference_coverage": 0.0,
            "reference_majority_fraction": pd.Series(np.nan, index=units.index, dtype="float64"),
        },
        index=units.index,
    )
    empty = WudaptMatch(
        n_polygons=int(len(polygons)),
        n_units=int(len(units)),
        n_units_labelled=0,
        n_units_multi_label=0,
        unit_area_m2=float(unit_area.sum()),
        labelled_area_m2=0.0,
        mean_coverage=0.0,
    )
    if polygons.empty or units.empty:
        return result, empty

    target = units.crs
    assert target is not None  # noqa: S101 - check_units already raised if it were
    labels = polygons.to_crs(target)[[class_column, "geometry"]]
    pieces = gpd.overlay(
        units.reset_index()[["unit_id", "geometry"]],
        labels,
        how="intersection",
        keep_geom_type=True,
    )
    if pieces.empty:
        return result, empty

    pieces = pieces.assign(area=pieces.geometry.area)
    pieces = pieces.loc[pieces["area"] > 0.0]
    if pieces.empty:
        return result, empty

    by_class = pieces.groupby(["unit_id", class_column], observed=True)["area"].sum()
    totals = by_class.groupby("unit_id").sum()
    # `idxmax` over the (unit, class) index rather than a pivot: the pivot would be 17 columns wide
    # over every unit, and the winner is a single positional lookup either way.
    winners = by_class.loc[by_class.groupby("unit_id").idxmax()]
    winner_class = winners.index.get_level_values(class_column)
    winner_area = pd.Series(winners.to_numpy(), index=winners.index.get_level_values("unit_id"))
    winner_code = pd.Series(winner_class, index=winner_area.index)

    labelled = totals.reindex(units.index)
    result["reference_lcz"] = winner_code.reindex(units.index).astype("Int8")
    result["reference_coverage"] = (
        labelled.div(unit_area.where(unit_area > 0)).fillna(0.0).clip(upper=1.0)
    )
    result["reference_majority_fraction"] = winner_area.reindex(units.index).div(
        labelled.where(labelled > 0)
    )

    assigned = result["reference_lcz"].dropna()
    return result, WudaptMatch(
        n_polygons=int(len(polygons)),
        n_units=int(len(units)),
        n_units_labelled=int(len(assigned)),
        n_units_multi_label=int((by_class.groupby("unit_id").size() > 1).sum()),
        unit_area_m2=float(unit_area.sum()),
        labelled_area_m2=float(pieces["area"].sum()),
        mean_coverage=float(result.loc[assigned.index, "reference_coverage"].mean())
        if len(assigned)
        else 0.0,
        class_counts={
            int(str(code)): int(count) for code, count in assigned.value_counts().items()
        },
    )

load_wudapt

load_wudapt(path: Path, bbox: BBox, *, crs: CRS, config: WudaptConfig | None = None) -> tuple[GeoDataFrame, WudaptSelection]

read_wudapt then prepare_wudapt: the whole file route in one call.

Split into two functions underneath so tests can exercise the cleaning against a committed fixture frame without DATA_DIR being set.

Source code in src/lczkit/validation/wudapt.py

def load_wudapt(
    path: Path,
    bbox: BBox,
    *,
    crs: CRS,
    config: WudaptConfig | None = None,
) -> tuple[gpd.GeoDataFrame, WudaptSelection]:
    """`read_wudapt` then `prepare_wudapt`: the whole file route in one call.

    Split into two functions underneath so tests can exercise the cleaning against a committed
    fixture frame without `DATA_DIR` being set.
    """
    settings = config or WudaptConfig()
    raw = read_wudapt(path, bbox, layer=settings.layer)
    return prepare_wudapt(raw, crs=crs, config=settings)

Uncertainty

Confidence intervals come from a spatial-block bootstrap — resampling contiguous blocks of cells rather than individual cells. So2Sat patches are 320 m across on a 100 m grid, so neighbouring cells frequently carry the same label from the same patch: a city's labelled cells are one correlated sheet, and resampling them individually would report an interval far too narrow.

lczkit.validation.uncertainty

Confidence intervals for the agreement figures, by spatial block bootstrap.

Why an interval at all. Every headline this package reports - overall_agreement, built_agreement, both axis lifts - is a point estimate, and the sixteen-city work compares them across cities, differences them between arms, and orders the next lever by which is larger. None of that is readable without knowing how much of the difference is noise. axis_reconciliation.py compares max_over_median_lift against a hard < 5.0 as though it were a test statistic.

Why blocks rather than units. The units are not independent draws. So2Sat patches are 320 m squares sampled on a 100 m stride, so neighbouring patches overlap about sevenfold and the labelled cells of a city are one contiguous, strongly autocorrelated sheet. Resampling cells independently would treat 9 627 near-duplicate observations as 9 627 independent ones and return an interval far too narrow - the classic effective-sample-size error. Resampling contiguous blocks keeps the local correlation inside the resampling unit, which is what a block bootstrap is for.

The block size is a judgement, not a measurement: it has to be larger than the correlation length and small enough to leave many blocks. The default is 1 km, comfortably above the 320 m patch width that generates most of the dependence, and it is reported in the output so a reader can see what it was rather than infer it.

What this does not do. It quantifies sampling variability under the labels this run was scored against. It says nothing about the label support mismatch (a 320 m patch label attributed to a 1 ha cell), which is a bias rather than a variance and does not shrink with more cells.

Interval

Bases: BaseModel

A point estimate with a percentile bootstrap interval around it.

width property

width: float

upper - lower, in the units of the figure — the interval's span, not a half-width.

BootstrapReport

Bases: BaseModel

Intervals for the figures the project actually compares across cities and arms.

n_blocks instance-attribute

n_blocks: int

Blocks resampled with replacement. This, not the unit count, is the effective sample size - and the gap between the two is the whole reason this module exists.

spatial_blocks

spatial_blocks(units: GeoDataFrame, block_size_m: float) -> Series

Assign each unit to a square block of side block_size_m, keyed by unit_id.

Blocks are anchored on the projected CRS origin rather than on the data's own bounding box, for the reason GridUnits uses the same anchor: two runs over overlapping extents then agree about which block a unit is in, so their intervals are computed over comparable partitions.

Uses a representative point rather than a centroid so a concave enclosure lands in a block it actually occupies.

Source code in src/lczkit/validation/uncertainty.py

def spatial_blocks(units: gpd.GeoDataFrame, block_size_m: float) -> pd.Series:
    """Assign each unit to a square block of side `block_size_m`, keyed by `unit_id`.

    Blocks are anchored on the projected CRS origin rather than on the data's own bounding box, for
    the reason `GridUnits` uses the same anchor: two runs over overlapping extents then agree about
    which block a unit is in, so their intervals are computed over comparable partitions.

    Uses a representative point rather than a centroid so a concave enclosure lands in a block it
    actually occupies.
    """
    assert_projected_crs(units, "units")
    if block_size_m <= 0:
        raise ValueError(f"block_size_m must be positive; got {block_size_m}")

    points = units.geometry.representative_point()
    col = np.floor(points.x.to_numpy() / block_size_m).astype("int64")
    row = np.floor(points.y.to_numpy() / block_size_m).astype("int64")
    return pd.Series(
        [f"block_{c}_{r}" for c, r in zip(col, row, strict=True)],
        index=units.index,
        name="block",
    )

bootstrap_agreement

bootstrap_agreement(predicted: Series, reference: Series, area_m2: Series, blocks: Series, *, coverage: Series | None = None, config: ValidationConfig | None = None, n_resamples: int = 200, confidence: float = 0.95, seed: int = 0, block_size_m: float = float('nan')) -> BootstrapReport

Percentile bootstrap over spatial blocks, resampled with replacement.

blocks is a unit_id-indexed block label, from spatial_blocks(). Every input is aligned on unit_id first, so a caller cannot silently pair a prediction with another unit's block.

seed is fixed by default: an interval that moves between two runs of the same data is not a property of the data, and this package treats reproducibility as a feature rather than an afterthought. block_size_m is carried through for the record only - the blocks themselves arrive already formed, so nothing here depends on it being right, but a reported interval whose block size a reader has to guess is not much of a record.

Source code in src/lczkit/validation/uncertainty.py

def bootstrap_agreement(
    predicted: pd.Series,
    reference: pd.Series,
    area_m2: pd.Series,
    blocks: pd.Series,
    *,
    coverage: pd.Series | None = None,
    config: ValidationConfig | None = None,
    n_resamples: int = 200,
    confidence: float = 0.95,
    seed: int = 0,
    block_size_m: float = float("nan"),
) -> BootstrapReport:
    """Percentile bootstrap over spatial blocks, resampled with replacement.

    `blocks` is a `unit_id`-indexed block label, from `spatial_blocks()`. Every input is aligned on
    `unit_id` first, so a caller cannot silently pair a prediction with another unit's block.

    `seed` is fixed by default: an interval that moves between two runs of the same data is not a
    property of the data, and this package treats reproducibility as a feature rather than an
    afterthought. `block_size_m` is carried through for the record only - the blocks themselves
    arrive already formed, so nothing here depends on it being right, but a reported interval whose
    block size a reader has to guess is not much of a record.
    """
    if not 0.0 < confidence < 1.0:
        raise ValueError(f"confidence must lie in (0, 1); got {confidence}")
    if n_resamples < 1:
        raise ValueError(f"n_resamples must be at least 1; got {n_resamples}")

    frame = pd.DataFrame(
        {
            "predicted": predicted,
            "reference": reference,
            "area": area_m2,
            "block": blocks,
            "coverage": 1.0 if coverage is None else coverage,
        }
    ).dropna(subset=["block"])

    def score(sample: pd.DataFrame) -> dict[str, float]:
        """The tracked agreement figures for one resample, as a flat name-to-value mapping."""
        # A resample repeats blocks, so the index carries duplicate unit_ids. `agreement()` joins
        # its inputs by position through a DataFrame constructor, but a duplicated label would make
        # any later reindex ambiguous, so the resample gets a fresh unique index.
        indexed = sample.reset_index(drop=True)
        indexed.index = pd.Index([f"r{i}" for i in range(len(indexed))], name="unit_id")
        return _statistics(
            agreement(
                indexed["predicted"],
                indexed["reference"],
                indexed["area"],
                coverage=indexed["coverage"],
                config=config,
            )
        )

    point = score(frame)
    members = {name: group for name, group in frame.groupby("block", sort=True)}
    names = list(members)
    rng = np.random.default_rng(seed)

    draws: list[dict[str, float]] = []
    for _ in range(n_resamples):
        picked = rng.choice(len(names), size=len(names), replace=True)
        draws.append(score(pd.concat([members[names[i]] for i in picked])))

    tail = 100.0 * (1.0 - confidence) / 2.0

    def interval(name: str) -> Interval:
        """The percentile interval for one figure across the draws, around its point estimate.

        Non-finite draws are dropped rather than propagated: a resample can miss a class
        entirely, which makes that class's figure undefined for that draw and not for the city.
        An interval with no usable draws is returned as NaN bounds rather than invented.
        """
        values = np.asarray([draw[name] for draw in draws], dtype="float64")
        usable = values[np.isfinite(values)]
        if usable.size == 0:
            return Interval(point=point[name], lower=float("nan"), upper=float("nan"))
        return Interval(
            point=point[name],
            lower=float(np.percentile(usable, tail)),
            upper=float(np.percentile(usable, 100.0 - tail)),
        )

    return BootstrapReport(
        n_resamples=n_resamples,
        confidence=confidence,
        block_size_m=block_size_m,
        n_blocks=len(names),
        n_units=int(len(frame)),
        overall_agreement=interval("overall_agreement"),
        built_agreement=interval("built_agreement"),
        built_natural_agreement=interval("built_natural_agreement"),
        height_lift=interval("height_lift"),
        compactness_lift=interval("compactness_lift"),
    )

Parameter ranges

lczkit.validation.ranges

Where a computed parameter actually falls, against the range Stewart & Oke published for it.

Agreement says whether a label was right. This says why it was wrong: if a parameter is systematically outside the published band for the class a unit genuinely belongs to, no weighting or threshold downstream can recover the label, because the prototype the unit should match is unreachable in that dimension.

Group by the reference class, not by the assigned one. Grouping by the assigned class asks whether the labelling is self-consistent, which it is almost by construction - the classifier placed those units near that prototype, so their parameter lands near its range. Grouping by the reference class asks whether the parameter reaches the published range on units of known type, which is the question. The two are easy to confuse and give opposite answers: on the Berlin fixture, building surface fraction grouped by assigned class sits inside LCZ 2's published 40-70%, while grouped by reference class it sits at 0.31.

"Reference" names a role, not a file, so the file has to travel with the numbers. Two different references can fill that role - lcz_v3, which is another model's estimate, and hand-labelled So2Sat patches, which are the primary reference wherever they exist - and they do not agree. Substituting one for the other has inverted a diagnosis outright (compactness 25.8% against 55.2% on the same cells), and grouping by lcz_v3 silently uses 91 158 Berlin cells where 9 627 carry labels. reference_file is therefore required to mean anything by a RangeReport: a table that does not name what it grouped by cannot be compared with another one, and the failure is invisible because both look like "the reference class".

Everything is area-weighted, for the reason agreement is: on enclosures, weighting by count lets a thousand street-margin slivers outvote a district.

ClassRange

Bases: BaseModel

One class: where the units it should contain actually sit, against its published range.

published_max instance-attribute

published_max: float | None

The prototype interval in the column's own unit. None is an open end - a class unbounded on that side, which no value can fall outside.

median instance-attribute

median: float | None

Area-weighted median of the parameter over this class' units.

share_in_range instance-attribute

share_in_range: float

Area share of this class' units whose parameter lies inside the published interval. The headline number: 1.0 means the published prototype is reachable for every unit that should match it, 0.0 means it is reachable for none of them.

RangeReport

Bases: BaseModel

One parameter, compared against its published range across every class present.

source instance-attribute

source: str

STEWART_OKE_2012 or lczkit - whether the range being tested against is published at all. Testing a measurement against a range this package invented proves nothing about the measurement, so the provenance travels with the numbers.

grouped_by instance-attribute

grouped_by: str

"reference", "ground_truth" or "assigned". See the module docstring: these answer different questions. "reference" and "ground_truth" are both "the class this unit really is" - they differ in who says so, which is why reference_file is not optional in practice.

reference_file class-attribute instance-attribute

reference_file: str | None = None

The file the grouping labels came from. None only where the caller genuinely has no file to name, as for "assigned", whose labels this package produced itself.

parameter_ranges

parameter_ranges(values: Series, labels: Series, area_m2: Series, *, column: str, grouped_by: str = 'reference', reference_file: str | None = None) -> RangeReport

Compare values against the published range of column, grouped by labels.

All three series are indexed by unit_id. Units with a null value or a null label are excluded - a unit the reference map does not reach cannot tell us anything about whether the parameter reaches a range - and classes absent from labels are omitted rather than reported as empty, since a class no unit belongs to has nothing to say about the parameter.

Pass reference_file wherever the labels came from one: two reports that do not name their reference cannot be told apart, and this module's caller grouped by the wrong one undetected for four phases.

Source code in src/lczkit/validation/ranges.py

def parameter_ranges(
    values: pd.Series,
    labels: pd.Series,
    area_m2: pd.Series,
    *,
    column: str,
    grouped_by: str = "reference",
    reference_file: str | None = None,
) -> RangeReport:
    """Compare `values` against the published range of `column`, grouped by `labels`.

    All three series are indexed by `unit_id`. Units with a null value or a null label are
    excluded - a unit the reference map does not reach cannot tell us anything about whether the
    parameter reaches a range - and classes absent from `labels` are omitted rather than reported
    as empty, since a class no unit belongs to has nothing to say about the parameter.

    Pass `reference_file` wherever the labels came from one: two reports that do not name their
    reference cannot be told apart, and this module's caller grouped by the wrong one undetected
    for four phases.
    """
    spec = property_of(column)
    frame = pd.DataFrame(
        {
            "value": pd.to_numeric(values, errors="coerce"),
            "label": pd.to_numeric(labels, errors="coerce"),
            "area": pd.to_numeric(area_m2, errors="coerce"),
        }
    ).dropna(subset=["value", "label", "area"])

    report = RangeReport(
        column=column,
        property_name=spec.name,
        source=spec.source,
        grouped_by=grouped_by,
        reference_file=reference_file,
        n_units=int(len(frame)),
    )
    for code in sorted(frame["label"].unique()):
        entry = lcz(int(code))
        group = frame[frame["label"] == code]
        lo, hi = ranges_for(entry.code).get(column, (None, None))
        value = group["value"].to_numpy(dtype="float64")
        weight = group["area"].to_numpy(dtype="float64")
        inside = (value >= (-np.inf if lo is None else lo)) & (
            value <= (np.inf if hi is None else hi)
        )
        total = float(weight.sum())
        report.per_class.append(
            ClassRange(
                code=entry.code,
                label=entry.label,
                name=entry.name,
                n=int(len(group)),
                area_m2=total,
                published_min=lo,
                published_max=hi,
                median=weighted_quantile(value, weight, 0.5),
                p10=weighted_quantile(value, weight, 0.1),
                p90=weighted_quantile(value, weight, 0.9),
                share_in_range=float(weight[inside].sum() / total) if total > 0 else 0.0,
            )
        )
    return report

weighted_quantile

weighted_quantile(values: ndarray, weights: ndarray, q: float) -> float | None

The q-quantile of values weighted by weights, or None for an empty input.

The lower of the two interpolation candidates rather than a blend: the quantity is reported beside a published interval, and an interpolated value is not one any unit actually has.

Source code in src/lczkit/validation/ranges.py

def weighted_quantile(values: np.ndarray, weights: np.ndarray, q: float) -> float | None:
    """The `q`-quantile of `values` weighted by `weights`, or `None` for an empty input.

    The lower of the two interpolation candidates rather than a blend: the quantity is reported
    beside a published interval, and an interpolated value is not one any unit actually has.
    """
    if values.size == 0:
        return None
    order = np.argsort(values, kind="stable")
    sorted_values, sorted_weights = values[order], weights[order]
    total = sorted_weights.sum()
    if total <= 0:
        return float(np.median(sorted_values))
    position = int(np.searchsorted(np.cumsum(sorted_weights), q * total))
    return float(sorted_values[min(position, sorted_values.size - 1)])