Visualization, Tiling & Web Delivery

A raster nobody can see is a raster nobody will use. This section covers the last mile: turning arrays of reflectance, indices and class codes into figures for a report, layers on an interactive map, and tile services that let someone pan across a continent without downloading it. The techniques are unglamorous and the failures are conspicuous — a scene that renders black, a colour ramp that invents boundaries, a map that stalls because the file has no overviews.

It builds directly on the file structure described in Core Raster Fundamentals & STAC Mapping, because almost everything here works well or badly depending on whether the underlying file is a properly tiled Cloud-Optimized GeoTIFF with overviews. Products generated by Satellite Processing Workflows & Index Pipelines and Raster Machine Learning & Model Inference are the usual subjects.

---

Two Delivery Paths from One File

There are only two things a viewer can be given: pixels that were rendered in advance, or an address from which pixels are rendered on demand. Everything else is a variation on those.

The second sentence of that last box is the one that causes arguments. A figure rendered with matplotlib defaults and a tile service rendered with its own defaults will show visibly different images of the same data, and the person comparing them will reasonably conclude that one of them is wrong. Fixing the rendering parameters once — as data, not as code in two places — is the cheapest thing in this whole section.

---

Key Components

---

Production Patterns

Pattern 1 — A percentile stretch computed once and reused

Nearly every “my image renders black” problem is a stretch problem. Reflectance stored as scaled integers occupies perhaps a twentieth of the int16 range, so a viewer that stretches across the full range shows almost nothing.

import numpy as np
import rasterio


def percentile_limits(path: str, band: int = 1, *,
                      lower: float = 2.0, upper: float = 98.0,
                      max_pixels: int = 4_000_000) -> tuple[float, float]:
    """Display limits from an overview, so the scan is cheap on a large scene."""
    with rasterio.open(path) as src:
        # Pick an overview level that keeps the sample under max_pixels
        factor = max(1, int(np.sqrt(src.width * src.height / max_pixels)))
        arr = src.read(band, out_shape=(src.height // factor, src.width // factor),
                       masked=True)
    valid = arr.compressed()
    return float(np.percentile(valid, lower)), float(np.percentile(valid, upper))

Reading an overview rather than the full band is what makes this practical: the limits are statistically indistinguishable and the read is a thousand times cheaper. Store the returned pair alongside the product so every renderer uses the same numbers — the detail is covered in plotting a GeoTIFF with correct extent and axes.

Pattern 2 — Rendering in map coordinates, not pixel coordinates

An image drawn with imshow and no extent has axes in pixels, which is useless for anyone trying to relate it to a map. Passing the dataset bounds fixes it in one argument.

import matplotlib.pyplot as plt
import rasterio


def plot_band(path: str, band: int = 1, cmap: str = "viridis"):
    with rasterio.open(path) as src:
        arr = src.read(band, masked=True)
        left, bottom, right, top = src.bounds
        crs = src.crs

    vmin, vmax = percentile_limits(path, band)
    fig, ax = plt.subplots(figsize=(8, 8))
    im = ax.imshow(arr, extent=(left, right, bottom, top),
                   cmap=cmap, vmin=vmin, vmax=vmax, interpolation="nearest")
    ax.set_xlabel(f"Easting ({crs.linear_units})")
    ax.set_ylabel(f"Northing ({crs.linear_units})")
    fig.colorbar(im, ax=ax, shrink=0.75)
    return fig, ax

interpolation="nearest" matters more than it looks. The default smooths between pixels, which is actively misleading for a class raster and subtly misleading for anything else, because it invents values at the display resolution.

Pattern 3 — Serving tiles rather than files

Once a product is a valid COG, a tile server can read the overview level that matches each requested zoom and return a 256-pixel PNG. Nothing is pre-generated and nothing is duplicated.

from rio_tiler.io import Reader

with Reader("https://example-bucket.s3.amazonaws.com/products/class.tif") as cog:
    # x, y, z are the usual XYZ tile indices
    img = cog.tile(x=1203, y=1544, z=12, tilesize=256)
    png = img.render(img_format="PNG", colormap=cog.colormap)

The whole approach depends on the file being read efficiently over HTTP, which is exactly what reading a COG over S3 without downloading describes. A striped GeoTIFF served this way reads full-width rows for every 256-pixel tile and will bankrupt both the latency budget and the egress bill.

---

Choosing a Colour Scheme

Three rules cover almost every case. Sequential for quantities that run from low to high, with the light end at low so the printed version behaves. Diverging for quantities with a meaningful zero, and then the zero must be pinned to the neutral colour — a diverging ramp whose midpoint drifts with the data range is worse than a sequential one, because it implies a sign change that is not there. Categorical for classes, with hues chosen to be distinguishable rather than pretty, and with the palette stored in the file so every renderer agrees.

The greyscale test is the fastest check available: convert the rendered image to luminance and see whether the structure is still legible. A perceptually ordered ramp survives it; a rainbow ramp dissolves into bands, which is precisely what it does for a viewer with a colour vision deficiency.

---

How a Tile Request Becomes Pixels

Understanding what happens between a map pane scrolling and a PNG arriving explains most of the performance advice in this section, and nearly all of the failures.

Step three deserves the most attention. The server computes the resolution the tile needs and chooses the overview whose resolution is closest without going under; if the file has no overviews there is only one level to choose, and a whole-country view reads full-resolution pixels for the entire visible extent. That single fact explains most reports of a “slow tile server” — the server is fine, the file is not.

Step five is the one people forget to budget for. A UTM product served to a Web Mercator map is reprojected on every request. The warp itself is fast, but it forces the reader to fetch a slightly larger source window than the tile, and near UTM zone edges it can mean reading from a region with heavy distortion. For a layer that will be viewed constantly, storing a Web Mercator copy alongside the analytical UTM one is usually worth the duplication; for a layer viewed occasionally it is not.

Caching applies at two levels and they behave differently. The tile server’s own cache holds decoded arrays and helps when many tiles come from the same source window, which is exactly what happens when someone pans. A CDN in front holds encoded PNGs and helps when many people look at the same place, which is what happens with a published product. Both are worth having, and neither substitutes for overviews.

---

Notebooks Are a Delivery Channel Too

It is easy to treat notebook visualization as scratch work, but for most remote sensing teams the notebook is the interface — it is where analysts check whether a mask is behaving, whether a composite has seams, whether a classification looks plausible before anyone computes a metric.

That has two practical consequences. The first is that a notebook map should point at the same tile service the web map uses rather than loading arrays into the browser: a leafmap or folium layer consuming an XYZ endpoint shows exactly what a stakeholder will see, at any zoom, without a memory limit. Loading a decimated array into the notebook shows something else, and the difference is where “it looked fine in my notebook” comes from.

The second is that interactive inspection is a debugging tool with no real substitute. Clicking a pixel and reading its value across bands and dates finds problems — an unscaled band, a shifted mask, a missing date — in seconds, where a summary statistic would hide them. Those workflows are covered in interactive raster exploration in Jupyter, and they pair naturally with the time-series work in temporal aggregation and time series analysis.

A split-screen comparison deserves a special mention because it answers the question stakeholders actually ask, which is not “what does the classification say” but “what changed”. Two synchronised panes showing before and after, sharing a stretch, communicate more in five seconds than a difference histogram does in a meeting — and the underlying arithmetic is the same as in change detection and differencing workflows.

---

Quicklooks as a Quality Gate

A quicklook is a small rendered image generated automatically for every product a pipeline emits. It costs a second of compute and it catches an entire class of defect that no assertion will.

import matplotlib
matplotlib.use("Agg")                     # no display on a batch worker

import matplotlib.pyplot as plt
import rasterio


def write_quicklook(src_path: str, png_path: str, *, width: int = 1024) -> None:
    """A decimated, stretched PNG thumbnail for pipeline output review."""
    with rasterio.open(src_path) as src:
        factor = max(1, src.width // width)
        arr = src.read(1, out_shape=(src.height // factor, src.width // factor),
                       masked=True)
        bounds = src.bounds

    vmin, vmax = percentile_limits(src_path)
    fig, ax = plt.subplots(figsize=(6, 6), dpi=110)
    ax.imshow(arr, extent=(bounds.left, bounds.right, bounds.bottom, bounds.top),
              vmin=vmin, vmax=vmax, cmap="viridis", interpolation="nearest")
    ax.set_axis_off()
    fig.savefig(png_path, bbox_inches="tight", pad_inches=0)
    plt.close(fig)

Reviewing a contact sheet of a few hundred quicklooks takes a couple of minutes and reliably surfaces the scene that came out black, the tile that is half nodata, the composite with a diagonal seam, and the classification that inverted its classes. None of those are caught by a schema check, and all of them are obvious to a human glancing at a grid of thumbnails.

Wiring quicklook generation into the pipeline as a mandatory final step — and publishing the PNGs next to the products — also makes the archive browsable without any tooling at all, which matters more than it sounds for anyone who inherits the data later.

---

Common Pitfalls & Failure Modes

Everything renders black or white. The display limits come from the dtype range rather than the data. Compute a percentile stretch and pass vmin/vmax explicitly, every time.

Nodata renders as a real value. A fill value of −9999 pulls the low end of the stretch down to −9999 and everything real is squeezed into the top pixel of the ramp. Read with masked=True so the fill is excluded from both the statistics and the render.

The map stalls at low zoom. The file has no overviews, so drawing a whole-country view reads the full-resolution data. Build overviews, as in adding internal overviews with the right resampling.

Class colours shift between viewers. The palette lives in each viewer’s configuration rather than in the file. Write a colour table into the raster and the problem disappears for every consumer at once.

A figure and a web map disagree. Two renderers with two sets of parameters. Store the stretch, ramp and nodata handling as data next to the product and have both read it.

Tiles are served in the wrong projection. Web maps expect Web Mercator, and a tile server asked to serve a UTM file reprojects per request. That is correct but slow; for a heavily used layer, store a Web Mercator copy as described in generating Web Mercator tile pyramids.

---

Performance and Scale

The performance of a visualization stack is almost entirely a property of the source file. A tile request that lands on a matching overview level reads a handful of internal tiles — tens of kilobytes — and renders in single-digit milliseconds. The same request against a file without overviews reads the full resolution for that area, and at low zoom that means the whole scene.

Three numbers are worth knowing. A 256-pixel tile from a well-formed COG costs roughly two to four HTTP range requests, which is why latency rather than bandwidth dominates and why co-locating the tile server with the storage matters more than either. A pyramid of overviews adds about a third to file size, which is the cheapest performance purchase available. And a rendered PNG tile is typically 5–30 KB, so a CDN in front of the tile server converts a compute cost into a cache hit for every viewer after the first.

For static figures the constraint is different: the bottleneck is reading more pixels than the figure can show. A 10,980-pixel scene rendered into an 8-inch figure at 150 dpi has 1,200 pixels across, so reading the full resolution wastes 99% of the work. Use out_shape to decimate at read time, which is exactly what the overview machinery exists for.

When many products must be rendered — a quicklook per scene across an archive — the job is embarrassingly parallel and belongs on the patterns in Cloud Execution & Orchestration. One scene per task, output written straight to object storage, no shared state.

---

What Each Topic Covers

The four topics below follow the order a product moves from a working array to something other people look at.

Rendering rasters with matplotlib is the foundation: correct extents, honest stretches, true-colour composites, colour ramps that suit the data, and figures that survive being printed. Everything else in the section borrows its rendering decisions from here, because a tile server is doing the same arithmetic on a smaller window.

Serving raster tiles with TiTiler turns a COG archive into an endpoint. It covers deployment, the rescale and colormap parameters that make a tile URL self-describing, caching strategy, and serving a whole archive as one layer through MosaicJSON rather than merging files.

Interactive raster exploration in Jupyter is where analysts actually work: displaying a COG on a map widget, drawing an area of interest and reading pixels from it, comparing two dates side by side, and clicking through to a time series at a single location.

Preparing rasters for the web handles the conversions that dynamic tiling cannot do for free — Web Mercator pyramids for heavily viewed layers, 8-bit conversion with a documented stretch, thumbnails, and publishing static tiles when a fixed basemap really is the right answer.

A practical reading order for someone starting from nothing: get the stretch right first, because everything downstream inherits it; make the file a valid COG second, because every delivery path depends on it; and only then choose between a figure, a tile service and a notebook map, since by that point all three are cheap.

---

Frequently Asked Questions

Q: Do I need to pre-generate tiles or can they be dynamic? Dynamic tiling from a Cloud-Optimized GeoTIFF is the default answer now. It removes the pyramid build, keeps one authoritative file, and lets the rendering parameters change without regenerating anything. Pre-generated tiles still win for a fixed basemap served at very high request rates.

Q: Why does my satellite image render almost black? Because reflectance occupies a small part of the int16 range and the viewer is stretching across the full range. Compute a percentile stretch — typically the 2nd to 98th percentile of actual valid pixels — and pass those limits explicitly.

Q: Which colour ramp should I use for NDVI? A sequential ramp if you care about magnitude, or a diverging ramp centred on zero if the sign matters, as it does for a difference image. Avoid rainbow ramps: they create boundaries in the data that do not exist and are unreadable to a substantial fraction of viewers.

Q: Is it worth serving 16-bit data or should everything be 8-bit for display? Serve the 16-bit product as the data and render 8-bit for display. Converting the archive to 8-bit throws away precision that analysis needs, while rendering to 8-bit per request costs almost nothing and keeps the stretch adjustable.

Q: How do I keep a legend consistent between a figure and a web map? Store the legend as data rather than as code. For a class raster that means the colour table written into the file itself, which both matplotlib and a tile server can read. For a continuous product it means a small JSON sidecar holding the stretch limits and the ramp name, published alongside the raster and read by every renderer. The moment a colour decision exists only inside one script, the figure in the report and the layer on the map begin to drift apart, and nobody notices until someone puts them on the same slide.

---

Related

• Rendering Rasters with Matplotlib — figures with correct extents, stretches and colour bars.
• Serving Raster Tiles with TiTiler — dynamic tiling, rescaling parameters, mosaics and caching.
• Interactive Raster Exploration in Jupyter — notebook maps, AOI drawing and pixel inspection.
• Preparing Rasters for the Web — Web Mercator pyramids, 8-bit conversion and quicklooks.
• Understanding Cloud-Optimized GeoTIFF Structure — why the file layout decides how well all of this performs.