My GeoPandas Choropleth Colours Look Wrong: How to Fix It

On this page

Problem statement

The map renders and the colours make no sense.

The whole thing is one shade of pale yellow with a few dark polygons. Or it is a rainbow patchwork with no visible pattern. Or the highest values are pale and the lowest dark, exactly backwards. Or whole regions are missing and the holes look like zeros. Or the "no change" colour on a diverging map sits somewhere other than zero, so growth and decline are not comparable.

Here is every US county coloured by deaths in 2020, from the Census Bureau's own population estimates (the setup is under The data used on this page):

import matplotlib
matplotlib.use("Agg")

counties.plot(column="deaths", cmap="YlOrRd", legend=True)

No error, no warning β€” and almost the whole country comes out one pale colour. The map is a faithful rendering of a set of decisions you did not know you were making, mostly about classification, outliers, missing values and colour direction.

There are seven distinct causes. Each has a specific tell, and most take one line to fix.

Quick answer

Look at the data before blaming the plot:

v = counties["deaths"]
print(f"n={len(v):,}  missing={v.isna().sum():,}  dtype={v.dtype}  skew={v.skew():+.2f}")
print(v.describe(percentiles=[.01, .25, .5, .75, .95, .99]).round(1).to_string())
n=3,234  missing=91  dtype=float64  skew=+11.45
count     3143.0
mean       977.0
std       2680.8
min          0.0
1%           7.0
25%        130.0
50%        303.0
75%        755.0
95%       4203.5
99%      11911.9
max      73513.0
Symptoms of a wrong-looking choropleth matched to their causes. the distribution names most of theseWhat the map showsWhat is actually wrongone colour, two dark polygonsan outlier is flattening the rampeverything in the first classequal_interval on skewed datacolours in a random orderthe column is a string, not anumberholes where polygons should benulls are not drawn β€” passmissing_kwdsthe pattern is just the citiesa count is being mapped, not arateOnly two of these are about colour at all spatialworkflow.io
Seven causes. The distribution and the dtype identify most of them.
What you see Cause Fix
one colour, a few extremes no scheme, or an outlier scheme="quantiles", or clip
everything in class 1 equal_interval on skewed data scheme="quantiles"
colours look random column is a string, not a number pd.to_numeric(...)
high values look pale reversed colormap cmap="YlOrRd" not "YlOrRd_r"
holes in the map nulls not drawn missing_kwds={...}
neutral colour is not at zero asymmetric diverging range vmin=-lim, vmax=lim
pattern follows population mapping a count, not a rate divide by a denominator
counties.plot(column="death_rate", scheme="quantiles", k=5, cmap="YlOrRd",
              legend=True, legend_kwds={"fmt": "{:,.1f}"},
              missing_kwds={"color": "#e2e8f0", "hatch": "///", "label": "no data"})

Outliers, the wrong scheme, text columns, reversed ramps, holes

1. One colour with a couple of extremes β€” an outlier is eating the ramp

The default is a continuous linear stretch from minimum to maximum. One extreme value compresses everything else into a sliver of the ramp:

print(f"p95 {v.quantile(0.95):,.0f}   max {v.max():,.0f}   ratio {v.max() / v.quantile(0.95):.1f}Γ—")
print(f"share of the colour range used by 95% of counties: {v.quantile(0.95) / v.max():.1%}")
print("the maximum:", counties.loc[v.idxmax(), "NAMELSAD"] + ", " + counties.loc[v.idxmax(), "STATE_NAME"])
p95 4,203   max 73,513   ratio 17.5Γ—
share of the colour range used by 95% of counties: 5.7%
the maximum: Los Angeles County, California

One county β€” the most populous in the country β€” pushes the top of the ramp so far out that 95% of counties share a sliver of the colours. Three fixes, in order of preference:

# (a) classify β€” the real answer
counties.plot(column="deaths", scheme="quantiles", k=5, cmap="YlOrRd", legend=True)

# (b) clip the continuous range at percentiles
lo, hi = v.quantile([0.02, 0.98])
counties.plot(column="deaths", cmap="YlOrRd", vmin=lo, vmax=hi, legend=True)

# (c) log scale, for genuinely multiplicative data
from matplotlib.colors import LogNorm
counties.plot(column="deaths", cmap="YlOrRd", norm=LogNorm(v[v > 0].min(), v.max()), legend=True)

Classification is the durable fix, because it makes the breaks explicit and puts them in the legend. Clipping hides the extremes rather than accounting for them, so say so in the caption when you use it.

2. Everything in the first class β€” the wrong scheme for the distribution

import mapclassify as mc

for scheme in ["equal_interval", "quantiles", "fisher_jenks"]:
    counts = [int(n) for n in mc.classify(v.dropna(), scheme, k=5).counts]
    print(f"{scheme:<15} counts {str(counts):34} largest class {100 * max(counts) / sum(counts):.0f}%")
equal_interval  counts [3123, 16, 2, 1, 1]                largest class 99%
quantiles       counts [630, 627, 633, 624, 629]          largest class 20%
fisher_jenks    counts [2752, 285, 82, 22, 2]             largest class 88%

equal_interval divides the value range evenly, which on right-skewed data β€” deaths, density, prices, case counts β€” puts nearly everything in the first class. A largest-class share above about 60% means the scheme is not separating your data; here even Fisher-Jenks, which follows the data, leaves most counties in one class. Full comparison in choropleth classification explained.

US counties per equal-interval class of deaths in 2020, against quantiles. equal interval on 3,143 counties: counties per classclass 1 (to 14,703)3,123 β€” 99%class 216class 32class 41class 5 (Los Angeles)103123 countiesQuantiles on the same numbers: 630, 627, 633, 624, 629 spatialworkflow.io
One outlier sets the width of every equal-interval class, so 99% of counties land in the first. Counts from step 2.

3. Colours look random β€” the column is not numeric

The same file holds each county's 2010 census population. Read it the ordinary way and look at the type:

print(est["CENSUS2010POP"].dtype, est["POPESTIMATE2020"].dtype)
print(est.loc[pd.to_numeric(est["CENSUS2010POP"], errors="coerce").isna(),
              ["STNAME", "CTYNAME", "CENSUS2010POP"]].to_string(index=False))
str int64
STNAME                  CTYNAME CENSUS2010POP
Alaska      Chugach Census Area             X
Alaska Copper River Census Area             X

Two Alaska census areas that did not exist in 2010 carry an X, and two letters in 3,143 rows turn the whole column into text. A string column is treated as categorical β€” each distinct value gets its own arbitrary colour and any ordering is coincidence β€” and arithmetic on it fails outright:

try:
    est["POPESTIMATE2020"] - est["CENSUS2010POP"]
except TypeError as err:
    print(f"TypeError: {err}")


est["CENSUS2010POP"] = pd.to_numeric(est["CENSUS2010POP"], errors="coerce")
print(f"{est['CENSUS2010POP'].isna().sum()} values could not be parsed")

# now it can join the map
counties = counties.merge(est[["GEOID", "CENSUS2010POP"]].rename(columns={"CENSUS2010POP": "pop2010"}),
                          on="GEOID", how="left")
TypeError: operation 'rsub' not supported for dtype 'str' with dtype 'int64'
2 values could not be parsed

errors="coerce" turns unparseable values into NaN rather than raising, which is right here β€” but print the count, because a silent 40% failure produces a map of the 60% that parsed. This is ordinary attribute cleaning, and it belongs before the map, not in it.

4. High values look pale β€” the colormap is reversed

Matplotlib's _r suffix reverses a colormap, and some ramps are counter-intuitive without it:

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
for ax, cmap in zip(axes, ["YlOrRd", "YlOrRd_r"]):
    counties.plot(column="death_rate", scheme="quantiles", k=5, cmap=cmap, ax=ax, legend=True)
    ax.set_title(cmap); ax.set_axis_off()

Beyond the suffix, check the ramp's direction matches the message. YlOrRd, Blues, viridis all run light β†’ dark as values increase, which reads as "more". RdYlGn runs red β†’ green, which carries a value judgement: green as "good" is right for income and wrong for pollution. RdYlGn_r flips it.

5. Holes in the map β€” missing values are not being drawn

missing = counties["deaths"].isna()
print(f"{missing.sum()} of {len(counties):,} features have no value:")
print(counties.loc[missing, "STATE_NAME"].value_counts().to_string())
91 of 3,234 features have no value:
STATE_NAME
Puerto Rico                                     78
American Samoa                                   5
Commonwealth of the Northern Mariana Islands     4
United States Virgin Islands                     3
Guam                                             1

The boundary file covers every county-equivalent the Census Bureau maps, including Puerto Rico's municipios and the island areas; the county estimates file covers the fifty states and DC. By default, features with a null value are simply not drawn, so Puerto Rico vanishes from the map β€” a hole a reader interprets as zero, or as an error in the boundaries.

counties.plot(column="death_rate", scheme="quantiles", k=5, cmap="YlOrRd", legend=True,
              missing_kwds={"color": "#e2e8f0", "edgecolor": "#cbd5e1",
                            "hatch": "///", "label": "no data"})

The hatch matters as much as the grey: a flat grey can be read as another class, while a hatch reads as "not applicable" at a glance. And "label" puts it in the legend, which is the part that makes it unambiguous. See how to handle missing and null values.

Also check for implicit missing values β€” a 0 or -999 standing in for "unknown" is worse than a null, because it plots as a real class. A placeholder shows up as one value far more common than any other:

print(counties["deaths"].value_counts().head(3).to_string())
deaths
110.0    13
23.0     13
73.0     12

6. A diverging map whose neutral colour is not at zero

A diverging ramp centred on the data midpoint compared with one centred on zero. matplotlib centres on the range, not on zeroNo explicit rangeβœ•neutral colour lands at +16.8%βœ•+10% growth is coloured as declineβœ•the two sides are not comparableβœ•nothing warns youvmin=-lim, vmax=+limβœ“neutral colour sits at zeroβœ“equal colour, equal magnitudeβœ“growth and decline read alikeβœ“clip at p98 for the tailTwoSlopeNorm keeps zero neutral when the tails are truly uneven spatialworkflow.io
Without a symmetric range, the neutral colour lands wherever the data's midpoint happens to be.
counties["change"] = (counties["pop2020"] - counties["pop2010"]) / counties["pop2010"] * 100
c = counties["change"]
print(f"population change 2010–2020: {c.min():+.1f}% … {c.max():+.1f}%; "
      f"{(c < 0).sum():,} counties shrank, {(c > 0).sum():,} grew")
print(f"midpoint of the range: {(c.min() + c.max()) / 2:+.1f}%")
population change 2010–2020: -33.3% … +139.7%; 1,695 counties shrank, 1,445 grew
midpoint of the range: +53.2%

With a RdBu ramp and no explicit range, matplotlib puts the neutral colour at the midpoint of the data range β€” well above zero here, because the fastest-growing counties grew by far more than the fastest-shrinking ones shrank. So a county that grew modestly is coloured as if it declined. Force symmetry:

import numpy as np

lim = np.abs(c).quantile(0.98)                 # ignore the extreme tail
counties.plot(column="change", cmap="RdBu", vmin=-lim, vmax=lim, legend=True,
              legend_kwds={"label": "Population change 2010–2020 (%)"})
print(f"symmetric range Β±{lim:.1f}%")
symmetric range Β±26.8%

Or use TwoSlopeNorm when the tails are genuinely asymmetric and clipping would lose information:

from matplotlib.colors import TwoSlopeNorm
norm = TwoSlopeNorm(vmin=c.min(), vcenter=0, vmax=c.max())
counties.plot(column="change", cmap="RdBu", norm=norm, legend=True)

TwoSlopeNorm keeps zero at the neutral colour while using the full ramp on both sides β€” at the cost that a given colour distance means different amounts above and below zero. Say which you used.

7. The pattern follows population β€” you are mapping a count

If the map's high values are exactly the cities, check what is being coloured:

print(f"deaths vs population:     r = {counties['deaths'].corr(counties['pop2020']):.2f}")
print(f"death rate vs population: r = {counties['death_rate'].corr(counties['pop2020']):.2f}")
deaths vs population:     r = 0.98
death rate vs population: r = -0.21

A correlation near 1 between the mapped column and population means the map is a population map: a choropleth colours areas, and a raw count is largest where the denominator is largest. The death rate β€” which the estimates file publishes, per 1,000 residents β€” tells a different story, and its weak negative correlation with population is the first hint of it: small rural counties are older.

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

This is the single highest-value correction in the whole article, and it is not a plotting fix at all.

US county estimates, seven ways

The data used on this page

The US Census Bureau's county population estimates, vintage 2020 (co-est2020-alldata.csv), joined to its 2020 cartographic boundary file of counties (cb_2020_us_county_5m) on the five-digit FIPS code. The estimates give the 2010 census count, the July 2020 estimate, deaths in the year to July 2020 and the death rate per 1,000. The CENSUS2010POP column is left as the file delivers it; step 3 converts it and joins it to the map.

import geopandas as gpd
import pandas as pd

POPEST = "https://www2.census.gov/programs-surveys/popest/datasets/2010-2020/counties/totals/co-est2020-alldata.csv"
est = pd.read_csv(POPEST, encoding="latin-1", dtype={"STATE": str, "COUNTY": str})
est = est[est["SUMLEV"] == 50].copy()                        # counties, not state totals
est["GEOID"] = est["STATE"] + est["COUNTY"]

counties = gpd.read_file("https://www2.census.gov/geo/tiger/GENZ2020/shp/cb_2020_us_county_5m.zip",
                         columns=["GEOID", "NAMELSAD", "STATE_NAME"])
counties = counties.merge(
    est[["GEOID", "POPESTIMATE2020", "DEATHS2020", "RDEATH2020"]].rename(columns={
        "POPESTIMATE2020": "pop2020", "DEATHS2020": "deaths", "RDEATH2020": "death_rate"}),
    on="GEOID", how="left").to_crs(5070)                         # Albers equal-area
print(f"{len(est):,} counties in the estimates, {len(counties):,} polygons in the boundary file")
3,143 counties in the estimates, 3,234 polygons in the boundary file

Example 1: a diagnostic that names the cause

def diagnose_choropleth(gdf, column, *, scheme=None, k=5, cmap="YlOrRd",
                        denominator_candidates=("population", "pop2020", "pop")):
    problems = []
    v = gdf[column]

    if not pd.api.types.is_numeric_dtype(v):
        sample = v.dropna().astype(str).head(3).tolist()
        return [f"'{column}' is {v.dtype}, not numeric β€” e.g. {sample}. "
                f"Use pd.to_numeric(..., errors='coerce') first."]

    n_missing = v.isna().sum()
    if n_missing:
        problems.append(f"{n_missing:,} of {len(gdf):,} features have no value β€” "
                        f"pass missing_kwds or they will not be drawn")

    vv = v.dropna()
    if vv.empty:
        return [f"'{column}' is entirely null"]

    for sentinel in (0, -1, -999, -9999):
        n = (vv == sentinel).sum()
        if n and n / len(vv) > 0.02:
            problems.append(f"{n:,} values equal {sentinel} β€” a placeholder for missing?")

    if vv.max() > 0 and vv.quantile(0.95) > 0:
        ratio = vv.max() / vv.quantile(0.95)
        if ratio > 3:
            problems.append(f"max is {ratio:.1f}x the 95th percentile β€” one outlier will "
                            f"flatten a continuous ramp; classify or clip")

    if scheme is None:
        problems.append("no scheme passed β€” the continuous default is the least readable "
                        "option; try scheme='quantiles'")
    else:
        counts = [int(n) for n in mc.classify(vv, scheme, k=k).counts]
        share = 100 * max(counts) / sum(counts)
        if share > 60:
            problems.append(f"scheme '{scheme}' puts {share:.0f}% of features in one "
                            f"class (counts {counts}) β€” try quantiles")

    if (vv < 0).any() and (vv > 0).any():
        if cmap.removesuffix("_r") not in {"RdBu", "PiYG", "BrBG", "coolwarm", "RdYlBu", "Spectral"}:
            problems.append(f"values cross zero but '{cmap}' is a sequential ramp β€” "
                            f"use a diverging one with symmetric vmin/vmax")

    for cand in denominator_candidates:
        if cand in gdf.columns and cand != column and pd.api.types.is_numeric_dtype(gdf[cand]):
            r = vv.corr(gdf.loc[vv.index, cand])
            if r > 0.8:
                problems.append(f"'{column}' correlates {r:.2f} with '{cand}' β€” this may "
                                f"be a count; map {column}/{cand} instead")
            break

    return problems or ["no obvious problem β€” check the colormap direction by eye"]


for column, scheme, cmap in [("deaths", "equal_interval", "YlOrRd"), ("change", None, "YlOrRd"),
                             ("death_rate", "quantiles", "YlOrRd")]:
    print(f"{column}, scheme={scheme}, cmap={cmap}:")
    for p in diagnose_choropleth(counties, column, scheme=scheme, cmap=cmap):
        print(f"  β†’ {p}")
deaths, scheme=equal_interval, cmap=YlOrRd:
  β†’ 91 of 3,234 features have no value β€” pass missing_kwds or they will not be drawn
  β†’ max is 17.5x the 95th percentile β€” one outlier will flatten a continuous ramp; classify or clip
  β†’ scheme 'equal_interval' puts 99% of features in one class (counts [3123, 16, 2, 1, 1]) β€” try quantiles
  β†’ 'deaths' correlates 0.98 with 'pop2020' β€” this may be a count; map deaths/pop2020 instead
change, scheme=None, cmap=YlOrRd:
  β†’ 93 of 3,234 features have no value β€” pass missing_kwds or they will not be drawn
  β†’ max is 7.7x the 95th percentile β€” one outlier will flatten a continuous ramp; classify or clip
  β†’ no scheme passed β€” the continuous default is the least readable option; try scheme='quantiles'
  β†’ values cross zero but 'YlOrRd' is a sequential ramp β€” use a diverging one with symmetric vmin/vmax
death_rate, scheme=quantiles, cmap=YlOrRd:
  β†’ 91 of 3,234 features have no value β€” pass missing_kwds or they will not be drawn

Four findings for the raw death count, four for the population change, and for the death rate only the missing values β€” the one problem left once the column is a rate, classified, on a sequential ramp. The count-versus-population finding is the important one β€” the others are about how the map looks, and that one is about what it says. The thresholds are heuristics chosen to be loud rather than precise; a false positive costs a moment's thought, and a false negative ships. (str.rstrip("_r") would be the tempting way to drop the reversal suffix; it strips characters, so it turns "Spectral" into "Spect". removesuffix is the method that means what it says.)

What the choropleth diagnostic reported for three US county columns. diagnose_choropleth on three county columnscolumnfindingswhat they weredeaths, equal_interval4nulls, outlier, one class, a countchange, no scheme4nulls, outlier, no scheme, crosseszerodeath_rate, quantiles1nulls onlyA rate, classified, on the right ramp leaves only the holes to label spatialworkflow.io
The output of example 1, summarised. Every finding on the first two rows was a real property of the county data, not a planted one.

Example 2: the same data, seven ways

Seeing the failure modes together makes them recognisable afterwards:

import matplotlib.pyplot as plt


def failure_gallery(gdf, column, denominator=None):
    v = gdf[column]
    lo, hi = v.quantile([0.02, 0.98])
    panels = [
        ("default β€” continuous", dict(column=column, cmap="YlOrRd")),
        ("equal_interval", dict(column=column, cmap="YlOrRd", scheme="equal_interval", k=5)),
        ("quantiles", dict(column=column, cmap="YlOrRd", scheme="quantiles", k=5)),
        ("fisher_jenks", dict(column=column, cmap="YlOrRd", scheme="fisher_jenks", k=5)),
        ("clipped 2–98%", dict(column=column, cmap="YlOrRd", vmin=lo, vmax=hi)),
        ("reversed ramp", dict(column=column, cmap="YlOrRd_r", scheme="quantiles", k=5)),
        ("nulls shown", dict(column=column, cmap="YlOrRd", scheme="quantiles", k=5,
                             missing_kwds={"color": "#e2e8f0", "hatch": "///"})),
    ]
    if denominator:
        gdf = gdf.assign(_rate=gdf[column] / gdf[denominator].replace(0, np.nan))
        panels.append(("as a rate", dict(column="_rate", cmap="YlOrRd", scheme="quantiles", k=5)))

    fig, axes = plt.subplots(2, 4, figsize=(17, 9))
    for ax, (title, kwargs) in zip(axes.ravel(), panels):
        gdf.plot(ax=ax, edgecolor="white", linewidth=0.05, **kwargs)
        ax.set_title(title, fontsize=9)
        ax.set_axis_off()
    for ax in axes.ravel()[len(panels):]:
        ax.set_axis_off()
    plt.tight_layout()
    return fig


conus = counties[~counties["STATE_NAME"].isin(["Alaska", "Hawaii", "Puerto Rico", "Guam", "American Samoa",
                                               "Commonwealth of the Northern Mariana Islands",
                                               "United States Virgin Islands"])]
failure_gallery(conus, "deaths", denominator="pop2020").savefig("gallery.png", dpi=70)

Compare the first panel with the third and the last. The first is nearly uniform; the third is legible but shows where people are; the last shows where deaths were high relative to population, which is what a reader assumes the map meant all along.

Example 3: a plotting function that refuses the common mistakes

DIVERGING = {"RdBu", "PiYG", "BrBG", "coolwarm", "RdYlBu", "Spectral", "PuOr"}


def safe_choropleth(gdf, column, *, denominator=None, per=1, scheme="quantiles",
                    k=5, cmap=None, crs=None, title=None, clip=None, ax=None):
    gdf = gdf.copy()

    if not pd.api.types.is_numeric_dtype(gdf[column]):
        raise TypeError(f"'{column}' is {gdf[column].dtype}; convert with pd.to_numeric first")

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

    if crs:
        gdf = gdf.to_crs(crs)
    elif gdf.crs is not None and gdf.crs.is_geographic:
        raise ValueError("geographic CRS β€” pass crs= (ideally equal-area) for a choropleth")

    v = gdf[col]
    crosses_zero = bool((v.dropna() < 0).any() and (v.dropna() > 0).any())
    cmap = cmap or ("RdBu" if crosses_zero else "YlOrRd")
    if crosses_zero and cmap.removesuffix("_r") not in DIVERGING:
        raise ValueError(f"values cross zero β€” '{cmap}' is sequential; use a diverging ramp")

    kwargs = dict(cmap=cmap, edgecolor="white", linewidth=0.05, legend=True,
                  missing_kwds={"color": "#e2e8f0", "edgecolor": "#cbd5e1", "hatch": "///",
                                "label": f"no data ({v.isna().sum():,})"})

    if crosses_zero:                      # symmetric range so neutral sits at zero
        lim = np.abs(v).quantile(clip or 0.98)
        kwargs.update(vmin=-lim, vmax=lim, legend_kwds={"label": col.replace("_", " ")})
    elif clip:
        lo, hi = v.quantile([clip, 1 - clip])
        kwargs.update(vmin=lo, vmax=hi, legend_kwds={"label": col.replace("_", " ")})
    else:
        kwargs.update(scheme=scheme, k=k,
                      legend_kwds={"loc": "lower right", "fmt": "{:,.1f}",
                                   "title": col.replace("_", " "), "fontsize": 9})

    if ax is None:
        _, ax = plt.subplots(figsize=(12, 8))
    gdf.plot(column=col, ax=ax, **kwargs)
    ax.set_title(title or col, fontsize=14, loc="left")
    ax.set_axis_off()
    return ax, cmap


for args in [dict(column="pop2010"), dict(column="change", cmap="YlOrRd")]:
    try:
        safe_choropleth(conus.assign(pop2010=conus["pop2010"].astype(str)) if args["column"] == "pop2010" else conus, **args)
    except (TypeError, ValueError) as err:
        print(f"{type(err).__name__}: {err}")

ax, chosen = safe_choropleth(conus, "deaths", denominator="pop2020", per=1_000,
                             title="Deaths per 1,000 residents, year to July 2020")
ax.figure.savefig("deaths_rate.png", dpi=70)
ax, chosen = safe_choropleth(conus, "change", title="Population change 2010–2020 (%)")
print(f"change map drawn with cmap={chosen!r}, centred on zero")
TypeError: 'pop2010' is str; convert with pd.to_numeric first
ValueError: values cross zero β€” 'YlOrRd' is sequential; use a diverging ramp
change map drawn with cmap='RdBu', centred on zero

Four refusals and two automatic choices. It raises on a non-numeric column and on a geographic CRS, both of which produce a plausible-looking wrong map. It picks a diverging ramp automatically when the values cross zero, and refuses a sequential one in that case. It always shows missing values, with a count in the legend. And when the range is symmetric it centres the neutral colour on zero rather than on the data's midpoint.

Raising rather than warning is deliberate. A warning in a notebook scrolls past; the map does not.

Three mappings, one of them failing

The three mappings a choropleth performs and the failure at each level. three mappings, three kinds of wrongvalue to numberthe datastrings, sentinels, countsnumber to classthe schemeskew versus a linear rampclass to colourthe rampdirection and midpointOnly the bottom layer is about colour spatialworkflow.io
Three mappings, each with its own failure. Only the last one is about colour.

A choropleth is a composition of three independent mappings, and "the colours look wrong" is always a failure in one of them.

Value β†’ number. The column has to be numeric and has to mean what you think. A string column silently becomes categorical; a -999 placeholder becomes a legitimate extreme; a count becomes a proxy for population. None of these are plotting problems, and none of them can be fixed by changing the colormap.

Number β†’ class. This is where most of the visible failures live, and the reason is the shape of real spatial data. Income, density, prices, case counts and property values are right-skewed almost without exception, and the default continuous ramp is a linear map from value to colour. Linear maps and skewed data do not mix: the top 1% of values claim most of the colour range and the other 99% share what is left. Classification is the fix because it replaces a linear map with one derived from the distribution β€” and different schemes derive it differently, which is why equal_interval and quantiles produce such different maps from identical numbers.

Class β†’ colour. The smallest of the three and still worth care. Sequential ramps encode magnitude and must run in the direction the reader expects. Diverging ramps encode deviation from a meaningful centre, and that centre must actually be at the centre β€” matplotlib puts the neutral colour at the midpoint of the range, which is only zero by coincidence. And a categorical palette on ordered data destroys the ordering that the data has.

Missing values are a fourth thing that behaves like a colour bug. GeoPandas does not draw features whose value is null, so they vanish. A reader seeing a hole in a choropleth reads it as zero, or as a gap in the boundaries β€” not as "unknown". missing_kwds with a hatch and a labelled legend entry is the difference between an honest absence and an accidental claim. The same reasoning applies to placeholder values masquerading as data, which is why they are worth checking for explicitly.

The deepest of the failures is the count-versus-rate one, because it is not about appearance at all. A choropleth encodes value as colour across an area, and human perception integrates colour over area whether or not you intend it. Mapping a count therefore produces a picture of where the denominator is large. The map is not badly rendered; it is showing a different variable from the one in its legend. Dividing by an appropriate denominator fixes more misleading choropleths than every other item in this article combined β€” the argument is developed in choropleth classification explained.

Distinct colours drawn for ten values in five quantile classes, depending on vmin and vmax. scheme="quantiles", k=5, ten values from 0 to 9,000distinct colourswhat happenedno vmin/vmaxβœ“5one per classvmin=-100, vmax=100!3classes 0–4 squeezed mid-rampvmin=0, vmax=9000βœ•1classes 0–4 all at the bottomWith a scheme, vmin and vmax apply to class numbers, not to values spatialworkflow.io
Measured on a ten-polygon test layer. Passing the data range to a classified plot β€” a natural thing to try β€” draws every class the same colour.

vmin with schemes, nullable integers and TwoSlopeNorm

  • scheme= needs mapclassify installed, or GeoPandas raises.

  • vmin/vmax do not mix with scheme. They are not ignored: with a scheme set they are applied to the class numbers (0 to kβˆ’1), not to your values. Passing the data's own range β€” vmin=0, vmax=9000 on values up to 9,000 β€” squeezes all five classes into the first sliver of the ramp and draws every polygon the same colour. Classification and continuous normalisation are alternatives, not layers.

  • _r reverses any colormap. YlOrRd_r makes high values pale.

  • missing_kwds also accepts hatch and label, and both help β€” grey alone can read as another class.

  • Integer columns with nulls become floats. Pandas' nullable Int64 keeps them integers if that matters for the legend.

  • TwoSlopeNorm keeps zero neutral with asymmetric tails, but colour distance then means different amounts on each side.

  • Categorical data needs categorical=True and a qualitative palette such as tab10; a sequential ramp implies an order that is not there.

  • Colour-blind readers cannot separate red from green. viridis, cividis and ColorBrewer's safe ramps can.

  • Very small polygons are invisible whatever the colour. Consider a cartogram or proportional symbols.

  • legend_kwds={"fmt": "{:,.0f}"} removes 50000.00000 from a classified legend.

FAQ

Why is my whole map one colour?

Either no classification scheme was passed and an outlier is flattening the continuous ramp, or equal_interval on skewed data has put nearly everything in the first class. Check the class counts.

Why are my colours in a random order?

The column is not numeric, so GeoPandas is treating it as categorical. Convert with pd.to_numeric(..., errors='coerce') and check how many values failed to parse.

Why are there holes in my map?

Features with null values are not drawn. Pass missing_kwds={"color": "#e2e8f0", "hatch": "///", "label": "no data"} so absence is visible and labelled.

Why is the neutral colour of my diverging map not at zero?

Matplotlib centres the ramp on the midpoint of the data range, not on zero. Pass symmetric vmin=-lim, vmax=lim, or use TwoSlopeNorm(vcenter=0).

My map just shows where the cities are.

You are mapping a count. A choropleth colours areas, so a count shows the denominator β€” usually population. Map a rate instead.

vmin and vmax are doing something strange.

With a scheme set they are applied to the class numbers, not to your values, so the data's own range squeezes every class into one colour. Classification replaces continuous normalisation β€” use one or the other.

Which colormap should I use?

Sequential (YlOrRd, Blues, viridis) for magnitude, diverging (RdBu, PiYG) for deviation from a centre, qualitative (tab10) for categories. Prefer colour-blind-safe ramps for anything published.