Choosing Jenks vs Quantile Classification for Choropleths

Choose Jenks natural breaks when the map must reveal the data’s own groupings and choose quantile classification when the map must support ranked comparison across regions: Jenks minimises within-class variance so its breaks track natural gaps in the distribution, but it is data-specific and runs in O(k·n^2), whereas quantile puts an equal number of features in every class so it is stable and reproducible across datasets — at the cost of splitting near-identical values and merging genuinely distinct ones. The rest of this page turns that trade-off into a distribution-and-purpose decision rule, a goodness-of-variance-fit comparison, class-count effects, and a reproducible chooser built on mapclassify.

Core Algorithm and Workflow

Both schemes assign every polygon a class index that maps to a colour, but they optimise for opposite things. Jenks (the Fisher-Jenks natural-breaks algorithm) searches for the set of k-1 break edges that minimises the sum of squared deviations of values from their class means — equivalently, it maximises the goodness of variance fit (GVF). Quantile ignores value spacing entirely and slices the sorted array into k equal-count segments, so break edges fall at the 1/k, 2/k, … percentile positions.

The consequence is visible in the diagram below: the same fifteen sample values, classified into four classes, produce different edges and different per-class counts under each scheme.

Read the two rows against the shared ticks. Jenks drops its edges into the empty gaps, so tightly grouped observations stay together and the sparse high tail becomes its own class of just a few features. Quantile forces four-per-class, so one of its edges lands inside the mid group, separating two observations that are almost equal in value — the classic quantile distortion. Conversely, Jenks assigns only two features to class 3, which can read as a fragile, cartographically noisy class if those two polygons are tiny.

That gives a decision rule keyed to distribution shape and map purpose:

1. Skewed or multimodal data, map reads absolute magnitude. Prefer Jenks. Population density, income, and pollutant concentration are heavy-tailed; natural breaks keep the dense low end from being over-partitioned and give the outliers their own honest class.
2. Near-uniform data, or the map must support ranked comparison. Prefer quantile. When readers ask “which regions are in the top fifth?”, equal-count classes answer directly and stay comparable if you re-run the map on next year’s data.
3. Many tied values (rates with lots of zeros, capped indices). Neither is automatic — quantile will collapse bins onto duplicate edges, and Jenks may strand a class. Reduce k, or move to a scheme designed for the shape.
4. A series of maps that must share a legend. Avoid Jenks: its breaks are recomputed per dataset, so the legend shifts between maps. Fix explicit breaks or use quantile with a pinned reference distribution.

This decision fits inside the broader Data-Driven Classification workflow, where the chosen break edges are handed to the colour ramp and the renderer.

Production-Ready Python Implementation

The script below is a complete, copy-pasteable chooser. It uses mapclassify==2.6.1 (the classification engine behind GeoPandas’ scheme= argument) and numpy==1.26.4. It computes GVF for both schemes at a fixed k, detects the two failure modes — quantile duplicate edges and Jenks undersized classes — and returns a recommendation with the break edges ready to bind to a colour ramp.

import numpy as np
import mapclassify
import geopandas as gpd


def goodness_of_variance_fit(values: np.ndarray, classifier) -> float:
    """
    Compute the Goodness of Variance Fit (GVF) for a fitted classifier.

    GVF = (SDAM - SDCM) / SDAM, where
      SDAM = sum of squared deviations from the array mean, and
      SDCM = sum of squared deviations from each class mean.
    Ranges from 0 (no fit) to 1 (perfect fit). Jenks maximises it by design.
    """
    yb = classifier.yb  # class index for each observation
    sdam = np.sum((values - values.mean()) ** 2)
    sdcm = 0.0
    for cls in np.unique(yb):
        members = values[yb == cls]
        sdcm += np.sum((members - members.mean()) ** 2)
    return float((sdam - sdcm) / sdam) if sdam > 0 else 0.0


def choose_scheme(
    values: np.ndarray,
    k: int = 5,
    min_class_share: float = 0.03,
    purpose: str = "magnitude",
) -> dict:
    """
    Recommend Jenks natural breaks vs quantile classification for a choropleth.

    Parameters
    ----------
    values          : 1-D array of the choropleth variable (NaNs removed upstream).
    k               : Number of classes. 4-7 keeps the legend readable.
    min_class_share : Minimum fraction of features a class may hold before it is
                      flagged as fragile (Jenks tail classes trip this).
    purpose         : "magnitude" (absolute reading) or "ranked" (comparison).

    Returns
    -------
    dict with per-scheme GVF, break edges, detected warnings, and a recommendation.
    """
    values = np.asarray(values, dtype=float)
    n = values.size

    jenks = mapclassify.NaturalBreaks(values, k=k)
    quant = mapclassify.Quantiles(values, k=k)

    gvf_jenks = goodness_of_variance_fit(values, jenks)
    gvf_quant = goodness_of_variance_fit(values, quant)

    warnings: list[str] = []

    # Failure mode 1: quantile collapses onto duplicate edges when values are tied.
    q_edges = np.round(quant.bins, 10)
    if len(np.unique(q_edges)) < len(q_edges):
        warnings.append(
            "quantile produced duplicate break edges (tied values); "
            "effective classes < k"
        )

    # Failure mode 2: Jenks strands a class with too few features to read cleanly.
    counts = np.bincount(jenks.yb, minlength=k)
    fragile = [int(c) for c in counts if c / n < min_class_share]
    if fragile:
        warnings.append(
            f"jenks has fragile class(es) below {min_class_share:.0%} "
            f"of features: sizes {fragile}"
        )

    # Decision rule: purpose dominates; GVF is a tie-breaker within a purpose.
    skew = float(
        np.mean(((values - values.mean()) / (values.std() + 1e-12)) ** 3)
    )
    if purpose == "ranked":
        recommended = "quantile"
    elif abs(skew) >= 1.0 and not fragile:
        recommended = "jenks"
    else:
        # roughly symmetric, or Jenks is fragile: prefer whichever fits better
        recommended = "jenks" if gvf_jenks - gvf_quant > 0.05 else "quantile"

    return {
        "n": n,
        "k": k,
        "skew": round(skew, 3),
        "gvf": {"jenks": round(gvf_jenks, 4), "quantile": round(gvf_quant, 4)},
        "bins": {
            "jenks": [round(float(b), 4) for b in jenks.bins],
            "quantile": [round(float(b), 4) for b in quant.bins],
        },
        "class_counts": {
            "jenks": counts.tolist(),
            "quantile": np.bincount(quant.yb, minlength=k).tolist(),
        },
        "warnings": warnings,
        "recommended": recommended,
    }


# ---------------------------------------------------------------------------
# Example usage
# ---------------------------------------------------------------------------
# gdf = gpd.read_file("counties_income.gpkg")
# col = "median_income"
# vals = gdf[col].dropna().to_numpy()
#
# report = choose_scheme(vals, k=5, purpose="magnitude")
# scheme = report["recommended"]
# breaks = report["bins"][scheme]
# print(scheme, "GVF", report["gvf"][scheme], "breaks", breaks)
#
# # Bind the chosen breaks directly to the renderer's colour ramp:
# gdf.plot(column=col, scheme="user_defined",
#          classification_kwds={"bins": breaks}, cmap="viridis", legend=True)

NaturalBreaks in mapclassify is the Fisher-Jenks optimiser; on large arrays it samples down before optimising, so for reproducible edges across runs pin the array or pass a seed via mapclassify.FisherJenksSampled. See the mapclassify classifiers reference for the full scheme list and their keyword arguments.

Performance Tuning and Cartographic Best Practices

• Mind the complexity gap. Exact Fisher-Jenks is O(k·n^2) in both time and memory because it fills a dynamic-programming matrix over every value pair. Quantiles is O(n log n) — just a sort and index. For n above ~10,000 features, switch from NaturalBreaks to FisherJenksSampled, which optimises on a random subset (default 10%) and applies the resulting breaks to the full array, or accept quantile if its distortion is tolerable for the map’s purpose.
• Fix k before comparing, not after. GVF rises monotonically with k for any scheme, so comparing a 7-class Jenks against a 4-class quantile is meaningless. Hold k between 4 and 7 — the readable range for a printed legend — and only then read the GVF delta. A Jenks-over-quantile GVF gain under ~0.05 rarely justifies the loss of comparability.
• Treat GVF as a tie-breaker, not a verdict. Jenks maximises GVF by construction, so it will almost always win the statistic. That does not make it the right map. Let purpose (absolute magnitude vs ranked comparison) decide first; use GVF only to separate two schemes that both suit the purpose.
• Guard against tied values before choosing quantile. Rates with many zeros, capped indices, and rounded counts create ties that collapse quantile edges onto duplicates, silently giving you fewer classes than requested. The chooser above flags this; when it fires, drop k or move to natural breaks.
• Pin explicit breaks for map series and time animation. Because Jenks recomputes per dataset, an animated or multi-panel choropleth will flicker as the legend shifts. Compute breaks once on a reference year (or a pooled distribution), then pass them as user_defined bins so every frame shares one legend — the same discipline that keeps a colour ramp stable across a series, discussed under Color Palette Generation for Thematic Maps.

Integration and Next Steps

The chooser returns break edges as a plain numeric list, which keeps it renderer-agnostic and easy to wire into a styling pipeline:

• GeoPandas static export. Pass the chosen list straight through classification_kwds={"bins": breaks} with scheme="user_defined", as in the example, so the map and any offline preview share identical edges.
• Legend synchronisation. Feed the same breaks and the class counts into the legend builder so the swatch labels read as real value ranges rather than raw indices — the equal-count property of quantile in particular needs explicit range labels to be honest, which pairs with the approach in Dynamic Legend Generation.
• Rule-based styling engines. Emit the edges as a JSON list of stop values and bind them to a colour ramp inside a declarative style so classification and symbology stay in one artefact, following the pattern in Rule-Based Styling Engines.
• CI validation. Run choose_scheme in a test that fails the build when warnings is non-empty for a published dataset, catching duplicate-edge and fragile-class regressions before a distorted choropleth ships.

For the next step, apply the chosen classification to a perceptually ordered ramp so that class order reads correctly, and validate the ramp’s contrast before publishing the map.

Frequently Asked Questions

Why do two choropleths of the same data look completely different under Jenks and quantile?

Jenks places break edges where the data has natural gaps, so class widths vary and the darkest class may hold only a handful of features. Quantile forces equal feature counts per class, so it spreads colour evenly across the map regardless of the underlying values. On a skewed variable the same region can land in the top class under one scheme and the middle under the other — which is exactly why the two schemes are not interchangeable and must be chosen against the map’s purpose.

Is a higher GVF always the better choice?

No. GVF rewards minimising within-class variance, so Jenks will almost always score higher than quantile at the same k. That makes GVF a valid tie-breaker only when both schemes already suit the map’s purpose. If the map is meant for ranked comparison across regions, a lower-GVF quantile scheme can still be correct because it keeps class sizes comparable across the whole map.

Why does mapclassify warn about duplicate or non-unique bins with quantile classification?

Quantile edges are computed from percentile positions. When many observations share the same value — a rate that is zero across half the units, for instance — several percentile edges collapse onto the same number, and you end up with fewer usable classes than requested. Switch to NaturalBreaks or reduce k until the bins are unique; a legend with repeated edges conveys no information.

Related

• Data-Driven Classification — the parent page covering how classification schemes drive choropleth colour, from equal-interval to head/tail breaks.
• Color Palette Generation for Thematic Maps — build the perceptually ordered ramp that the chosen break edges bind to.
• Dynamic Legend Generation — turn the chosen breaks and class counts into honest, readable value-range legends.

Back to Programmatic Map Styling and Label Automation