How to Make a Map Series with Consistent Symbology

Problem statement

A report needs one map per region, or one per year, or one per scenario. The obvious loop produces maps that cannot be compared:

import mapclassify

for year, subset in series.groupby("year"):
    subset.plot(column="rate", cmap="YlGnBu", scheme="Quantiles", k=5, legend=True)   # wrong
    plt.savefig(f"change-{year}.png")
    plt.close()
    darkest = mapclassify.Quantiles(subset["rate"], k=5).bins[-2]
    print(f"{year}: darkest class starts at {darkest:+.2f}%")

2021: darkest class starts at +1.25%
2022: darkest class starts at +1.53%
2023: darkest class starts at +1.84%
2024: darkest class starts at +1.49%
2025: darkest class starts at +1.02%

This is the annual population change of Great Britain's 350 local authority districts, one map per year from 2021 to 2025 (the data is set up under Code examples). The darkest blue means "above +1.25%" in 2021, "above +1.84%" in 2023 and "above +1.02%" in 2025.

Every panel gets its own colour scale, computed from its own data, so the darkest colour means something different on every map. Readers compare the colours anyway, because that is what colours are for.

The same happens with class breaks, legend order, extent, aspect ratio, label density and the position of the apparatus. A map series is defined by everything that stays the same.

Quick answer

Compute the shared decisions once, before the loop, and pass them in:

import mapclassify


def series_style(gdf, column, k=5, scheme="Quantiles"):
    """Everything that must be identical across the panels."""
    classifier = getattr(mapclassify, scheme)(gdf[column].dropna(), k=k)
    return {
        "bins": [round(float(b), 2) for b in classifier.bins],
        "vmin": float(gdf[column].min()),
        "vmax": float(gdf[column].max()),
        "cmap": "YlGnBu",
        "figsize": (170 / 25.4, 120 / 25.4),
    }


style = series_style(series, "rate")             # once, on the whole dataset
print(style["bins"])

[0.15, 0.57, 0.99, 1.46, 38.34]

for year, subset in series.groupby("year"):
    plot_panel(subset, style, title=year)        # many times, identically

The classifier is fitted to all the data, not to each panel. That single change is what makes the panels comparable.

Share the breaks, the extent and the panel function

1. List what must be shared

For a comparable series, all of these are computed once:

Extent is the interesting exception: a series of regions usually needs each panel zoomed to its own region, while a series of years usually needs one fixed extent.

2. Fit the classifier on the whole dataset

mapclassify fits breaks to whatever you give it. Give it everything:

import numpy as np

classifier = mapclassify.Quantiles(series["rate"].dropna(), k=5)
bins = classifier.bins
for year, subset in series.groupby("year"):
    shared = np.digitize(subset["rate"], bins, right=True)
    own = np.digitize(subset["rate"], mapclassify.Quantiles(subset["rate"], k=5).bins, right=True)
    print(f"{year}: {(shared != own).mean():5.1%} of districts change class; "
          f"shared-bin counts {np.bincount(shared, minlength=5).tolist()}")

2021: 52.1% of districts change class; shared-bin counts [134, 79, 44, 51, 53]
2022: 18.6% of districts change class; shared-bin counts [51, 73, 77, 80, 80]
2023: 70.9% of districts change class; shared-bin counts [13, 53, 78, 93, 124]
2024: 23.5% of districts change class; shared-bin counts [46, 63, 91, 84, 77]
2025: 62.9% of districts change class; shared-bin counts [115, 91, 68, 51, 25]

Then apply the same bins to each panel with scheme="UserDefined". Quantiles fitted per panel are the most common way a series becomes misleading, precisely because quantiles are the most common default: here they put between a fifth and seven-tenths of the districts in a different class from the shared scale, and they hide what the shared counts show — 2023 had 124 districts in the fastest-growing class and 2025 had 25.

3. Decide the extent policy deliberately

Three defensible policies, and they answer different questions:

• Shared extent — every panel shows the same window. Right for a time series; makes spatial change obvious.
• Per-panel extent, shared scale — each panel is centred on its own region but at the same map scale, so a bigger region simply extends beyond the frame or leaves more white space. Right for comparing regions of similar size.
• Per-panel extent, per-panel scale — each region fills its frame. Right when the shapes matter more than the sizes, and it must be said in the caption because readers assume a shared scale.

4. Build one panel function and call it in a loop

The loop should contain no styling decisions. Everything is in the style dictionary or in the function.

This is also what makes the series maintainable: a change to the legend is one edit, not one edit per panel.

5. Handle the empty and the extreme panels

Two cases that break a naive loop:

• A panel with no data in some classes. The legend must still show every class, or the reader will think the class does not exist.
• A panel whose values are all in one class. It will be a single flat colour, which is correct and looks broken. A note in the panel — "all districts in class 3" — prevents the misreading.

6. Render the legend once if the panels share it

For a grid of small multiples in one figure, one legend serves all panels and saves the space that repeating it would consume. For separate files, each needs its own, because they will be seen apart.

Six years of UK population on one scale

The data used on this page

The UK's local authority districts (ONS, December 2024, ultra-generalised) and the ONS mid-year population estimates for 2020 to 2025 from NOMIS. Each district gets one row per year from 2021 to 2025 with rate, the percentage change on the year before. Northern Ireland has no mid-2025 value in the table, so 2025 has 350 districts and the other years 361.

import os

import geopandas as gpd
import matplotlib
import pandas as pd
import requests

matplotlib.use("Agg")
import matplotlib.pyplot as plt

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Local_Authority_Districts_December_2024_Boundaries_UK_BUC/FeatureServer/0/query")
NOMIS = ("https://www.nomisweb.co.uk/api/v01/dataset/NM_2002_1.data.csv?geography=TYPE424&gender=0"
         "&c_age=200&measures=20100&select=date_name,geography_code,obs_value&date=")
if not os.path.exists("districts.gpkg"):
    reply = requests.get(ONS, timeout=120, params={
        "where": "1=1", "outFields": "LAD24CD,LAD24NM", "outSR": 27700, "f": "geojson"})
    reply.raise_for_status()
    gpd.read_file(reply.text).to_file("districts.gpkg")
    pd.concat([pd.read_csv(NOMIS + str(year)) for year in range(2020, 2026)]).to_csv("population.csv", index=False)

districts = gpd.read_file("districts.gpkg")
population = pd.read_csv("population.csv").pivot(index="GEOGRAPHY_CODE", columns="DATE_NAME", values="OBS_VALUE")
series = gpd.GeoDataFrame(pd.concat([
    districts.merge(((population[y] / population[y - 1] - 1) * 100).rename("rate"),
                    left_on="LAD24CD", right_index=True).assign(year=y)
    for y in range(2021, 2026)]), crs=districts.crs).dropna(subset=["rate"])
print(series.groupby("year")["rate"].agg(["count", "min", "median", "max"]).round(2))

      count   min  median    max
year                            
2021    361 -2.41    0.40  11.87
2022    361 -0.50    0.86  38.34
2023    361 -0.90    1.19  20.75
2024    361 -3.33    0.90   6.77
2025    350 -3.52    0.44   2.39

Example 1 — the shared style object

from dataclasses import dataclass, field
import mapclassify
import matplotlib.pyplot as plt

MM = 1 / 25.4


@dataclass
class SeriesStyle:
    column: str
    bins: list
    cmap: str = "YlGnBu"
    vmin: float | None = None
    vmax: float | None = None
    figsize: tuple = (90 * MM, 90 * MM * 0.8)
    fontsize: float = 7
    boundary_width: float = 0.35
    missing_colour: str = "#e2e8f0"
    extent: tuple | None = None            # None = per-panel extent
    labels: list = field(default_factory=list)

    @classmethod
    def fit(cls, gdf, column, k=5, scheme="Quantiles", **kwargs):
        values = gdf[column].dropna()
        classifier = getattr(mapclassify, scheme)(values, k=k)
        bins = list(classifier.bins)
        labels = []
        lower = values.min()
        for upper in bins:
            labels.append(f"{lower:,.1f} – {upper:,.1f}")
            lower = upper
        return cls(column=column, bins=bins, vmin=float(values.min()),
                   vmax=float(values.max()), labels=labels, **kwargs)


def plot_panel(gdf, style: SeriesStyle, title=None, ax=None, show_legend=True):
    fig = None
    if ax is None:
        fig, ax = plt.subplots(figsize=style.figsize)

    gdf.plot(column=style.column, cmap=style.cmap, ax=ax,
             scheme="UserDefined", classification_kwds={"bins": style.bins},
             edgecolor="white", linewidth=style.boundary_width,
             legend=show_legend,
             legend_kwds={"fontsize": style.fontsize - 0.5, "frameon": False,
                          "loc": "lower left", "labels": style.labels},
             missing_kwds={"color": style.missing_colour, "label": "no data"})

    if style.extent:
        ax.set_xlim(style.extent[0], style.extent[2])
        ax.set_ylim(style.extent[1], style.extent[3])
    ax.set_aspect("equal")
    ax.set_axis_off()
    if title:
        ax.set_title(title, fontsize=style.fontsize + 1, loc="left")
    return fig, ax

style = SeriesStyle.fit(series, "rate", k=5)
print(style.labels)
fig, ax = plot_panel(series[series["year"] == 2023], style, title="2023")
fig.savefig("change-2023.png", dpi=150)

['-3.5 – 0.1', '0.1 – 0.6', '0.6 – 1.0', '1.0 – 1.5', '1.5 – 38.3']

Example 2 — small multiples in one figure with one shared legend

import math
import matplotlib.pyplot as plt


def small_multiples(gdf, style, facet_col, ncols=4, share_extent=True):
    facets = sorted(gdf[facet_col].dropna().unique())
    nrows = math.ceil(len(facets) / ncols)

    fig, axes = plt.subplots(nrows, ncols,
                             figsize=(style.figsize[0] * ncols,
                                      style.figsize[1] * nrows))
    axes = axes.ravel() if hasattr(axes, "ravel") else [axes]

    extent = tuple(gdf.total_bounds) if share_extent else None
    for ax, facet in zip(axes, facets):
        subset = gdf[gdf[facet_col] == facet]
        panel_style = style.__class__(**{**style.__dict__, "extent": extent})
        plot_panel(subset, panel_style, title=str(facet), ax=ax, show_legend=False)
        classes = np.digitize(subset[style.column].dropna(), style.bins, right=True)
        if len(subset) and np.unique(classes).size == 1:
            ax.text(0.5, 0.02, "all areas in one class", transform=ax.transAxes,
                    ha="center", fontsize=style.fontsize - 1.5, color="#64748b")

    for ax in axes[len(facets):]:
        ax.set_visible(False)

    handles, labels = build_class_handles(style)
    fig.legend(handles, labels, loc="lower center", ncols=len(labels),
               frameon=False, fontsize=style.fontsize - 0.5,
               bbox_to_anchor=(0.5, -0.01))
    return fig, axes


def build_class_handles(style):
    import matplotlib.patches as mpatches
    cmap = plt.get_cmap(style.cmap, len(style.bins))
    handles = [mpatches.Patch(facecolor=cmap(i), edgecolor="white", linewidth=0.5)
               for i in range(len(style.bins))]
    return handles, style.labels

fig, axes = small_multiples(series, style, "year", ncols=5)
fig.savefig("change-2021-2025.png", dpi=120, bbox_inches="tight")
print(len(fig.axes), "axes;", sum(ax.get_visible() for ax in axes), "panels drawn with one legend")

5 axes; 5 panels drawn with one legend

Example 3 — verifying the series is actually consistent

def check_series(paths_or_figs, style):
    """Catch the panels that drifted."""
    import matplotlib.pyplot as plt

    problems = []
    sizes, extents = set(), set()

    for i, fig in enumerate(paths_or_figs):
        sizes.add(tuple(round(float(v), 3) for v in fig.get_size_inches()))
        for ax in fig.axes:
            if ax.get_images() or ax.collections:
                extents.add((round(ax.get_xlim()[1] - ax.get_xlim()[0], 1),
                             round(ax.get_ylim()[1] - ax.get_ylim()[0], 1)))
        for ax in fig.axes:                                  # only text a reader sees:
            legend = ax.get_legend()                         # titles, annotations, legend
            for text in [ax.title, *ax.texts, *(legend.get_texts() if legend else [])]:
                if text.get_visible() and text.get_text().strip():
                    if abs(text.get_fontsize() - style.fontsize) > 1.5:
                        problems.append(f"panel {i}: '{text.get_text()}' at {text.get_fontsize()} pt")

    if len(sizes) > 1:
        problems.append(f"{len(sizes)} different figure sizes: {sizes}")
    if len(extents) > 1 and style.extent is not None:
        problems.append(f"{len(extents)} different extents despite a shared extent")

    print(f"{len(problems)} consistency problem(s)")
    for p in problems[:10]:
        print("  !", p)
    return problems

shared = SeriesStyle.fit(series, "rate", k=5, extent=tuple(districts.total_bounds))
figures = [plot_panel(subset, shared, title=str(year))[0] for year, subset in series.groupby("year")]
check_series(figures, shared)

figures[2].set_size_inches(5, 4)                   # one panel exported at another size
figures[4].axes[0].set_title("2025", fontsize=12)  # and one with a bigger title
check_series(figures, shared)

0 consistency problem(s)
2 consistency problem(s)
  ! panel 4: '2025' at 12.0 pt
  ! 2 different figure sizes: {(5.0, 4.0), (3.543, 2.835)}

Run it on the figures just before they are saved. It catches the drift that creeps in when one panel is edited by hand — a resized figure, a restyled title — which is how a series stops being one. Check only text the reader sees: an earlier version of this function walked every Text object in the figure and reported 125 problems on five identical panels, because the axes keep their tick-label text objects, invisible and at the default 10 pt, even when switched off.

Why per-panel classification misleads

Why per-panel classification is the central error

Quantiles, natural breaks and equal intervals are all fitted to the data they are given. Fitted per panel, they guarantee that each map uses the full ramp — which is exactly what makes the panels incomparable, because "dark" now means the top fifth of that panel.

Readers do not know this and would not expect it. A series of choropleths is read as a single map cut into pieces, and the colour is assumed to be the shared language.

The fix is one line — fit on the whole dataset, apply UserDefined breaks per panel — and it is the difference between a series and a collection.

Why the extent policy has to be explicit

A series of years with a shared extent shows change; a series of regions with per-panel extents shows each region. Both are correct, and mixing them silently is not.

The case that misleads is per-panel extent with per-panel scale, where a small region is enlarged to fill the frame. It looks like the other panels, and readers assume the scale is shared. That series needs the scale stated on each panel or in the caption.

Why the legend must show empty classes

A panel where no area falls in the top class will, by default, produce a legend without that class. The reader then compares a four-class panel with a five-class one and concludes the classification changed.

Building the legend from the shared style rather than from the panel's data keeps every class visible, with the empty ones simply unused.

Why a style object beats a function with many arguments

A series accumulates shared decisions: breaks, ramp, sizes, widths, the missing-data colour, the legend labels. Passed as keyword arguments, they drift — one call gets a different linewidth and nobody notices until the panels are side by side.

A single object, constructed once and passed down, makes drift impossible for anything it covers, and makes the shared decisions reviewable in one place.

Shared breaks in GeoPandas, diverging and categorical series

• scheme="UserDefined" with classification_kwds={"bins": bins} is how you apply shared breaks in GeoPandas.
• missing_kwds ensures no-data areas are drawn and labelled rather than left blank.
• Diverging series need a shared, symmetric norm anchored at the midpoint, not per-panel limits.
• Categorical series need a shared category order, or the legend reorders between panels.
• Outliers in one panel widen the shared scale for everybody. Consider a top class of "and over".
• Label placement cannot be shared — it is computed in pixels per panel.
• Render one panel first and check it before looping over sixty.
• File names should sort in the order the panels should be read.

Internal links

• Choropleth classification explained — the class breaks being shared
• How to build a reusable map style module — where the style object lives
• How to generate map images in batch — running the loop unattended
• Colour on maps explained: sequential, diverging and qualitative — choosing the shared ramp
• How to build a print-ready map layout in Matplotlib — the panel layout
• How to make a choropleth map in GeoPandas — the single-panel case
• Fixing choropleth colours that look wrong — when the classification is the problem
• How to label map features without overlaps in Matplotlib — per-panel labelling

FAQ

Why do my map panels use different colour scales?

Because the classifier was fitted per panel. Fit it once on the whole dataset and apply the same breaks to every panel with scheme="UserDefined".

Should every panel show the same extent?

For a time series, yes. For a series of regions, per-panel extents are usual — but keep the scale shared, or state in the caption that it is not.

How do I keep the legend the same across panels?

Build it from the shared style object rather than from each panel's data, so classes with no members in a panel still appear.

What if one panel's values are all in one class?

It will be a flat colour, which is correct and looks like an error. Add a short note in the panel saying so.

How do I stop the panels drifting?

Put every shared decision in one style object constructed before the loop, and check the rendered panels for differing figure sizes, extents and type sizes.

Can the panels share one legend?

In a small-multiples figure, yes — one legend for the whole figure saves the space that repetition would consume. Separate files each need their own.