Diverging Colormap Selection for Bivariate Thematic Maps

For a single variable that diverges from a meaningful reference value, pick a luminance-balanced diverging map such as colorcet CET-D and pin its neutral midpoint to that reference with matplotlib.colors.TwoSlopeNorm; for two variables at once, bin each into three classes and blend two hue ramps into a 3x3 bivariate color matrix, then render the grid itself as the legend. The two failure modes to avoid are an unbalanced ramp — where one arm is darker and falsely reads as more extreme — and a midpoint that floats to the data mean instead of the value that actually matters.

Core Algorithm and Workflow

A diverging colormap encodes deviation in two directions from a centre. It only tells the truth when two conditions hold: the centre is anchored to a value the reader cares about (zero change, the national average, a break-even threshold), and the two arms are perceptually symmetric so equal deviations look equally strong. A bivariate map extends this to two variables by treating colour as a two-dimensional signal — one hue axis per variable — which is why a one-dimensional colour bar cannot describe it and a 3x3 grid must.

3x3 bivariate color grid with two variable axes A three by three matrix of nine cells. Variable A increases in fill opacity from left to right along the horizontal axis; Variable B increases in fill opacity from bottom to top along the vertical axis. The top-right cell, where both variables are high, is the most saturated. Arrows label both axes. Variable B: vacancy rate high mid low Variable A: median income low mid high 3x3 bivariate color matrix (cell index = 3*B + A)

The workflow splits by how many variables you are mapping:

  1. Single-variable divergence. Choose a perceptually uniform, luminance-symmetric diverging map — colorcet’s CET-D family is built for this — and anchor its centre with TwoSlopeNorm(vcenter=anchor). The anchor is the value that separates “above” from “below,” never the data mean unless the mean is itself meaningful.
  2. Two-variable maps. Classify each variable into three ordered bins, combine them into a single index 3 * bin_B + bin_A running 0…8, blend two hue ramps into a 3x3 color matrix, and colour each feature by its cell. The matching legend is the grid itself, drawn with both axes labelled.

This page sits under Color Theory for GIS; the sequential-scale companion, Generating WCAG-Compliant Sequential Colormaps with Colorcet, covers the single-hue case where the data does not diverge.

Production-Ready Python Implementation

The script below has two independent halves. The first renders a single diverging variable with an anchored midpoint. The second builds a 3x3 bivariate matrix by blending two hue ramps in linear RGB and assigns every feature a cell colour. It uses matplotlib 3.9, colorcet 3.1, numpy 2.0, and geopandas 1.0.

import colorcet as cc
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import LinearSegmentedColormap, TwoSlopeNorm, to_rgb


# ---------------------------------------------------------------------------
# 1. Single variable: luminance-balanced diverging map, anchored midpoint
# ---------------------------------------------------------------------------
def plot_diverging(gdf: gpd.GeoDataFrame, value_col: str, anchor: float = 0.0):
    """
    Render a diverging choropleth whose neutral centre is pinned to `anchor`.

    CET-D1 is perceptually uniform and luminance-symmetric: both arms reach
    equal lightness contrast against the neutral midpoint, so equal deviations
    read as equally strong. `TwoSlopeNorm` keeps the neutral colour locked to
    the anchor even when the two arms span different data ranges.
    """
    cmap = cc.cm["CET_D1"]  # blue<->neutral<->red, balanced luminance
    vmin, vmax = gdf[value_col].min(), gdf[value_col].max()

    # Guard: TwoSlopeNorm requires vmin < vcenter < vmax strictly.
    if not (vmin < anchor < vmax):
        raise ValueError(
            f"anchor={anchor} must lie strictly between "
            f"data bounds ({vmin}, {vmax})"
        )

    norm = TwoSlopeNorm(vmin=vmin, vcenter=anchor, vmax=vmax)

    fig, ax = plt.subplots(figsize=(9, 9))
    gdf.plot(column=value_col, cmap=cmap, norm=norm, ax=ax,
             edgecolor="0.85", linewidth=0.2, legend=True,
             legend_kwds={"label": f"{value_col} (anchored at {anchor})",
                          "shrink": 0.6})
    ax.set_axis_off()
    return fig


# ---------------------------------------------------------------------------
# 2. Two variables: 3x3 bivariate color matrix from two blended hue ramps
# ---------------------------------------------------------------------------
def build_bivariate_matrix(hue_a: str = "#2166ac",
                           hue_b: str = "#b2182b",
                           n: int = 3) -> np.ndarray:
    """
    Construct an (n*n, 3) RGB matrix by blending two hue ramps.

    Variable A ramps from white toward `hue_a` across columns; variable B ramps
    from white toward `hue_b` across rows. The two contributions are combined
    in linear-light RGB (gamma-decoded) so the mid-cells do not turn muddy, the
    classic failure of naive sRGB averaging.

    Returns rows ordered by index = n * bin_b + bin_a, so cell (a, b) is
    matrix[n * b + a].
    """
    def _srgb_to_linear(c: np.ndarray) -> np.ndarray:
        return np.where(c <= 0.04045, c / 12.92, ((c + 0.055) / 1.055) ** 2.4)

    def _linear_to_srgb(c: np.ndarray) -> np.ndarray:
        c = np.clip(c, 0.0, 1.0)
        return np.where(c <= 0.0031308, c * 12.92,
                        1.055 * c ** (1 / 2.4) - 0.055)

    rgb_a = _srgb_to_linear(np.array(to_rgb(hue_a)))
    rgb_b = _srgb_to_linear(np.array(to_rgb(hue_b)))
    white = np.ones(3)

    matrix = np.zeros((n * n, 3))
    for b in range(n):        # variable B -> rows
        for a in range(n):    # variable A -> columns
            fa = a / (n - 1)   # 0..1 strength of A
            fb = b / (n - 1)   # 0..1 strength of B
            # Interpolate each hue from white toward its endpoint, then average
            lin_a = white + fa * (rgb_a - white)
            lin_b = white + fb * (rgb_b - white)
            lin_mix = (lin_a + lin_b) / 2.0
            matrix[n * b + a] = _linear_to_srgb(lin_mix)
    return matrix


def classify_bivariate(gdf: gpd.GeoDataFrame, col_a: str, col_b: str,
                       n: int = 3) -> np.ndarray:
    """
    Bin two columns into n tertile classes each and return the combined
    cell index (0..n*n-1) with index = n * bin_b + bin_a.
    """
    # qcut-style tertiles; use fixed breaks instead if classes must be stable
    # across multiple maps in an atlas.
    quantiles = np.linspace(0, 1, n + 1)
    breaks_a = np.quantile(gdf[col_a].to_numpy(), quantiles)
    breaks_b = np.quantile(gdf[col_b].to_numpy(), quantiles)
    # np.digitize returns 1..n on interior; clip the top edge back into 0..n-1.
    bin_a = np.clip(np.digitize(gdf[col_a], breaks_a[1:-1]), 0, n - 1)
    bin_b = np.clip(np.digitize(gdf[col_b], breaks_b[1:-1]), 0, n - 1)
    return n * bin_b + bin_a


def plot_bivariate(gdf: gpd.GeoDataFrame, col_a: str, col_b: str, n: int = 3):
    matrix = build_bivariate_matrix(n=n)
    cell = classify_bivariate(gdf, col_a, col_b, n=n)
    colors = matrix[cell]

    fig, ax = plt.subplots(figsize=(9, 9))
    gdf.plot(color=colors, ax=ax, edgecolor="0.85", linewidth=0.2)
    ax.set_axis_off()

    # Legend is the grid itself: a one-dimensional colour bar cannot describe
    # two variables. Inset a small axes and paint the 3x3 matrix.
    legend = fig.add_axes([0.70, 0.12, 0.22, 0.22])
    grid = matrix.reshape(n, n, 3)  # row 0 = low B at the bottom after origin flip
    legend.imshow(grid, origin="lower")
    legend.set_xlabel(col_a, fontsize=8)
    legend.set_ylabel(col_b, fontsize=8)
    legend.set_xticks([])
    legend.set_yticks([])
    return fig


# ---------------------------------------------------------------------------
# Example usage
# ---------------------------------------------------------------------------
# tracts = gpd.read_file("census_tracts.gpkg").to_crs("EPSG:5070")
#
# # Single variable: percentage change in population, anchored at zero change.
# fig1 = plot_diverging(tracts, "pop_change_pct", anchor=0.0)
# fig1.savefig("pop_change_diverging.png", dpi=300, bbox_inches="tight")
#
# # Two variables: income vs vacancy rate on a 3x3 bivariate grid.
# fig2 = plot_bivariate(tracts, col_a="median_income", col_b="vacancy_rate")
# fig2.savefig("income_vacancy_bivariate.png", dpi=300, bbox_inches="tight")

Performance Tuning and Cartographic Best Practices

  • Anchor the midpoint to meaning, not to the mean. TwoSlopeNorm(vcenter=0.0) for percentage change, vcenter=break_even for profit/loss, vcenter=national_average for relative comparison. Letting the centre drift to data.mean() makes half the features look “below” purely because of the sample, and the map changes meaning every time the input changes. The choice of anchor is a classification decision as much as a colour one, and pairs directly with the break strategy discussed in Data-Driven Classification.
  • Verify luminance symmetry before trusting a diverging ramp. Convert both endpoints to CIELAB and compare L*. If the two arms differ by more than a few units, the darker side reads as more extreme at equal data magnitude. The CET-D maps are pre-balanced; hand-built ramps from two arbitrary hues almost never are. This is the same perceptual-uniformity requirement that governs the sequential case in Color Palette Generation for Thematic Maps.
  • Blend bivariate cells in linear-light RGB. Averaging two sRGB values directly darkens and muddies the mid-cells because sRGB is gamma-encoded. Gamma-decode to linear, average, then re-encode — the build_bivariate_matrix helper does exactly this, which keeps the centre cell a clean neutral rather than a dull grey.
  • Cap the bivariate grid at 3x3 for most output. Nine colours is the practical ceiling for a legend a reader can decode cell by cell. Reserve 4x4 for large-format print only, and never push both variables past three classes on a screen map.
  • Precompute the matrix once, index many times. build_bivariate_matrix is cheap but constant across an atlas; build it once and reuse it so every sheet shares an identical legend. Feature colour assignment is then a single matrix[cell] fancy-index over the whole GeoDataFrame, which is vectorised and scales to hundreds of thousands of features without a Python loop.

Integration and Next Steps

Both halves emit standard matplotlib figures and per-feature RGB arrays, so they slot into existing render pipelines without special handling:

  • Static print export. The plot_diverging and plot_bivariate figures save at 300 DPI directly; the anchored TwoSlopeNorm and the fixed matrix guarantee the legend is reproducible across a multi-sheet run.
  • Dynamic legends. The 3x3 matrix and its tertile breaks are the exact inputs a legend generator needs. Feed build_bivariate_matrix output and the breaks_a / breaks_b arrays into the workflow described in Dynamic Legend Generation to emit the two-axis swatch grid automatically rather than drawing it by hand.
  • Web tiles. Serialise the per-feature hex colours as a style property and reference them from your tile renderer, so placement and classification are resolved upstream and the client only paints pre-computed colours.

For validation, always run the finished palette through a color-vision simulation pass: a red/blue diverging map is generally safe for deuteranopia, but a red/green bivariate blend is not. Pin the two hues after checking them, not before.

Frequently Asked Questions

Why does TwoSlopeNorm distort my colorbar tick spacing?

TwoSlopeNorm maps vmin..vcenter to the lower half of the colormap and vcenter..vmax to the upper half independently. If the data is asymmetric — say -20 to +80 around a vcenter of 0 — each data unit occupies a different amount of colour space on each side, so the colorbar ticks appear non-linearly spaced. That is by design: it keeps the neutral colour pinned to the anchor. If you need linear ticks, use a symmetric range with vmin=-vmax, accepting that the shorter arm wastes colour range.

Why does my custom diverging map look brighter on one side?

Most hand-picked ramps are not luminance-symmetric. A saturated blue and a saturated red can differ by 15-25 CIELAB L* units, so the darker arm reads as “more extreme” even at equal data magnitude. Match the two endpoints’ L* values in CIELAB, or use a pre-balanced map such as colorcet CET-D1 whose arms reach equal lightness contrast against the neutral centre.

How many classes should each variable have in a bivariate map?

Three per variable — a 3x3 grid, nine colours — is the practical maximum, because readers must distinguish every cell from its four neighbours. A 4x4 grid (sixteen colours) is legible only on large-format output with a generous legend. If you need finer resolution on one variable, keep the other at two classes (a 2x4 grid) rather than pushing both to four.

Back to Color Theory for GIS