Rendering Rasters with Matplotlib

Most of what goes wrong when a satellite raster is plotted comes down to three defaults that are reasonable for a photograph and wrong for geospatial data: the axes are in pixels, the colour limits come from the array’s extremes, and the image is interpolated between samples. Fixing those three takes four arguments, and the result is a figure that can be trusted. This topic is part of Visualization, Tiling & Web Delivery.

Beyond the basics, the interesting decisions are about honesty: which stretch to apply, whether a colour ramp is implying structure the data does not have, and how to render a nodata region so nobody mistakes it for a value. Those decisions carry through to every other delivery path, because a tile server is doing the same arithmetic on a 256-pixel window.

---

Prerequisites

pip install "matplotlib>=3.8" "rasterio>=1.3.0" "numpy>=1.23"

You need a raster whose nodata value is set correctly — the check is in extracting nodata and dtype from a GeoTIFF — and ideally one with overviews, so the decimated read in step one is cheap rather than merely correct.

---

Step-by-Step Workflow

Step 1 — Read only as many pixels as the figure can show

A figure has a fixed pixel budget: width in inches times dpi. Reading more than that is pure waste, and on a large scene it is the difference between a figure that renders instantly and one that exhausts memory.

import numpy as np
import rasterio


def read_for_figure(path: str, band: int = 1, *, target_px: int = 1600):
    """Read a decimated, masked band plus the map extent for plotting."""
    with rasterio.open(path) as src:
        factor = max(1, int(max(src.width, src.height) / target_px))
        arr = src.read(
            band,
            out_shape=(max(1, src.height // factor), max(1, src.width // factor)),
            masked=True,                       # nodata becomes a masked array
        )
        b = src.bounds
    return arr, (b.left, b.right, b.bottom, b.top)

masked=True is doing two jobs. It keeps the fill value out of the percentile calculation, and it makes matplotlib draw those cells as transparent rather than as the bottom of the colour ramp. Without it, a scene with a -9999 fill renders as a uniform block with one bright corner.

Step 2 — Derive display limits from the data

import numpy as np


def percentile_limits(arr: np.ma.MaskedArray, lower: float = 2.0,
                      upper: float = 98.0) -> tuple[float, float]:
    """Display limits from valid pixels only."""
    valid = arr.compressed() if np.ma.isMaskedArray(arr) else arr[np.isfinite(arr)]
    if valid.size == 0:
        raise ValueError("no valid pixels to derive a stretch from")
    return float(np.percentile(valid, lower)), float(np.percentile(valid, upper))

The choice of percentiles is a judgement about what to sacrifice. Two and ninety-eight clips a little detail at each end and gives good contrast across the body of the distribution; one and ninety-nine preserves more extremes at the cost of contrast. For a product where the extremes are the point — a burn severity map, say — use fixed limits chosen from the physics rather than percentiles from this particular scene, or two scenes of the same fire will be rendered on different scales.

Step 3 — Draw in map coordinates

import matplotlib.pyplot as plt


def plot_raster(path: str, band: int = 1, *, cmap: str = "viridis",
                title: str | None = None):
    arr, extent = read_for_figure(path, band)
    vmin, vmax = percentile_limits(arr)

    fig, ax = plt.subplots(figsize=(8, 8))
    im = ax.imshow(arr, extent=extent, cmap=cmap, vmin=vmin, vmax=vmax,
                   interpolation="nearest")           # never invent intermediate values
    ax.set_xlabel("Easting (m)")
    ax.set_ylabel("Northing (m)")
    ax.ticklabel_format(style="plain")                # no 1e6 offset on the axis
    if title:
        ax.set_title(title)
    cb = fig.colorbar(im, ax=ax, shrink=0.75, pad=0.02)
    cb.set_label("reflectance")
    return fig, ax

ticklabel_format(style="plain") is a small fix with a large payoff: by default matplotlib factors out a common offset and prints +3.9e5 under the axis, which turns readable UTM coordinates into a puzzle.

Step 4 — True-colour composites

A true-colour image is three bands stacked, each stretched on its own statistics. Stretching them together preserves the raw colour balance, which for satellite reflectance is dominated by atmospheric scattering in the blue.

import numpy as np
import rasterio


def true_colour(path: str, bands: tuple[int, int, int] = (3, 2, 1),
                *, target_px: int = 1600, lower: float = 2.0, upper: float = 98.0):
    """Stack red, green and blue with an independent stretch per band."""
    with rasterio.open(path) as src:
        factor = max(1, int(max(src.width, src.height) / target_px))
        shape = (max(1, src.height // factor), max(1, src.width // factor))
        planes = []
        for b in bands:
            a = src.read(b, out_shape=shape, masked=True).astype("float32")
            lo, hi = percentile_limits(a, lower, upper)
            planes.append(np.clip((a - lo) / max(hi - lo, 1e-6), 0, 1))
        extent = (src.bounds.left, src.bounds.right,
                  src.bounds.bottom, src.bounds.top)
    rgb = np.ma.dstack(planes)
    return rgb, extent

For false-colour composites the same function serves — pass near infrared, red and green as the band triple and vegetation renders red, which is the standard way to make crop structure legible. The band conventions differ per sensor, and getting them from band descriptions rather than positions avoids the usual mix-up, as discussed in reading and writing GDAL tags and band descriptions.

Step 5 — Export for the destination

fig.savefig("scene_ndvi.png", dpi=200, bbox_inches="tight", pad_inches=0.05)
fig.savefig("scene_ndvi.pdf", bbox_inches="tight")     # vector axes, raster image

A PDF keeps the axes, labels and colour bar as vectors while the image stays a raster, which is exactly what a publication wants. For a slide deck, a PNG at 150–200 dpi is enough; going higher only increases the file size, since the underlying array has no more detail to give.

---

Rendering Categorical Data

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import BoundaryNorm, ListedColormap
from matplotlib.patches import Patch

CLASSES = {1: ("forest", "#267300"), 2: ("cropland", "#a8cc54"),
           3: ("built", "#c4281b"), 4: ("water", "#1c5ca8")}


def plot_classes(arr, extent):
    codes = sorted(CLASSES)
    cmap = ListedColormap([CLASSES[c][1] for c in codes])
    norm = BoundaryNorm([codes[0] - 0.5] + [c + 0.5 for c in codes], cmap.N)

    fig, ax = plt.subplots(figsize=(8, 8))
    ax.imshow(np.ma.masked_equal(arr, 0), extent=extent, cmap=cmap, norm=norm,
              interpolation="nearest")
    ax.legend(handles=[Patch(facecolor=CLASSES[c][1], label=CLASSES[c][0])
                       for c in codes],
              loc="lower right", framealpha=0.9)
    return fig, ax

BoundaryNorm with half-integer edges is what pins each class code to exactly one colour; without it, a four-class raster drawn with a four-colour map still interpolates at the boundaries when the figure is resampled. Masking class 0 rather than colouring it keeps unclassified area transparent, which is almost always what you want over a basemap.

---

Diverging Data and the Zero Problem

Change images, anomalies and differences all share a property that trips up naive rendering: the sign matters as much as the magnitude. A field that gained vegetation and one that lost it are opposite results, and a sequential ramp renders them as merely different shades.

import matplotlib.pyplot as plt
from matplotlib.colors import TwoSlopeNorm

lo, hi = percentile_limits(diff, 2, 98)
norm = TwoSlopeNorm(vmin=min(lo, -1e-6), vcenter=0.0, vmax=max(hi, 1e-6))

fig, ax = plt.subplots(figsize=(8, 8))
im = ax.imshow(diff, extent=extent, cmap="RdBu", norm=norm, interpolation="nearest")
fig.colorbar(im, ax=ax, shrink=0.75).set_label("NDVI change")

TwoSlopeNorm maps the negative side and the positive side with different slopes so the centre lands exactly on zero. The colour bar then shows an uneven spacing, which looks odd the first time and is exactly right: the data really does extend further in one direction, and hiding that by symmetrising the range would either clip real change or waste half the ramp on values that never occur.

The alternative — forcing a symmetric range around zero — is defensible when several images must be compared on one scale, because then every map shares a legend. Choose deliberately, state which you used, and never let the choice vary between figures in the same document. The underlying arithmetic these figures render is covered in computing NDVI difference between two dates.

---

Parameter Reference

---

Annotation That Earns Its Place

A figure destined for anything other than a quick check needs three annotations, and rarely more: a scale bar, a north arrow where the orientation is not obvious, and a note of the source and date. Each answers a question a reader will otherwise ask.

from matplotlib.patches import Rectangle


def add_scale_bar(ax, length_m: float = 5000, *, label: str = "5 km") -> None:
    """A simple scale bar in data coordinates, anchored bottom-left."""
    x0, x1 = ax.get_xlim()
    y0, y1 = ax.get_ylim()
    pad_x, pad_y = 0.06 * (x1 - x0), 0.06 * (y1 - y0)
    bx, by = x0 + pad_x, y0 + pad_y

    ax.add_patch(Rectangle((bx, by), length_m, 0.008 * (y1 - y0),
                           facecolor="#0f3460", edgecolor="white", linewidth=0.8))
    ax.text(bx + length_m / 2, by + 0.018 * (y1 - y0), label,
            ha="center", va="bottom", fontsize=9, color="#0f3460")

The scale bar works in data coordinates, which is only meaningful because the axes are in map units — another payoff from passing extent in step three. In a geographic CRS the same code would be wrong, since a degree is not a fixed distance; reproject to a projected CRS before drawing anything with a scale bar on it.

What to leave out matters as much. Gridlines over imagery obscure the data and rarely help; a colour bar on a true-colour composite means nothing, since there is no single quantity being shown; and a title repeating what the caption already says wastes the space a legend could use. The discipline is the same one that applies to the data itself: every element should answer a question someone actually has.

---

Verification & Testing

The useful checks are visual, but three of them can be automated.

import numpy as np

arr, extent = read_for_figure("scene.tif")
assert extent[0] < extent[1] and extent[2] < extent[3], "extent order wrong"

valid = arr.compressed()
assert valid.size > 0.05 * arr.size, "over 95% of the scene is nodata"

vmin, vmax = percentile_limits(arr)
assert vmax > vmin, "degenerate stretch: the band may be constant"
spread = (vmax - vmin) / max(abs(vmax), 1e-6)
assert spread > 0.01, "stretch is extremely narrow — check the nodata value"

The narrow-stretch assertion catches the commonest silent failure: a fill value that was not declared, so it survives into compressed() and drags one percentile to an extreme. If it fires, print the unique values at the tails and the culprit will be obvious.

For a human check, render a quicklook of every product in a batch and view them as a contact sheet, as described on the section overview. Defects that no assertion catches — an inverted class map, a diagonal seam, a scene that is half cloud — are instantly visible in a grid of thumbnails.

---

Troubleshooting

The image is upside down

extent was passed with top and bottom swapped, or origin="lower" was set without adjusting the extent. Take the order straight from src.bounds as (left, right, bottom, top) and leave origin alone.

The whole figure is one flat colour

Either the band is genuinely constant, or a fill value is dominating the stretch. Read with masked=True and check arr.count() against arr.size.

Axis labels show +3.9e5

Matplotlib factored out an offset. Call ax.ticklabel_format(style="plain"), and consider dividing coordinates by 1000 and labelling in kilometres for a cleaner axis.

The figure takes minutes to save

The full-resolution array reached imshow. Read with out_shape sized to the figure; a 10,980-pixel scene in an 8-inch figure is 99% wasted work.

Colours differ between two scenes of the same product

Percentile limits were computed per scene. For a comparable series, fix the limits once — from a representative scene or from the physics — and reuse them everywhere, exactly as a tile service must.

---

Frequently Asked Questions

Q: Why is my raster upside down in matplotlib? Because imshow puts row zero at the top while a north-up raster has its highest northing in row zero. Passing extent with bottom and top in the right order fixes both the orientation and the axis labels; setting origin=‘lower’ without swapping the extent flips the image without fixing the coordinates.

Q: How do I make a true-colour composite that does not look washed out? Stretch each band independently on its own percentiles, then stack. A single stretch applied to all three bands preserves the colour balance of the raw data, which for satellite reflectance is heavily blue-biased by atmospheric scattering.

Q: Should I plot the full resolution raster? Never. A figure eight inches wide at 150 dpi shows 1,200 pixels, so reading a 10,980 pixel scene wastes almost all the work and the memory. Read with out_shape sized to the figure and let the overviews do the decimation.

Q: Can I overlay vector data on the same axes? Yes, and it works cleanly as long as both are in the same CRS. Reproject the vector layer to the raster’s CRS first, then plot it onto the same axes — the extent is already in map coordinates, so the geometries land in the right place without further work.

---

Related

• Plotting a GeoTIFF with Correct Extent and Axes — the extent and orientation question in full.
• Building True Colour Composites with Contrast Stretch — per-band stretching and colour balance.
• Choosing Colour Ramps for Continuous Raster Data — sequential, diverging and the ramps to avoid.
• Exporting Publication Quality Map Figures — scale bars, north arrows, dpi and vector export.
• Serving Raster Tiles with TiTiler — the same rendering decisions applied per tile.