Histogram Matching Across Scenes

To make two adjacent scenes look continuous, fit a gain and offset from the pixels they share and apply it to one of them:

import numpy as np

# a and b are masked arrays over the overlap, cloud already excluded
gain = float(a.std() / b.std())
offset = float(a.mean() - gain * b.mean())

b_matched = b * gain + offset       # b now shares a's brightness and contrast

This is the radiometric half of the problem in Seamless Mosaicking and Edge Blending; feathering is the geometric half.

Why This Arises in Remote Sensing Workflows

Two scenes acquired on different days differ before any surface change is involved. The sun is at a different elevation, so the same slope reflects differently. The atmosphere holds different aerosol, so contrast differs. The sensor may be a different platform in the same constellation with slightly different calibration. Put those two scenes side by side and the join is visible as a brightness step, even when both are correct.

Feathering does not fix this. Blending across the overlap turns a step into a ramp, which is less obvious but still visible, and it is worse for analysis because the ramp is now a spatial gradient in the data. The fix has to be radiometric: bring the two scenes onto a common brightness scale first, then blend the residual.

The critical constraint is that the correction must be estimated from ground the two scenes share. Statistics computed over each scene’s full extent describe different landscapes — one may contain a city and the other a lake — so matching them makes the difference worse, not better.

Fitting the correction on shared ground Over the overlap, scene B's histogram is shifted brighter and slightly narrower than scene A's. Fitting a gain and offset on those shared pixels and applying it to B brings the two distributions together without changing the shape of either scene's land-cover signal. overlap histograms, before after gain + offset on B scene A scene B a visible step at the join distributions coincide; the join disappears Only pixels inside the overlap are used for the fit — statistics from the full scenes describe different landscapes. The correction is two numbers per band, which is small enough to record with the product.

Environment & Setup

Package Version Why
rasterio ≥1.3.0 Reading the overlap windows
numpy ≥1.23 Fitting and applying the correction
scikit-image ≥0.22 Optional: full histogram matching
shapely ≥2.0 Computing the overlap geometry 
pip install "rasterio>=1.3.0" "numpy>=1.23" "scikit-image>=0.22" "shapely>=2.0"

Complete Working Example

This function computes the overlap between two scenes, fits a per-band gain and offset on the shared valid pixels, and applies it — returning the coefficients so they can be recorded.

What to exclude from the overlap before fitting A gain and offset fitted on contaminated pixels compensates for weather rather than radiometry. Excluding cloud, shadow, water and saturated pixels from the overlap sample is what makes the fit describe the sensor difference rather than the scene content. Pixels that must not enter the fit exclude why effect if left in cloud and shadow not surface reflectance gain fitted to weather water near-zero NIR skews the spread implausible gain, often > 1.5 saturated pixels clipped at the top of the range compresses the fitted contrast nodata on either date not observed fit on undefined values changed ground not a radiometric difference the correction erases real change Fit on shared, clear, unchanged ground — everything else is a different quantity. 
import numpy as np
import rasterio
from rasterio.windows import from_bounds
from shapely.geometry import box


def overlap_windows(src_a, src_b):
    """Matching windows over the shared extent of two datasets, or None."""
    inter = box(*src_a.bounds).intersection(box(*src_b.bounds))
    if inter.is_empty or inter.area == 0:
        return None, None
    b = inter.bounds
    return (from_bounds(*b, transform=src_a.transform),
            from_bounds(*b, transform=src_b.transform))


def fit_gain_offset(
    path_reference: str,
    path_adjust: str,
    *,
    band: int = 1,
    valid_range: tuple[float, float] = (0.0, 1.0),
    sample: int = 512,
) -> dict:
    """Fit a linear radiometric correction for `path_adjust` onto `path_reference`."""
    with rasterio.open(path_reference) as ref, rasterio.open(path_adjust) as adj:
        win_ref, win_adj = overlap_windows(ref, adj)
        if win_ref is None:
            raise ValueError("scenes do not overlap — match to a common reference instead")

        a = ref.read(band, window=win_ref, out_shape=(1, sample, sample), masked=True)[0] \
            if False else ref.read(band, window=win_ref, out_shape=(sample, sample), masked=True)
        b = adj.read(band, window=win_adj, out_shape=(sample, sample), masked=True)

    good = (~a.mask) & (~b.mask)
    good &= (a >= valid_range[0]) & (a <= valid_range[1])
    good &= (b >= valid_range[0]) & (b <= valid_range[1])
    if good.sum() < 1000:
        raise ValueError(f"only {int(good.sum())} shared valid pixels — correction unreliable")

    av, bv = np.asarray(a[good], dtype="float64"), np.asarray(b[good], dtype="float64")
    gain = float(av.std() / bv.std())
    offset = float(av.mean() - gain * bv.mean())

    residual_before = float(np.mean(av - bv))
    residual_after = float(np.mean(av - (bv * gain + offset)))
    return {"band": band, "gain": round(gain, 5), "offset": round(offset, 6),
            "n_pixels": int(good.sum()),
            "mean_residual_before": round(residual_before, 5),
            "mean_residual_after": round(residual_after, 6)}


def apply_correction(path_in: str, path_out: str, coeffs: dict[int, dict]) -> None:
    """Apply per-band gain and offset, recording them in the output's tags."""
    with rasterio.open(path_in) as src:
        profile = src.profile | {"dtype": "float32", "driver": "COG", "compress": "DEFLATE"}
        with rasterio.open(path_out, "w", **profile) as dst:
            for band in range(1, src.count + 1):
                arr = src.read(band, masked=True).astype("float32")
                c = coeffs.get(band)
                if c:
                    arr = arr * np.float32(c["gain"]) + np.float32(c["offset"])
                dst.write(arr.filled(src.nodata if src.nodata is not None else 0), band)
                if c:
                    dst.update_tags(band, gain=c["gain"], offset=c["offset"])


if __name__ == "__main__":
    coeffs = {b: fit_gain_offset("scene_a.tif", "scene_b.tif", band=b) for b in (1, 2, 3, 4)}
    for b, c in coeffs.items():
        print(b, c["gain"], c["offset"], "residual:", c["mean_residual_before"], "→",
              c["mean_residual_after"])
    apply_correction("scene_b.tif", "scene_b_matched.tif", coeffs)

The residual before and after is the check that matters: a correction that does not reduce the mean difference over the overlap has not worked, usually because the overlap contains cloud or water that dominates the statistics.

Variant Patterns

1. Full histogram matching, and when it is too much

from skimage.exposure import match_histograms

b_matched = match_histograms(b_array, a_array)      # reshapes B's distribution onto A's

This forces B’s histogram to equal A’s exactly, which removes the seam completely — and also removes any genuine difference in land-cover composition between the two scenes. It is appropriate for visual products where continuity is the goal, and inappropriate for analytical products where the values must remain measurements.

2. Matching a chain of scenes

Chained corrections drift; anchored corrections do not Matching each scene to its neighbour in sequence compounds every fitting error, so the last scene in a long chain can drift far from the first. Matching every scene to a single anchor scene, or to a common coarse reference, keeps each error independent and bounded. Five scenes, two correction topologies chained A B C D E four fits compound: E drifts from A by the sum of the errors anchored anchor A B C D E Anchoring needs every scene to overlap the anchor — otherwise use a coarse common reference instead.

3. Matching without an overlap

When scenes do not touch, fit on pseudo-invariant features instead: surfaces whose reflectance should be identical in both, such as deep water, bare rock, or paved areas. Extract their statistics from each scene independently and fit the correction between those, rather than between the full distributions.

pif = (ndvi < 0.15) & (ndwi < -0.1)          # bare, non-water, non-vegetated
gain = float(a_values[pif_a].std() / b_values[pif_b].std())

The assumption is that the two scenes contain comparable amounts of such surface, which is worth checking before trusting the result.

Verifying the Correction

Two checks are enough, and both operate on the overlap.

The mean residual over shared pixels should fall by an order of magnitude after correction. If it barely moves, the fit was contaminated — usually by cloud in the overlap, or by water, whose near-zero reflectance in the near-infrared skews the standard deviation.

A transect across the join should be flat. Extract a line of pixels crossing the seam in the corrected mosaic and plot it: any remaining step is the residual the blending step will have to hide, and its size tells you how wide the feather needs to be, as discussed in Removing Seams in Multi-Scene Mosaics with Feathering.

It is also worth checking what the correction did outside the overlap. A gain far from 1.0 — beyond roughly 0.8 to 1.25 — usually indicates a bad fit rather than a genuinely different scene, and applying it will distort the whole raster to fix a strip.

Common Errors

The gain comes out enormous or negative

Water or cloud dominates the overlap sample. Mask them out and restrict the fit to a plausible reflectance range before fitting.

The join looks worse after matching

The correction was fitted on full-scene statistics rather than on the overlap, so it aligned two different landscapes. Always fit on shared pixels.

Indices change after matching

Expected, and a reason for care: a gain and offset applied to reflectance changes normalised indices unless the same correction is applied to both bands. Apply per-band corrections consistently, and prefer to compute indices before matching if the indices are the deliverable.

Frequently Asked Questions

Q: Linear gain-offset or full histogram matching? Linear first. A gain and offset fitted on the overlap fixes most illumination and atmospheric differences while preserving the shape of the distribution. Full histogram matching also reshapes the distribution, which can erase real land-cover differences between scenes.

Q: What if the scenes do not overlap? Match both to a common reference instead — a coarse-resolution product, or an earlier mosaic — using pseudo-invariant features such as bare soil and deep water that should look the same in both.

Q: Should I match before or after masking? After. Cloud and shadow pixels in the overlap will dominate the fitted statistics and produce a correction that compensates for weather rather than for radiometry.