Color Theory for GIS: Automated Cartographic Design & Export Workflows

In automated cartographic pipelines, color is not an aesthetic preference — it is a precision data encoding channel. Every assignment of hue, value, or saturation to a spatial attribute carries implicit claims about data magnitude, categorical distinction, or bipolar deviation. When those assignments are made ad hoc, they introduce perceptual distortion, fail accessibility standards, and produce rendering inconsistencies across web tiles, print layouts, and agency deliverables. A reproducible, code-driven color workflow eliminates those failure modes by anchoring every decision in perceptual mathematics, validated against the same contrast and profile requirements that govern professional cartographic publishing.

This page establishes that workflow end-to-end: from statistical classification through palette selection, color space transformation, contrast validation, and ICC-tagged export. It builds on the architectural patterns in Automated Cartographic Design Fundamentals and connects to adjacent concerns — projection integrity and scale-dependent layout — that affect how color is perceived in the final output.

Color Encoding Pipeline

The diagram below shows the five-stage flow from raw attribute values to a publication-ready raster with embedded ICC profile. Each stage has a discrete responsibility; coupling any two stages causes the kind of rendering drift that is hard to trace in production.

Five-stage color encoding pipeline A horizontal flow diagram showing five stages: Classify & Normalize, Palette Selection, Color Space Transform, Contrast Validation, and ICC Export. Arrows connect each stage left to right. Classify & Normalize Palette Selection Color Space Transform Contrast Validation ICC-Tagged Export ① mapclassify ② colorcet / palettable ③ colorspacious ④ WCAG 2.2 ⑤ Pillow / cairosvg

Prerequisites and Environment Configuration

The stack below covers all five stages. Pin versions to prevent silent palette behavior changes between releases — matplotlib in particular has changed default colormap normalization behavior across minor versions.

python==3.11
geopandas==0.14.4
mapclassify==2.6.1
colorspacious==1.1.2
colorcet==3.1.0
palettable==3.3.3
matplotlib==3.8.4
numpy==1.26.4
Pillow==10.3.0
cairosvg==2.7.1

CRS requirements: Color assignment is attribute-driven, not projection-driven, but the downstream visual density of a choropleth is projection-dependent. Ensure your GeoDataFrame is in the intended display CRS before classification — equal-area projections prevent area-weighted visual bias in choropleth rendering. Consult Projection Selection Algorithms for programmatic CRS selection logic.

Conceptual Foundation: Perceptual Color Spaces

The engineering core of this workflow is the distinction between device-dependent and device-independent color spaces.

sRGB is device-dependent: the same (R, G, B) triplet renders differently on uncalibrated hardware. CIE XYZ and CIELAB are device-independent: they encode color as human perception experiences it, independent of the display. CAM02-UCS (Color Appearance Model 2002, Uniform Color Space) extends CIELAB to account for chromatic adaptation, surround luminance, and other perceptual non-linearities.

The practical consequence for GIS: sequential colormaps must have uniform perceptual step sizes, not uniform numerical step sizes. The jet colormap has near-zero perceptual contrast in its green band, which creates false visual boundaries in elevation or temperature rasters. viridis, plasma, and cividis were designed to be perceptually monotonic in CAM02-UCS space and are the correct defaults for continuous spatial data.

For categorical data, hue separation must also be measured in a perceptual space — palette entries that look distinct on a calibrated monitor may collapse under deuteranopia simulation. The Color Palette Generation for Thematic Maps page covers algorithmic hue maximization for categorical palettes.

Step-by-Step Implementation

Step 1 — Statistical Classification and Value Normalization

Raw attribute values must be binned before they can be mapped to a colormap index. The choice of classification method is not cosmetic — it changes which data features are emphasized and how legend breaks will be labeled.

  • Natural breaks (Jenks): minimizes within-class variance; best for data with natural clustering (income, population density).
  • Quantile: equal counts per class; prevents empty classes but can mask outliers.
  • Equal interval: constant bin width; appropriate only when the value range is meaningful in absolute terms.
  • Standard deviation: highlights deviation from the mean; correct for normalized anomaly data.
import numpy as np
import geopandas as gpd
import mapclassify

gdf = gpd.read_file("county_income.gpkg")

# Natural breaks into 5 classes — explicitly specify k to prevent
# mapclassify from choosing k automatically based on n.
classifier = mapclassify.NaturalBreaks(gdf["median_income"], k=5)
gdf["color_class"] = classifier.yb  # integer 0..k-1

# Normalize class index to [0, 1] for direct colormap lookup.
# Subtract 1 from k in the denominator so class k-1 maps to exactly 1.0.
gdf["normalized"] = gdf["color_class"] / (classifier.k - 1)

For diverging data anchored at a meaningful midpoint (e.g., income deviation from the national median), pin the midpoint to 0.5 before normalizing:

midpoint = 65_000  # national median, USD
gdf["delta"] = gdf["median_income"] - midpoint
max_abs = gdf["delta"].abs().max()
# Maps [-max_abs, 0, max_abs] → [0.0, 0.5, 1.0]
gdf["normalized"] = (gdf["delta"] / (2 * max_abs)) + 0.5

When this map will appear across multiple zoom levels, coordinate classification thresholds with Scale Mapping for Web and Print to verify that class boundaries remain visually meaningful at tile zoom 8 through 14.

Step 2 — Perceptually Uniform Palette Selection

Load palettes from validated libraries. Never construct sequential palettes by hand-picking hex values without measuring perceptual linearity.

import colorcet as cc
import numpy as np

# colorcet.b_linear_blue_95_50_c20 is a perceptually linear sequential palette.
# Access it as a list of 256 sRGB hex strings.
palette_hex = cc.b_linear_blue_95_50_c20  # 256 entries

def hex_to_rgb_float(h):
    h = h.lstrip("#")
    return tuple(int(h[i:i+2], 16) / 255.0 for i in (0, 2, 4))

palette_rgb = np.array([hex_to_rgb_float(h) for h in palette_hex])

To validate perceptual linearity, compute pairwise distances in CAM02-UCS space. The minimum inter-step distance should be roughly uniform across the palette — a large spike or dip indicates a perceptual discontinuity.

import colorspacious as cspace

# Convert the full palette from sRGB to CAM02-UCS
ucs_palette = cspace.convert(palette_rgb, "sRGB1", "CAM02-UCS")

# Compute sequential step distances (not pairwise — just adjacent steps)
step_distances = np.linalg.norm(np.diff(ucs_palette, axis=0), axis=1)
print(f"Min step: {step_distances.min():.4f}  Max step: {step_distances.max():.4f}")
# A well-designed palette will have a small max/min ratio (< 3x).

Step 3 — Color Space Transformation and Gamma Correction

Apply gamma encoding before any raster write. Matplotlib’s colormaps return linear float RGB; writing those values directly to PNG produces a washed-out render on gamma-corrected displays.

import numpy as np
import matplotlib.pyplot as plt

def linear_to_srgb(linear_rgb: np.ndarray) -> np.ndarray:
    """
    Apply the IEC 61966-2-1 sRGB transfer function to linear float RGB.
    Input: float array in [0, 1], shape (..., 3).
    Output: gamma-encoded float array in [0, 1].
    """
    out = np.where(
        linear_rgb <= 0.0031308,
        12.92 * linear_rgb,
        1.055 * np.power(np.clip(linear_rgb, 0, 1), 1.0 / 2.4) - 0.055,
    )
    return np.clip(out, 0.0, 1.0)

cmap = plt.get_cmap("viridis")
# cmap returns (N, 4) RGBA in linear space
linear_rgba = cmap(gdf["normalized"].values)
srgb_rgb = linear_to_srgb(linear_rgba[:, :3])  # drop alpha for now

Note the exponent 1/2.4, not 1/2.2. The sRGB standard (IEC 61966-2-1) specifies a two-piece function with exponent 2.4 in the power segment — 1/2.2 is an approximation that introduces measurable error in shadows.

Step 4 — Contrast Validation

Every foreground/background pair must be validated before rendering. The WCAG 2.2 relative luminance formula uses a linearization threshold of 0.04045 — the earlier 0.03928 value from WCAG 2.0 errata is incorrect and produces wrong results near that boundary.

import numpy as np

def relative_luminance(srgb_0_255: tuple) -> float:
    """
    Compute WCAG 2.2 relative luminance from an 8-bit sRGB tuple.
    """
    srgb = np.array(srgb_0_255, dtype=float) / 255.0
    linear = np.where(
        srgb <= 0.04045,
        srgb / 12.92,
        ((srgb + 0.055) / 1.055) ** 2.4,
    )
    return float(0.2126 * linear[0] + 0.7152 * linear[1] + 0.0722 * linear[2])

def contrast_ratio(fg: tuple, bg: tuple) -> float:
    l1 = relative_luminance(fg)
    l2 = relative_luminance(bg)
    lighter, darker = max(l1, l2), min(l1, l2)
    return (lighter + 0.05) / (darker + 0.05)

# Validate label color against each class fill color
label_color = (30, 30, 30)
for i, hex_val in enumerate(palette_hex[::51]):  # sample 5 classes from 256-entry palette
    fill = tuple(int(hex_val.lstrip("#")[j:j+2], 16) for j in (0, 2, 4))
    ratio = contrast_ratio(label_color, fill)
    status = "PASS" if ratio >= 4.5 else "FAIL"
    print(f"Class {i}: {ratio:.2f}:1 ({status})")

For full batch validation across all map layers, including transparency compositing and dynamic basemap backgrounds, see WCAG Contrast Checking for Map Layers, which covers raster compositing and the special case of semi-transparent fills over satellite imagery.

Step 5 — Export with ICC Profile Embedding

Raster outputs must carry an embedded ICC profile. Without it, browsers, OS color managers, and print RIPs will assume a default profile — often sRGB — but will not know whether the encoding is correct, leading to rendering inconsistencies across workflows.

import matplotlib.pyplot as plt
import geopandas as gpd
from PIL import Image
import io

fig, ax = plt.subplots(figsize=(12, 8), dpi=150)
gdf.assign(
    facecolor=list(map(tuple, (srgb_rgb * 255).astype(int)))
).plot(
    ax=ax,
    color=[tuple(c / 255 for c in row) for row in (srgb_rgb * 255).astype(int)],
    linewidth=0.3,
    edgecolor="white",
)
ax.set_axis_off()

# Write to in-memory buffer first; avoids a temp-file on disk.
buf = io.BytesIO()
fig.savefig(buf, format="png", dpi=300, bbox_inches="tight", transparent=False)
plt.close(fig)
buf.seek(0)

# Embed sRGB IEC61966-2-1 ICC profile.
icc_path = "/usr/share/color/icc/sRGB_v4_ICC_preference.icc"
with open(icc_path, "rb") as f:
    icc_data = f.read()

img = Image.open(buf)
img.save("county_income_map.png", "PNG", icc_profile=icc_data)

For SVG/PDF vector exports, declare the color space explicitly in the file header rather than relying on viewer defaults. cairosvg respects the color_mode parameter:

import cairosvg

# Convert a styled SVG to PDF/X-3 with embedded color profile
cairosvg.svg2pdf(
    url="styled_map.svg",
    write_to="county_income_map.pdf",
    dpi=300,
)

Complete Working Example

The function below wires all five stages into a single callable that accepts a GeoDataFrame, a numeric column name, a classification method string, and an output path. It returns the augmented GeoDataFrame with contrast metrics attached.

import io
import numpy as np
import geopandas as gpd
import mapclassify
import colorspacious as cspace
import colorcet as cc
import matplotlib.pyplot as plt
from PIL import Image


def automated_color_export(
    gdf: gpd.GeoDataFrame,
    value_col: str,
    output_path: str,
    k: int = 5,
    method: str = "NaturalBreaks",
    icc_path: str = "/usr/share/color/icc/sRGB_v4_ICC_preference.icc",
) -> gpd.GeoDataFrame:
    """
    Classify, assign perceptually uniform colors, validate contrast, and export
    a GeoDataFrame as a 300-DPI PNG with an embedded sRGB ICC profile.

    Parameters
    ----------
    gdf         : Input GeoDataFrame (must be in display CRS before calling).
    value_col   : Column name to map to color.
    output_path : Destination PNG file path.
    k           : Number of classification classes.
    method      : mapclassify classifier name.
    icc_path    : Path to sRGB ICC profile on disk.

    Returns
    -------
    gdf with added columns: color_class, normalized, srgb_hex, contrast_vs_white.
    """
    # --- Stage 1: Classify ---
    classifier = getattr(mapclassify, method)(gdf[value_col], k=k)
    gdf = gdf.copy()
    gdf["color_class"] = classifier.yb
    gdf["normalized"] = gdf["color_class"] / (classifier.k - 1)

    # --- Stage 2: Palette (perceptually linear blue sequential) ---
    palette_hex = cc.b_linear_blue_95_50_c20
    def hex_to_rgb_float(h):
        h = h.lstrip("#")
        return tuple(int(h[i:i+2], 16) / 255.0 for i in (0, 2, 4))
    palette_rgb = np.array([hex_to_rgb_float(h) for h in palette_hex])

    # --- Stage 3: Gamma encode ---
    def linear_to_srgb(arr):
        return np.where(
            arr <= 0.0031308,
            12.92 * arr,
            1.055 * np.power(np.clip(arr, 0, 1), 1.0 / 2.4) - 0.055,
        )

    cmap = plt.get_cmap("cividis")
    linear_rgba = cmap(gdf["normalized"].values)
    srgb_rgb = np.clip(linear_to_srgb(linear_rgba[:, :3]), 0, 1)
    gdf["srgb_hex"] = [
        "#{:02x}{:02x}{:02x}".format(*row)
        for row in (srgb_rgb * 255).astype(np.uint8)
    ]

    # --- Stage 4: Contrast validation (labels vs white background) ---
    def rel_lum(srgb_01):
        lin = np.where(srgb_01 <= 0.04045, srgb_01 / 12.92,
                       ((srgb_01 + 0.055) / 1.055) ** 2.4)
        return 0.2126 * lin[0] + 0.7152 * lin[1] + 0.0722 * lin[2]

    white_lum = 1.0
    gdf["contrast_vs_white"] = [
        (white_lum + 0.05) / (rel_lum(row) + 0.05)
        for row in srgb_rgb
    ]

    # --- Stage 5: Render and export with ICC ---
    fig, ax = plt.subplots(figsize=(12, 8))
    gdf.plot(
        ax=ax,
        color=list(map(tuple, srgb_rgb)),
        linewidth=0.3,
        edgecolor="white",
    )
    ax.set_axis_off()

    buf = io.BytesIO()
    fig.savefig(buf, format="png", dpi=300, bbox_inches="tight")
    plt.close(fig)
    buf.seek(0)

    with open(icc_path, "rb") as f:
        icc_data = f.read()

    img = Image.open(buf)
    img.save(output_path, "PNG", icc_profile=icc_data)
    print(f"Exported: {output_path}  (contrast range {gdf['contrast_vs_white'].min():.2f}–{gdf['contrast_vs_white'].max():.2f}:1)")
    return gdf

Performance Optimization Patterns

Vectorize colormap lookups. Calling cmap() once with the full normalized array is O(n); calling it row-by-row is O(n·overhead). The working example above always passes the full gdf["normalized"].values array in one call.

Cache ICC file reads. open(icc_path, "rb").read() on every export in a batch loop adds measurable I/O overhead. Load icc_data once outside the loop and pass it in — at 500KB per profile, a 200-map batch saves 100MB of reads.

Reuse the classifier for incremental data. mapclassify classifiers are not thread-safe, but they are stateless after fitting. Fit once on the full dataset, then apply classifier.yb (or classifier(new_series)) to subsets in parallel batch workers.

Pre-validate contrast before rendering. Computing contrast ratios on the normalized float array (Stage 4) takes microseconds. Rendering a full 300-DPI figure (Stage 5) takes seconds. Check that all classes pass before committing to the render — fail fast to avoid wasted CPU time in CI.

Common Pitfalls and Debugging

Gamma applied twice. If your pipeline calls both linear_to_srgb() and Matplotlib’s internal gamma (which savefig() applies when format="png" and the colormap returns linear values), the output will appear too dark with flattened shadows. Fix: either disable Matplotlib’s gamma management by working in already-encoded sRGB throughout, or apply gamma yourself and write raw bytes to Pillow directly without passing through savefig().

Wrong linearization threshold in WCAG formula. Using 0.03928 instead of 0.04045 produces a luminance error of roughly 0.3% near the threshold. For most colors this does not change the pass/fail outcome, but it will cause a borderline 4.5:1 pair to flip. Always use 0.04045 from WCAG 2.2.

mapclassify.NaturalBreaks fails with k >= n unique values. When a column has fewer unique values than requested classes, the Jenks algorithm raises a ValueError. Guard with k = min(k, gdf[value_col].nunique()) before fitting.

colorcet palette indexing produces banding. cc.b_linear_blue_95_50_c20 is a 256-entry list. Indexing it with a normalized float (e.g., palette_hex[gdf["normalized"]]) will raise a TypeError — you must convert to integer index first: int(norm * 255). Better: use Matplotlib’s LinearSegmentedColormap.from_list() to convert the colorcet palette to a proper callable colormap.

ICC profile path not found on headless build servers. /usr/share/color/icc/sRGB_v4_ICC_preference.icc is present on most Debian/Ubuntu installs with icc-profiles but absent on minimal CI images. Bundle the profile in your repository under assets/icc/sRGB_v4_ICC_preference.icc and reference it by relative path from the project root — this makes the export step reproducible across environments.

Choropleth color density changes after reprojection. If you run color assignment before reprojecting, visual polygon areas will shift after gdf.to_crs() is called, making the thematic encoding look different than intended. Always reproject to display CRS before classification. See Projection Selection Algorithms for the correct order of operations.

Conclusion

A rigorous color workflow treats every pipeline stage — classification, palette selection, space transformation, contrast validation, and ICC export — as a discrete, testable engineering concern. By anchoring palette choices in CAM02-UCS perceptual distances, enforcing the correct WCAG 2.2 linearization threshold, and embedding ICC profiles at export time, automated map pipelines produce output that is visually accurate, accessibility-compliant, and reproducible across rendering environments. The automated_color_export() function above is a drop-in starting point for any batch generation system; extend it by loading palette configuration from a YAML manifest so the same function can drive sequential, diverging, and qualitative outputs without code changes.

Back to Automated Cartographic Design Fundamentals