Converting Rasters to 8-bit for Display

Map valid values onto 1–255 with a fixed stretch, and keep 0 for nodata:

import numpy as np

def to_uint8(arr: np.ma.MaskedArray, lo: float, hi: float, gamma: float = 1.0) -> np.ndarray:
    scaled = np.clip((arr.astype("float32") - lo) / (hi - lo), 0, 1) ** gamma
    return (scaled * 254 + 1).filled(0).astype("uint8")     # 0 = nodata, 1-255 = data

The display copy is a derivative. The analytical file is never modified and remains the only thing anyone measures. This page belongs to preparing rasters for the web in Visualization, Tiling & Web Delivery.

What the Conversion Throws Away

Clip, quantise, reserve A sixteen-bit reflectance range is reduced in three ways. Values below the lower stretch limit are clipped to level one and values above the upper limit to level 255. Everything between is quantised into 254 steps. Level zero is kept for nodata only. The result is a picture, not a measurement. Three lossy steps in one conversion int16 source: 0 to 10,000 reflectance units clipped low stretched: 300 to 3,500 into 254 steps clipped high 0 uint8 output: 1 to 255 for data nodata About 12 source units collapse into each display level — fine for looking, useless for measuring.

Each of the three losses is irreversible. Clipped values are gone; quantisation merges roughly a dozen source units into each display level; and gamma, if applied, redistributes those levels non-linearly. Recording the parameters makes it possible to say what a display value means approximately, but never to recover the original.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Reading masked bands and writing tags
numpy >=1.23 Stretch, gamma and quantisation 
pip install "rasterio>=1.3.0" "numpy>=1.23"

Complete Working Example

import json

import numpy as np
import rasterio
from rasterio.windows import Window


def collection_limits(paths: list[str], band: int = 1, *, sample_px: int = 2_000_000,
                      pct: tuple[float, float] = (2, 98)) -> tuple[float, float]:
    """One pair of limits for a whole collection, from decimated samples."""
    samples = []
    per_file = max(1, sample_px // max(len(paths), 1))
    for p in paths:
        with rasterio.open(p) as src:
            f = max(1, int(np.sqrt(src.width * src.height / per_file)))
            a = src.read(band, out_shape=(src.height // f, src.width // f), masked=True)
            samples.append(a.compressed())
    allv = np.concatenate(samples)
    return float(np.percentile(allv, pct[0])), float(np.percentile(allv, pct[1]))


def write_display_copy(src_path: str, dst_path: str,
                       limits: dict[int, tuple[float, float]],
                       *, gamma: float = 1.0, block: int = 1024) -> None:
    with rasterio.open(src_path) as src:
        profile = src.profile | {"dtype": "uint8", "nodata": 0, "compress": "deflate",
                                 "tiled": True, "blockxsize": 512, "blockysize": 512}
        profile.pop("photometric", None)
        with rasterio.open(dst_path, "w", **profile) as dst:
            for row in range(0, src.height, block):
                for col in range(0, src.width, block):
                    win = Window(col, row, min(block, src.width - col),
                                 min(block, src.height - row))
                    for b in range(1, src.count + 1):
                        lo, hi = limits[b]
                        a = src.read(b, window=win, masked=True).astype("float32")
                        s = np.clip((a - lo) / max(hi - lo, 1e-9), 0, 1) ** gamma
                        dst.write((s * 254 + 1).filled(0).astype("uint8"), b, window=win)
            for b, (lo, hi) in limits.items():
                dst.update_tags(b, stretch_min=f"{lo:g}", stretch_max=f"{hi:g}")
            dst.update_tags(display_copy="true", gamma=f"{gamma:g}",
                            source=src_path,
                            decode="value = lo + ((v - 1) / 254) ** (1/gamma) * (hi - lo)")

The windowed loop keeps memory flat regardless of scene size, and computing the collection limits from decimated samples keeps that step cheap too — percentile estimates from a few million pixels are indistinguishable from those over the full archive. Writing a decode expression into the tags is a courtesy to future readers: it says exactly what a display value approximately means.

Gamma and Where the Levels Go

Gamma reallocates the 254 display levels With gamma of one, display levels are spread evenly across the stretched range. With gamma of 0.7 the curve rises faster at the low end, allocating more of the 254 levels to dark values where most land surfaces sit and fewer to bright values such as cloud and bare sand. Stretched input against display level gamma 1.0 gamma 0.7 stretched input, 0 to 1 display level more levels for dark land, fewer for cloud

Satellite reflectance is concentrated at the dark end — vegetation, water and soil all sit well below the bright tail formed by cloud, snow and sand. A linear conversion spends half its levels on that bright tail, where there is little to see. A gamma around 0.7 moves those levels down into the range where the landscape actually is, and the image gains visible detail without any change to the stretch limits.

For a single-band index product rendered through a colour ramp, leave gamma at 1.0: the ramp already controls how values map to colours, and adding a gamma makes the colour bar lie about the value at each colour.

Keeping Display and Analysis Apart

The safest arrangement is to make it structurally impossible to confuse the two files. Name them so the difference is in the filename — ndvi_analysis.tif and ndvi_display.tif — and set a display_copy=true tag on the derivative, as the example does. A pipeline step that computes statistics can then refuse any input carrying that tag, which turns a convention into a guard.

It is also worth publishing them from different places: the analytical file under the product path that a catalog item points at, the display file under a path that only the tile service and the web viewer use. People find files through catalogs and links, and if the only link to the 8-bit copy is from a map, nobody will pick it up by mistake. The catalog side of that separation is covered in attaching model metadata to STAC items.

Verification

A healthy display histogram Display values spread across most of the 1 to 255 range with small spikes at level one and level 255 from clipped pixels. Large spikes at either end mean the stretch limits are too tight; a narrow cluster in the middle means they are too loose. Clip fractions should be small at both ends 1 255 Orange: clipped pixels. Aim for a few per cent at each end, not more. 
import numpy as np
import rasterio

with rasterio.open("ndvi_display.tif") as src:
    v = src.read(1)
    tags = src.tags(1)
valid = v[v > 0]
low_clip = float((valid == 1).mean())
high_clip = float((valid == 255).mean())
print(f"clipped low {low_clip:.1%}, high {high_clip:.1%}")
assert low_clip < 0.05 and high_clip < 0.05, "stretch limits are too tight"
assert "stretch_min" in tags and "stretch_max" in tags, "stretch not recorded"

The clip fractions are the direct measure of whether the limits suit the data. With 2nd and 98th percentile limits from the same collection they should sit near 2% at each end; much higher means the limits came from a different product or the scene is unusual, and much lower on both ends means the limits are loose and contrast is being wasted.

Common Errors

Nodata areas render as black

Zero was written but the file’s nodata was not set to 0. Set nodata=0 in the profile so viewers treat it as transparent.

Scene boundaries show in a mosaic of display copies

Limits were computed per scene. Compute one set for the collection and apply it everywhere.

The display copy has banding in smooth gradients

Quantisation to 254 levels is visible across a very wide stretch. Narrow the limits, or accept that 8-bit cannot show fine gradients across a large range.

Values near zero disappear

The valid data included zero and was not offset to start at 1. Map onto 1–255, never 0–255.

Frequently Asked Questions

Q: Why reserve zero for nodata? Because an 8-bit band has no spare value outside its range. If valid pixels can be zero, nodata and the darkest real value become indistinguishable. Mapping valid data onto 1 to 255 costs one grey level and keeps the footprint unambiguous.

Q: Can the 8-bit copy be used for analysis? No. It is quantised to 256 levels, clipped at both ends of the stretch, and possibly gamma-corrected, so arithmetic on it gives wrong answers. It is a picture of the data, and the analytical file remains the only thing to measure.

Q: Should every scene get its own stretch? Not if the scenes will ever be seen together. Per-scene stretches make each image look good and make a mosaic of them patchy. Choose limits from a representative sample of the collection and apply them everywhere.

Q: Is a dynamic tiler a substitute for a display copy? Usually, yes — a tiler applies the stretch per request, so the analytical file can be served directly. A display copy earns its place when the layer is viewed constantly, when static tiles are needed, or when the rendering must be frozen for a publication.