Building True Colour Composites with Contrast Stretch

Stretch each band on its own percentiles, then stack:

import numpy as np


def stretch(band: np.ndarray, lo: float, hi: float) -> np.ndarray:
    return np.clip((band.astype("float32") - lo) / max(hi - lo, 1e-9), 0, 1)

rgb = np.dstack([stretch(red, *lim_r), stretch(green, *lim_g), stretch(blue, *lim_b)])

A single shared stretch preserves the raw colour balance of the data, which for satellite reflectance is heavily blue-biased by the atmosphere. This page belongs to rendering rasters with matplotlib in Visualization, Tiling & Web Delivery.

Why Per-Band Stretching Looks Right

Rayleigh scattering is far stronger at short wavelengths, so the blue band of an uncorrected scene carries a substantial additive haze that the red band does not. The three bands therefore occupy different value ranges for the same surface, and treating them identically renders that atmospheric difference as colour.

Three bands, three ranges The blue band's distribution sits noticeably higher than the green, which sits higher than the red, because atmospheric scattering adds a wavelength-dependent offset. A shared stretch therefore renders every surface too blue. Stretching each band onto the same output range removes the offset and restores natural colour. Reflectance distributions by band blue green red a shared stretch spans all three and keeps the offsets per-band stretch each band mapped onto the full 0-1 output range the haze offset cancels This is an empirical haze correction, not a physical one — it removes the appearance, not the signal. For quantitative work, correct properly instead of stretching.

It is worth being explicit that this is a display trick. The per-band stretch makes an image look natural by removing an additive offset empirically; it does not recover surface reflectance, and a composite built this way should never be measured. Proper correction is the subject of applying dark object subtraction in Python.

Environment & Setup

Package Version pin Used for
rasterio >=1.3.0 Reading bands by description, decimated and masked
numpy >=1.23 Percentiles, stretching, gamma and stacking
matplotlib >=3.8 Displaying and exporting the composite 
pip install "rasterio>=1.3.0" "numpy>=1.23" "matplotlib>=3.8"

Complete Working Example

import matplotlib.pyplot as plt
import numpy as np
import rasterio


def composite(path: str, band_names: tuple[str, str, str] = ("B04", "B03", "B02"),
              *, target_px: int = 1600, percentiles: tuple[float, float] = (2, 98),
              gamma: float = 0.75,
              fixed_limits: dict[str, tuple[float, float]] | None = None):
    """Build an RGB composite with a per-band stretch and gamma correction."""
    with rasterio.open(path) as src:
        names = list(src.descriptions)
        missing = [n for n in band_names if n not in names]
        if missing:
            raise ValueError(f"missing bands {missing}; file has {names}")

        factor = max(1, int(max(src.width, src.height) / target_px))
        shape = (max(1, src.height // factor), max(1, src.width // factor))

        planes, used = [], {}
        for name in band_names:
            arr = src.read(names.index(name) + 1, out_shape=shape,
                           masked=True).astype("float32")
            if fixed_limits and name in fixed_limits:
                lo, hi = fixed_limits[name]
            else:
                lo, hi = np.percentile(arr.compressed(), percentiles)
            used[name] = (float(lo), float(hi))
            scaled = np.clip((arr - lo) / max(hi - lo, 1e-9), 0, 1)
            planes.append(np.power(scaled, gamma))     # gamma brightens mid-tones

        extent = (src.bounds.left, src.bounds.right,
                  src.bounds.bottom, src.bounds.top)

    rgb = np.ma.dstack(planes)
    return rgb.filled(0), extent, used


if __name__ == "__main__":
    rgb, extent, limits = composite("S2A_36NYF_20260614_stack.tif")
    print("stretch used:", limits)          # record these for reuse across scenes

    fig, ax = plt.subplots(figsize=(9, 9))
    ax.imshow(rgb, extent=extent, interpolation="nearest")
    ax.ticklabel_format(style="plain")
    ax.set_aspect("equal")
    fig.savefig("truecolour.png", dpi=180, bbox_inches="tight")

Returning the limits that were used, and printing them, is the habit that makes a sequence of scenes comparable. Run the function once on a representative scene, take the limits it reports, and pass them as fixed_limits for everything else in the collection.

False Colour and Other Band Triples

Four band triples worth keeping Red green blue gives a natural image. Near infrared with red and green renders vegetation bright red and is the classic choice for crop structure. Shortwave infrared with near infrared and red separates burned ground and moisture. Shortwave infrared with near infrared and blue is the standard agriculture composite. Same function, different bands B04 B03 B02 true colour what a person would see for reports and context B08 B04 B03 colour infrared vegetation renders red crop structure, vigour B12 B08 B04 shortwave composite burns and bare soil stand out fire, drought, geology B11 B08 B02 agriculture composite crop vigour and moisture field-level monitoring Band numbers are Sentinel-2; the same wavelengths carry different numbers on Landsat. Select by description, never by position, or a cross-sensor script produces confident nonsense.

Because the function takes band names rather than indices, switching composites is a one-argument change:

rgb_cir, extent, _ = composite("stack.tif", band_names=("B08", "B04", "B03"))
rgb_swir, extent, _ = composite("stack.tif", band_names=("B12", "B08", "B04"))

The shortwave bands arrive at 20 m on Sentinel-2, so a composite mixing them with 10 m visible bands needs the resampling described in resampling Sentinel-2 20 m bands to 10 m before the arrays can be stacked at all.

Verification

import numpy as np

assert rgb.ndim == 3 and rgb.shape[2] == 3, rgb.shape
assert rgb.min() >= 0 and rgb.max() <= 1, "values outside the display range"

for i, name in enumerate(("red", "green", "blue")):
    mean = float(rgb[..., i].mean())
    print(f"{name:<6} mean {mean:.3f}")
assert 0.15 < rgb.mean() < 0.75, "composite is too dark or too bright overall"
Channel means as a balance check A correctly stretched composite has red, green and blue means within about a tenth of each other. An unbalanced one, where the per-band stretch was skipped, shows a blue mean far above the other two, which is exactly the blue cast a viewer sees. balanced shared stretch R G B spread 0.04 — natural colour R G B spread 0.29 — the blue cast

Comparing the three channel means is the quickest colour-balance check available. For a mixed land scene they should be within roughly 0.1 of each other; a blue mean far above the others means the per-band stretch was skipped, and a red mean far above the others usually means a colour-infrared triple was used by mistake.

The overall brightness bound catches the two failures that make a composite unusable: a scene dominated by cloud, which pushes every channel toward one, and a scene dominated by water or shadow, which pushes them toward zero. Neither is a bug in the code, and both are worth knowing about before the image reaches a report.

Common Errors

The composite is entirely white

The scene is mostly cloud, so the percentile limits sit in the bright tail and everything else saturates. Mask cloud before computing the stretch, using masking clouds with the Sentinel-2 SCL band.

Colours look inverted or lurid

The band order is wrong — most often blue, green, red instead of red, green, blue. Select by description and print the names actually used.

Nodata areas render as black rather than transparent

filled(0) was used to produce a plain array. Keep the masked array and pass it to imshow, or build an RGBA array with alpha from the mask.

Frequently Asked Questions

Q: Why does a raw satellite composite look blue and washed out? Atmospheric scattering raises the blue band far more than the red, so the raw ratio between bands is genuinely blue-biased. Stretching each band on its own percentiles removes that bias, which is why a per-band stretch produces a natural image and a shared stretch does not.

Q: What gamma value should I use? Around 0.7 to 0.8 for satellite reflectance, which brightens the mid-tones where most land surfaces sit. Below about 0.6 the image starts to look flat and washed out; above 1.0 it darkens, which is almost never what you want for a scene that is already dark.

Q: How do I keep composites comparable across dates? Fix the stretch limits once for the collection and reuse them, rather than recomputing percentiles per scene. Per-scene stretching makes every image look good on its own and makes any sequence of them flicker.

Q: Should the stretch be linear or histogram-equalised? Linear for anything that will be compared, because it is reproducible and invertible. Histogram equalisation maximises local contrast and is useful for spotting faint features in a single image, but it makes two scenes incomparable and exaggerates noise in flat regions.