Choropleth Classification Explained: Quantiles, Equal Interval and Natural Breaks #

Problem statement #

The same data, two lines of code, two completely different maps:

gdf.plot(column="income", scheme="quantiles", k=5, cmap="YlOrRd", legend=True)

gdf.plot(column="income", scheme="equal_interval", k=5, cmap="YlOrRd", legend=True)

The first shows a patchwork with roughly equal numbers of wards in each colour. The second is almost entirely pale yellow with three dark red wards. Nothing about the data changed. Neither map is wrong. They are answers to different questions, and a reader shown one of them will not know the other exists.

Worse is the default:

gdf.plot(column="income", cmap="YlOrRd", legend=True)

No classification at all — a continuous colour ramp stretched linearly from minimum to maximum. One outlier compresses everything else into a narrow band of indistinguishable colour, and the map shows nothing but where the outlier is.

Classification is the step that turns numbers into colours, and it is where a choropleth acquires most of its editorial content.

Quick answer #

import geopandas as gpd
import mapclassify

gdf.plot(column="income", scheme="quantiles", k=5,
         cmap="YlOrRd", legend=True, edgecolor="white", linewidth=0.3)

Always pass a scheme. The default continuous ramp is the one option that reliably hides the data.

# rates and ratios, never counts
gdf["rate"] = gdf["cases"] / gdf["population"] * 100_000
gdf.plot(column="rate", scheme="quantiles", k=5, legend=True)

Step-by-step solution #

1. Map a rate, not a count #

Before choosing any scheme: a choropleth colours areas, and area is proportional to size, not to population. Colouring a raw count means the map shows where people are, and only incidentally what you meant.

# ❌ this map shows population distribution, whatever the column is called
gdf.plot(column="covid_cases", scheme="quantiles", k=5)

# ✅ a rate is comparable between a city ward and a rural district
gdf["cases_per_100k"] = gdf["covid_cases"] / gdf["population"] * 100_000
gdf.plot(column="cases_per_100k", scheme="quantiles", k=5)

Anything divided by an appropriate denominator works: per capita, per square kilometre, per household, as a percentage. This single change fixes more misleading choropleths than every classification decision combined.

2. Look at the distribution before choosing #

The right scheme follows from the shape of the data:

import matplotlib.pyplot as plt
import numpy as np

def profile(gdf, column):
    v = gdf[column].dropna()
    print(f"{column}: n={len(v):,}  missing={gdf[column].isna().sum():,}")
    print(f"  min {v.min():,.1f}  median {v.median():,.1f}  max {v.max():,.1f}")
    print(f"  mean {v.mean():,.1f}  skew {v.skew():+.2f}")
    for q in [0.01, 0.25, 0.5, 0.75, 0.95, 0.99]:
        print(f"  p{int(q*100):<3} {v.quantile(q):>12,.1f}")
    fig, ax = plt.subplots(figsize=(8, 3))
    ax.hist(v, bins=60, color="#0ea5e9")
    ax.set_title(f"{column} — skew {v.skew():+.2f}")
    return fig

profile(gdf, "income")

income: n=8,436  missing=112
  min 8,204.0  median 31,880.0  max 412,004.0
  mean 36,118.4  skew +4.82
  p1        12,204.0
  p25       24,880.0
  p50       31,880.0
  p75       41,220.0
  p95       78,404.0
  p99      184,002.0

A skew of +4.8 and a p99 more than twice p95 says this data is heavily right-tailed. equal_interval on it will put 99% of wards in the first class. quantiles will spread them out. That is the decision made, in one histogram.

3. Understand what each scheme does #

import mapclassify as mc

v = gdf["income"].dropna()
for name, cls in [("quantiles", mc.Quantiles),
                  ("equal_interval", mc.EqualInterval),
                  ("natural_breaks", mc.NaturalBreaks),
                  ("std_mean", mc.StdMean)]:
    c = cls(v, k=5) if name != "std_mean" else cls(v)
    print(f"{name:<16} breaks {[f'{b:,.0f}' for b in c.bins]}")
    print(f"{'':<16} counts {list(c.counts)}")

quantiles        breaks ['24,880', '29,104', '35,220', '46,880', '412,004']
                 counts [1688, 1687, 1687, 1687, 1687]
equal_interval   breaks ['88,964', '169,724', '250,484', '331,244', '412,004']
                 counts [8218, 174, 32, 8, 4]
natural_breaks   breaks ['38,204', '68,880', '124,004', '228,880', '412,004']
                 counts [6412, 1604, 336, 68, 16]
std_mean         breaks ['-31,204', '4,457', '67,780', '103,441', '412,004']
                 counts [0, 412, 7658, 344, 22]

Read the counts, not the breaks. equal_interval puts 8,218 of 8,436 wards in one class — the map will be one colour. quantiles splits them evenly by construction. natural_breaks finds a middle position that follows the data's own structure.

Quantiles guarantee equal counts, so every class is visible and the map is about rank. The cost: two wards on either side of a break can be nearly identical in value, and a class can span a huge range at the tail. A quantile map always looks varied, even when the underlying differences are trivial.

Equal interval divides the value range into equal parts, so the legend is intuitive and comparable across maps. It only works when the data is reasonably evenly spread. On skewed data it produces a map of one colour.

Natural breaks (Jenks) minimises within-class variance, putting boundaries where the data has genuine gaps. It follows the data's structure, and it is data-dependent — two maps of different years get different breaks and cannot be compared.

Standard deviations classify by distance from the mean, which is the right frame for a diverging colour scheme showing above and below average. It assumes something roughly normal; on our skewed income data the first class is empty.

4. Choose the number of classes #

for k in [3, 5, 7, 9]:
    c = mc.Quantiles(v, k=k)
    print(f"k={k}  ADCM {c.adcm:>12,.0f}")

Five to seven is the usual range. Below four the map loses detail; above about eight, readers cannot reliably match a polygon's colour to a legend swatch — a limit of perception rather than of the data.

mapclassify can compare schemes numerically:

print(mc.classify(v, "quantiles", k=5).adcm)        # absolute deviation about class medians

Lower ADCM means classes that group similar values more tightly. It is a useful tiebreaker between schemes at the same k, and it is not a substitute for the question of what the map is for — a low-ADCM scheme nobody can compare across years may still be the wrong choice.

5. Match the colour scheme to the data type #

Classification decides the breaks; the colour ramp decides whether they read correctly.

# diverging data must have its midpoint at the meaningful zero
import numpy as np
gdf["change_pct"] = (gdf["pop_2026"] - gdf["pop_2016"]) / gdf["pop_2016"] * 100
lim = np.abs(gdf["change_pct"]).quantile(0.98)
gdf.plot(column="change_pct", cmap="RdBu", vmin=-lim, vmax=lim, legend=True)

The symmetric vmin/vmax is what makes a diverging map honest. Without it, matplotlib centres the ramp on the midpoint of the data range rather than on zero, so "no change" gets a colour from one side and growth and decline are not visually comparable.

Code examples #

Example 1: comparing every scheme at once #

Before publishing, look at the alternatives:

import matplotlib.pyplot as plt
import geopandas as gpd
import mapclassify as mc

SCHEMES = ["quantiles", "equal_interval", "natural_breaks", "std_mean",
           "fisher_jenks", "boxplot"]

def compare_schemes(gdf, column, k=5, cmap="YlOrRd"):
    fig, axes = plt.subplots(2, 3, figsize=(18, 11))
    v = gdf[column].dropna()
    for ax, scheme in zip(axes.ravel(), SCHEMES):
        try:
            gdf.plot(column=column, scheme=scheme, k=k, cmap=cmap, ax=ax,
                     edgecolor="white", linewidth=0.2,
                     missing_kwds={"color": "#e2e8f0", "hatch": "///"})
            classifier = mc.classify(v, scheme, k=k)
            counts = list(classifier.counts)
            biggest = 100 * max(counts) / sum(counts)
            ax.set_title(f"{scheme}\ncounts {counts} · largest class {biggest:.0f}%",
                         fontsize=9)
        except Exception as exc:
            ax.set_title(f"{scheme}: {exc}", fontsize=8)
        ax.set_axis_off()
    plt.tight_layout()
    return fig

compare_schemes(gpd.read_file("wards.gpkg").to_crs(27700), "income")

The "largest class %" in each title is the number to scan. Anything above about 60% means most of the map is one colour and the scheme is not separating the data. On the income example, equal_interval reports 97% and can be dismissed immediately.

missing_kwds is worth keeping as a habit: without it, features with a null value are simply not drawn, leaving holes a reader will interpret as zero rather than unknown. See how to handle missing and null values.

Example 2: fixed breaks so maps compare across years #

Data-driven schemes recompute their breaks per dataset, so a 2016 map and a 2026 map classified with natural_breaks are not comparable — the same colour means different values.

import geopandas as gpd
import mapclassify as mc
import matplotlib.pyplot as plt

years = [2016, 2021, 2026]
layers = {y: gpd.read_file(f"wards_{y}.gpkg") for y in years}

# derive one set of breaks from the pooled data, then reuse it
pooled = gpd.pd.concat([g[["income"]] for g in layers.values()])["income"].dropna()
breaks = list(mc.Quantiles(pooled, k=5).bins)
print("shared breaks:", [f"{b:,.0f}" for b in breaks])

fig, axes = plt.subplots(1, len(years), figsize=(6 * len(years), 6))
for ax, y in zip(axes, years):
    layers[y].plot(column="income", scheme="user_defined",
                   classification_kwds={"bins": breaks},
                   cmap="YlOrRd", ax=ax, edgecolor="white", linewidth=0.2,
                   legend=(y == years[-1]))
    ax.set_title(str(y))
    ax.set_axis_off()
plt.tight_layout()

shared breaks: ['24,204', '29,880', '36,112', '48,004', '412,004']

scheme="user_defined" with explicit bins is what makes a series of maps a series rather than three unrelated pictures. Deriving those bins from the pooled data means the breaks suit the whole period rather than any one year.

A legend on only the last panel avoids repeating the same key three times.

Example 3: an honest choropleth, end to end #

import geopandas as gpd
import matplotlib.pyplot as plt
import mapclassify as mc
import numpy as np

def choropleth(gdf, column, *, denominator=None, per=100_000, scheme="quantiles",
               k=5, cmap="YlOrRd", crs=None, title=None, note=None):
    gdf = gdf.copy()

    if denominator:
        col = f"{column}_per_{per:,}"
        gdf[col] = gdf[column] / gdf[denominator].replace(0, np.nan) * per
    else:
        col = column

    if crs:
        gdf = gdf.to_crs(crs)
    elif gdf.crs and gdf.crs.is_geographic:
        raise ValueError("geographic CRS — pass an equal-area crs= for a choropleth")

    missing = gdf[col].isna().sum()
    fig, ax = plt.subplots(figsize=(10, 11))
    gdf.plot(column=col, scheme=scheme, k=k, cmap=cmap, ax=ax,
             edgecolor="white", linewidth=0.25,
             legend=True, legend_kwds={"loc": "lower right", "fontsize": 9,
                                       "title": col.replace("_", " ")},
             missing_kwds={"color": "#e2e8f0", "hatch": "///",
                           "label": f"no data ({missing})"})

    classifier = mc.classify(gdf[col].dropna(), scheme, k=k)
    counts = list(classifier.counts)
    ax.set_title(title or col, fontsize=14, loc="left")
    ax.set_axis_off()
    ax.annotate(
        (note or "") + f"\n{scheme}, k={k} · class counts {counts}"
        f" · {len(gdf) - missing:,} of {len(gdf):,} areas mapped",
        xy=(0.01, 0.01), xycoords="axes fraction", fontsize=7.5, color="#64748b")
    return fig, classifier

fig, cls = choropleth(
    gpd.read_file("wards.gpkg"),
    "covid_cases", denominator="population", per=100_000,
    scheme="quantiles", k=5, crs=27700,
    title="COVID-19 cases per 100,000 residents",
    note="Source: UKHSA, 2026-08-01. Rate, not count.",
)

Four things this function refuses to let you get wrong. It computes a rate when given a denominator, rather than trusting you to remember. It raises on a geographic CRS instead of silently drawing a distorted map — the reasoning is in choosing a map projection for display. It draws missing values in a distinct hatched grey with a count, so absent data cannot be mistaken for low data. And it prints the scheme, k and class counts on the figure itself, so a reader can see how the colours were assigned rather than having to trust them.

The replace(0, np.nan) on the denominator turns a division by zero into a missing value, which is honest, rather than into an infinity, which plots as the maximum class.

Explanation #

A choropleth performs two mappings: from a value to a class, and from a class to a colour. The second is a design decision. The first is an analytical one, and it determines what the map is capable of showing.

The reason schemes differ so much is that they optimise different things. Quantiles optimise for visual balance — equal counts per class guarantee every colour appears, which makes a map that looks informative regardless of whether the underlying differences are meaningful. Equal interval optimises for legend legibility — even, round-ish ranges that a reader can hold in their head and compare with another map. Natural breaks optimise for within-class homogeneity, placing boundaries where the data has genuine gaps.

These goals conflict. On skewed data — which is most real spatial data, since income, population density, house prices and disease counts are all right-tailed — the conflict is stark. Equal interval collapses to one colour. Quantiles produce a balanced map whose top class spans from £48,000 to £412,000, so wildly different places share a shade. Neither is dishonest, and both need to be labelled.

This is why the class counts belong on the map, and why a published choropleth should state its scheme. A reader looking at a quantile map is looking at ranks; a reader looking at an equal-interval map is looking at magnitudes. Without the label they cannot tell which, and the two support different conclusions from identical data.

The count-versus-rate error dominates everything else. A choropleth encodes value as colour over an area, and the reader integrates colour across area whether or not you intend it. So a raw count mapped this way produces a picture of where the denominator is large — usually population. This is not a subtle bias; it is the map showing a different variable from the one in the legend. Dividing by an appropriate denominator is the single highest-value correction available.

And the classification interacts with the projection. Both are ways the map's visual weight can diverge from the data's actual weight: an unequal-area projection makes northern regions look larger, and a poor classification makes them look more extreme. Together they compound. An equal-area projection and a scheme suited to the distribution are the two halves of a choropleth that says what its legend claims.

Finally, note the limit of the form. A choropleth cannot show within-area variation, and it inherits whatever the areal units impose — the modifiable areal unit problem means the same underlying data aggregated to different boundaries can produce opposite patterns. When that matters, a dot density map, a cartogram, or mapping the underlying points is more honest than a better classification of the aggregate.

Edge cases or notes #

• scheme= requires mapclassify. pip install mapclassify, or GeoPandas raises when you pass one.
• Ties break quantiles. If many features share a value, quantile classes cannot have equal counts and mapclassify warns.
• missing_kwds is essential. Without it, null-valued features are not drawn at all and read as zero.
• k above about 8 exceeds what readers can match to a legend, however much detail the data has.
• Jenks is slow on large datasets — natural_breaks samples above 1,000 features; fisher_jenks is exact and O(n²).
• Diverging ramps need symmetric vmin/vmax, or the neutral colour lands somewhere other than zero.
• std_mean assumes an approximately normal distribution. On skewed data its first class is often empty.
• user_defined bins are upper bounds, and the last one must be at least the data maximum.
• Colour-blind readers cannot separate red from green; viridis, cividis and ColorBrewer's colour-blind-safe ramps can.
• legend_kwds={"fmt": "{:,.0f}"} stops a legend showing eight decimal places.

Internal links #

• How to create a choropleth map in Python with GeoPandas — the practical build
• Choosing a map projection for display — the other half of an honest choropleth
• My GeoPandas choropleth colours look wrong — when the map does not match the data
• How to plot maps in Python with GeoPandas and Matplotlib — plotting fundamentals
• How to handle missing and null values in spatial datasets — what missing_kwds is for
• How to aggregate spatial data by region in GeoPandas — producing the values you classify
• How to add legends, labels and a scale bar to a GeoPandas map — making the classes readable
• Spatial data quality: the six dimensions that matter — what the numbers are worth before you map them

FAQ #

Which classification scheme should I use? #

quantiles for showing rank and guaranteeing every class is visible; equal_interval when the data is evenly spread and the legend must be intuitive; natural_breaks to find real groupings. Always pass one — the default continuous ramp hides everything behind outliers.

Why is my whole map one colour? #

Almost certainly equal_interval on skewed data, where the top few values stretch the range so far that nearly every feature falls in the first class. Check the class counts and switch to quantiles or natural_breaks.

Should I map counts or rates? #

Rates, essentially always. A choropleth colours areas, so a raw count produces a map of where the denominator is large — usually population — regardless of what the legend says.

How many classes? #

Five to seven. Fewer loses detail; more exceeds what a reader can reliably match between a polygon and a legend swatch.

Why do my maps for different years look inconsistent? #

Data-driven schemes recompute breaks per dataset. Derive one set of breaks from the pooled data and pass them with scheme="user_defined".

What is Jenks natural breaks actually doing? #

Minimising variance within classes and maximising it between them, so boundaries land where the data has genuine gaps. It is data-dependent, which is why it does not compare across datasets.

How should I show missing data? #

With missing_kwds — a distinct grey, ideally hatched, and a count in the legend. Undrawn polygons read as zero, which is a different claim from unknown.