Data-Driven Classification and Choropleth Break Selection

A choropleth map is only as honest as its class breaks. The same county-level dataset can be rendered to imply a crisis, a gradient, or near-uniformity depending on where the boundaries between colour classes fall — and in an automated pipeline that decision cannot be left to a cartographer nudging sliders. When you generate hundreds of thematic maps from heterogeneous columns, break selection has to be a deterministic function of the data distribution: profile the shape, choose a classification scheme that respects it, compute the breaks reproducibly, and score the result so the pipeline can reject a bad classification before it ships. This guide builds that function on top of mapclassify, the PySAL classification library that also powers geopandas choropleths.

Prerequisites and Environment Configuration

Classification is a numeric operation, so the environment is light — but pin the versions, because mapclassify’s scheme names and return attributes have shifted across minor releases.

  • Python 3.10+ with venv or conda isolation
  • mapclassify>=2.6 — the classification engine; provides FisherJenks, NaturalBreaks, Quantiles, EqualInterval, StdMean, HeadTailBreaks, and the classify convenience wrapper
  • geopandas>=0.14 — reads the vector data and applies the resulting breaks in GeoDataFrame.plot(scheme=...)
  • numpy>=1.24 — distribution profiling (np.percentile, variance decomposition) and edge de-duplication
  • scipy>=1.11 (optional) — scipy.stats.skew and scipy.stats.kurtosis for distribution profiling; numpy alone can approximate these if you want to avoid the dependency
  • Data assumption: classification operates on a single numeric attribute column, not on geometry. Reproject geometry as needed for rendering, but the break computation is CRS-independent — it only reads the value array.

A design note that governs everything below: classification drives both symbology and the legend from one source of truth. The bin edges you compute here are consumed twice — once to assign each feature a fill colour, and once to label the swatches in the legend. If those two consumers ever recompute breaks independently, they will disagree on ties and rounding. Compute the edges once, and pass the same array to your rule-based styling engine and your dynamic legend generator.

Conceptual Foundation: Classification Schemes and the Variance They Preserve

A classification scheme partitions a continuous value array into k ordered classes. Every scheme is an answer to a different question about what the map should emphasise:

  1. Equal interval divides the value range into k equal-width bins: edge_i = min + i·(max − min)/k. It preserves the magnitude of differences, which is ideal for data the reader already understands in absolute units (temperature, elevation). Its fatal weakness is sensitivity to range: a single outlier stretches the interval so the bulk of the data collapses into one class.
  2. Quantile places an equal count of observations in each class. Every colour is guaranteed to appear on roughly the same number of features, which produces a visually balanced map — but it can put nearly identical values in different classes and wildly different values in the same class, exaggerating small differences and hiding large ones.
  3. Natural breaks (Jenks / Fisher-Jenks) minimise within-class variance and maximise between-class variance, finding the “natural” gaps in the distribution. Fisher-Jenks is the exact dynamic-programming optimum; the classic Jenks is an iterative approximation of it. This is the default choice for multimodal or clustered data because it lets the data’s own groupings define the breaks.
  4. Standard deviation places breaks at multiples of the standard deviation around the mean. It only makes sense for roughly symmetric, near-normal distributions, where it communicates “how many sigma from average” — inherently a diverging story that pairs with a diverging colour palette.
  5. Head/tail breaks targets heavy-tailed, power-law data (city sizes, network degrees) by recursively splitting around the mean of the “head.” It produces few classes across the sparse upper tail and fine resolution among the extreme values.

The metric that lets a pipeline compare these objectively is the Goodness of Variance Fit (GVF). It decomposes the sum of squared deviations. Let SDAM be the Sum of squared Deviations from the Array Mean (total variance) and SDCM be the Sum of squared Deviations from each Class Mean (residual within-class variance). Then:

GVF = (SDAM − SDCM) / SDAM = 1 − SDCM / SDAM

GVF ranges from 0 (classes explain nothing) to 1 (each class is internally identical). Fisher-Jenks maximises GVF by construction for a given k; the value it reaches tells you whether the data actually has k separable groups. The diagram below shows the same right-skewed distribution classified three ways, making visible why the scheme choice, not the data, determines the story:

Three Break Placements on One Skewed Distribution A single right-skewed histogram is drawn once. Three rows of vertical break lines are overlaid beneath it: equal-interval breaks land at even value positions leaving the tail bins nearly empty; quantile breaks cluster tightly in the dense left region to equalise counts; natural-breaks fall in the gaps between the data's clusters. min max Right-skewed value distribution (one dataset) Equal interval tail classes nearly empty Quantile breaks crowd the dense head to equalise counts Natural breaks breaks land in the low-density gaps Same data, three stories — the scheme, not the values, decides where the colour boundaries fall

Step-by-Step Implementation

Step 1: Profile the Distribution Before Choosing a Scheme

Never pick a scheme before you have measured the shape. Skewness and kurtosis are cheap and decisive: they tell you whether the mean is a meaningful centre (standard deviation is viable) or a lie told by a heavy tail (it is not).

import numpy as np

def profile_distribution(values: np.ndarray) -> dict:
    """
    Summarise the shape of a value array to drive scheme selection.
    Uses moment-based skewness/kurtosis so scipy is optional.
    """
    x = np.asarray(values, dtype="float64")
    x = x[np.isfinite(x)]                 # drop NaN / inf before profiling
    n = x.size
    mean = x.mean()
    std = x.std()
    # Population skewness and excess kurtosis via central moments
    skew = ((x - mean) ** 3).mean() / std**3 if std > 0 else 0.0
    kurt = ((x - mean) ** 4).mean() / std**4 - 3.0 if std > 0 else 0.0
    return {
        "n": n,
        "unique": int(np.unique(x).size),
        "skew": float(skew),
        "kurtosis": float(kurt),
        "zero_frac": float((x == 0).mean()),
    }

The unique count matters as much as the moments: if a column has fewer distinct values than your target class count, no scheme can produce k non-empty classes, and quantile in particular will emit duplicate edges.

Step 2: Map the Shape to Candidate Schemes

Turn the profile into a shortlist rather than a single answer — the pipeline should compute breaks for two or three candidates and let the GVF score arbitrate. A right-skew above roughly 1.0 rules out standard deviation and equal interval; heavy zero-inflation favours schemes that tolerate ties.

def recommend_scheme(profile: dict) -> list[str]:
    """
    Return an ordered shortlist of mapclassify scheme names given a profile.
    The first entry is the primary recommendation.
    """
    skew = abs(profile["skew"])
    if profile["unique"] < 5:
        return ["Quantiles"]                       # too few values for optimisation
    if profile["zero_frac"] > 0.30 or skew > 2.5:
        return ["HeadTailBreaks", "FisherJenks"]   # heavy tail / power-law
    if skew > 1.0:
        return ["FisherJenks", "Quantiles"]        # moderate right skew
    if skew < 0.5:
        return ["StdMean", "FisherJenks"]          # near-symmetric
    return ["FisherJenks", "Quantiles"]            # default: let the data speak

This mirrors the design logic behind choosing Jenks vs quantile classification for choropleths, which walks through the trade-off between the count-balanced and variance-optimal families in depth.

Step 3: Compute Breaks with mapclassify

mapclassify.classify is the uniform entry point. Every scheme returns an object exposing .bins (the upper edge of each class), .yb (the class index per observation), and .counts (observations per class).

import mapclassify

def compute_breaks(values: np.ndarray, scheme: str, k: int):
    """
    Run a single mapclassify scheme and return the classifier object.
    StdMean ignores k (breaks are fixed at sigma multiples); guard it.
    """
    x = np.asarray(values, dtype="float64")
    x = x[np.isfinite(x)]
    if scheme == "StdMean":
        return mapclassify.classify(x, "StdMean")
    if scheme == "HeadTailBreaks":
        return mapclassify.classify(x, "HeadTailBreaks")
    return mapclassify.classify(x, scheme, k=k)

Step 4: Evaluate GVF and Select the Break Count

The GVF is not exposed uniformly across every scheme, so compute it directly from the class assignments. This makes the score comparable across schemes and independent of mapclassify version quirks.

def gvf(values: np.ndarray, class_index: np.ndarray) -> float:
    """
    Goodness of Variance Fit = 1 - SDCM/SDAM.
    SDAM: squared deviations from the global mean.
    SDCM: squared deviations from each class's own mean.
    """
    x = np.asarray(values, dtype="float64")
    sdam = ((x - x.mean()) ** 2).sum()
    if sdam == 0:
        return 1.0
    sdcm = 0.0
    for c in np.unique(class_index):
        members = x[class_index == c]
        sdcm += ((members - members.mean()) ** 2).sum()
    return 1.0 - sdcm / sdam

On break count: 5 to 7 classes is the working range for sequential choropleths. Below five, real structure is lost; above seven, most readers cannot reliably distinguish adjacent swatches on a sequential ramp, so extra classes add GVF on paper but no communicative value on the map. Let the GVF curve flatten — the “elbow” where adding a class stops materially raising GVF is the honest stopping point.

Complete Working Code Example

The script below ties the pieces into one reusable API: classify_column returns bins plus a GVF score, and auto_classify recommends a scheme from the distribution and picks the best-scoring candidate.

"""
choropleth_classify.py

Data-driven choropleth classification with GVF scoring.
    profile -> recommend -> classify -> score -> select
"""

from __future__ import annotations
from dataclasses import dataclass

import numpy as np
import mapclassify


@dataclass
class Classification:
    scheme: str
    k: int
    bins: list[float]        # upper edge of each class
    counts: list[int]        # observations per class
    gvf: float               # goodness of variance fit, 0..1
    class_index: np.ndarray  # per-observation class assignment


def _clean(values) -> np.ndarray:
    x = np.asarray(values, dtype="float64")
    return x[np.isfinite(x)]


def _gvf(x: np.ndarray, class_index: np.ndarray) -> float:
    sdam = ((x - x.mean()) ** 2).sum()
    if sdam == 0:
        return 1.0
    sdcm = sum(
        ((x[class_index == c] - x[class_index == c].mean()) ** 2).sum()
        for c in np.unique(class_index)
    )
    return float(1.0 - sdcm / sdam)


def classify_column(values, scheme: str = "FisherJenks", k: int = 5) -> Classification:
    """
    Classify one numeric column and report bins plus a GVF score.

    Raises ValueError when the column has fewer distinct values than k,
    which would force empty or duplicate-edge classes.
    """
    x = _clean(values)
    if np.unique(x).size < k:
        raise ValueError(
            f"{np.unique(x).size} unique values cannot fill {k} classes; "
            f"reduce k or choose a tie-tolerant scheme."
        )

    if scheme in ("StdMean", "HeadTailBreaks"):
        clf = mapclassify.classify(x, scheme)           # k is intrinsic
    else:
        clf = mapclassify.classify(x, scheme, k=k)

    bins = [float(b) for b in clf.bins]
    # Collapse duplicate edges (zero-inflated data can produce them)
    dedup = sorted(set(bins))
    if len(dedup) < len(bins):
        bins = dedup

    return Classification(
        scheme=scheme,
        k=len(bins),
        bins=bins,
        counts=[int(c) for c in clf.counts],
        gvf=_gvf(x, clf.yb),
        class_index=clf.yb,
    )


def recommend_scheme(values) -> list[str]:
    """Ordered shortlist of candidate schemes from the distribution shape."""
    x = _clean(values)
    mean, std = x.mean(), x.std()
    skew = abs(((x - mean) ** 3).mean() / std**3) if std > 0 else 0.0
    zero_frac = (x == 0).mean()
    uniq = np.unique(x).size

    if uniq < 5:
        return ["Quantiles"]
    if zero_frac > 0.30 or skew > 2.5:
        return ["HeadTailBreaks", "FisherJenks"]
    if skew > 1.0:
        return ["FisherJenks", "Quantiles"]
    if skew < 0.5:
        return ["StdMean", "FisherJenks"]
    return ["FisherJenks", "Quantiles"]


def auto_classify(values, k: int = 5, min_gvf: float = 0.80) -> Classification:
    """
    Recommend candidates, classify each, and return the highest-GVF result.
    Falls back to Quantiles when no candidate clears min_gvf without empties.
    """
    x = _clean(values)
    best: Classification | None = None
    for scheme in recommend_scheme(x):
        try:
            result = classify_column(x, scheme=scheme, k=k)
        except ValueError:
            continue
        if 0 in result.counts:            # reject empty classes outright
            continue
        if best is None or result.gvf > best.gvf:
            best = result
        if best is not None and best.gvf >= min_gvf:
            break

    if best is None:                      # last-resort, always-valid fallback
        best = classify_column(x, scheme="Quantiles", k=k)
    return best


if __name__ == "__main__":
    rng = np.random.default_rng(42)
    # Right-skewed synthetic column: median income by tract, one wealthy outlier cluster
    income = np.concatenate([
        rng.lognormal(mean=10.8, sigma=0.35, size=480),
        rng.normal(loc=210_000, scale=8_000, size=20),
    ])

    chosen = auto_classify(income, k=6, min_gvf=0.85)
    print(f"scheme={chosen.scheme}  k={chosen.k}  GVF={chosen.gvf:.3f}")
    print(f"bins={[round(b) for b in chosen.bins]}")
    print(f"counts={chosen.counts}")

The synthetic column is deliberately right-skewed with an outlier cluster, so recommend_scheme shortlists Fisher-Jenks ahead of quantile, and auto_classify reports a GVF near the 0.85 threshold with no empty classes — exactly the reproducible, self-validating output a batch pipeline needs.

Performance Optimization Patterns

1. Sample Fisher-Jenks for large n (O(k·n²) → tractable). Exact Fisher-Jenks solves a dynamic program that is quadratic in the number of observations; on arrays past ~10,000 values it becomes the pipeline bottleneck. mapclassify.FisherJenksSampled computes breaks on a random subsample, then applies them to the full array. Fix the seed so the sample — and therefore the breaks — are reproducible across runs.

def breaks_for_scale(values: np.ndarray, k: int = 6, threshold: int = 10_000):
    """Exact Fisher-Jenks under the threshold, sampled above it (seeded)."""
    x = np.asarray(values, dtype="float64")
    x = x[np.isfinite(x)]
    if x.size <= threshold:
        return mapclassify.FisherJenks(x, k=k)
    frac = threshold / x.size
    return mapclassify.FisherJenksSampled(x, k=k, pct=frac, truncate=True)

2. Reuse classifier objects instead of recomputing edges. Once classify_column returns bins, cache the Classification keyed by (column_name, scheme, k). Rendering the symbology and the legend both read from that cached object — this is the mechanism that keeps map and legend in lockstep. Recomputing per consumer is wasted O(k·n²) work and the primary source of map/legend drift.

3. Assign classes with np.digitize, not a Python loop. To apply cached bins to a fresh or streaming array, np.digitize(new_values, bins) is a vectorised O(n log k) binary search — orders of magnitude faster than iterating features, which matters when the same breaks are reapplied to millions of geometries in a tiling job.

4. Profile once, classify many. When several columns share a distribution family (e.g. a panel of yearly income snapshots), compute the profile on the pooled array once and reuse the recommended scheme, rather than re-profiling every column.

Common Pitfalls and Debugging

1. An outlier empties every equal-interval class but one. Symptom: a five-class equal-interval map is 96% one colour. One extreme value stretched the range so four bins sit above the data mass. Fix: profile skewness first (Step 1); for skew above 1.0, recommend_scheme will steer you to Fisher-Jenks or quantile. If equal interval is mandated, clip the outlier into a dedicated top class before classifying.

2. Quantile emits duplicate bin edges on zero-inflated data. Symptom: clf.bins contains the same number twice and the legend shows two identical-range swatches. When more than n/k observations share a value (very common with columns full of zeros), consecutive quantile edges collapse. Fix: de-duplicate edges with sorted(set(bins)) as classify_column does, or drop the number of classes until each is reachable.

3. Non-deterministic Jenks breaks across runs. Symptom: the same input yields different bins on re-run, so cached tiles never match freshly rendered ones. The culprit is a sampled or approximate Jenks with an unseeded RNG, or ties in variance reduction resolved by iteration order. Fix: use exact FisherJenks for reproducibility, or seed the RNG and cache the sample for FisherJenksSampled.

4. Empty classes silently break the legend. Symptom: the legend lists a colour that appears on no feature. A 0 in clf.counts means an empty class — visually it advertises a category that does not exist. Fix: check if 0 in result.counts (as auto_classify does) and reject the candidate, or lower k until every class is populated.

5. Map and legend computed breaks independently. Symptom: the legend says the top class starts at 180,000 but the map fills tracts above 178,500. Two code paths ran classification on slightly different arrays (one included NaNs, one dropped them). Fix: compute the Classification once and pass the identical bins array to both the fill logic and the legend — never re-derive breaks downstream.

Conclusion

Automated choropleth classification is the discipline of turning a subjective cartographic judgement into a reproducible, testable function. Profile the distribution so the mean’s honesty is known, shortlist schemes that respect the shape, compute breaks with mapclassify, and score them with GVF so the pipeline can reject a classification that blurs real structure — all while holding the break count in the 5-to-7 range where colour classes stay perceptually distinct. The single most important invariant is that the bins computed here are the only bins: feed the same array to symbology and legend, and the two can never disagree. From here, the natural next steps are wiring these breaks into a rule-based styling engine that maps each class to a fill, pairing them with a perceptually uniform ramp from color theory for GIS, and returning to the parent programmatic map styling and label automation workflow that consumes classified output end to end.

FAQ

Why does Fisher-Jenks return different breaks on identical data across runs? Fisher-Jenks itself is deterministic, but the older jenkspy and the sampled approximation are not. When two candidate break edges yield the same variance reduction, the tie is resolved by iteration order, and sampling a large array changes which values are seen. Use mapclassify.FisherJenks on the full array for reproducibility, or fix the numpy seed and cache the sample when using FisherJenksSampled.

How do I detect empty or duplicate-edge classes after classification? Inspect the counts attribute of the mapclassify result; any zero entry is an empty class. Duplicate bin edges appear when more than k values are identical (common with zero-inflated data), collapsing two breaks to the same number. Both conditions produce a legend with unreachable colours. Drop duplicate edges with numpy.unique or switch from quantile to a scheme that tolerates ties.

What GVF value indicates an acceptable classification? The Goodness of Variance Fit ranges from 0 to 1, where 1 means classes capture all variance. In practice a GVF above 0.80 at five to seven classes indicates the breaks separate natural groupings well. Below 0.70 the classes blur real clusters; chasing values above 0.95 usually just means adding classes past the point of perceptual usefulness.

Why do equal-interval classes look empty on skewed economic data? Equal interval divides the value range into equal-width bins regardless of density. A single outlier (one very high income tract) stretches the range so that four of five bins fall above the bulk of the data, leaving them empty while every observation crowds into the first class. Use quantile or Fisher-Jenks for right-skewed data, or clip the outlier into a dedicated top class.

Back to Programmatic Map Styling and Label Automation