My Composite Has Holes, Stripes or Grey Patches #

Problem statement #

A cloud-free composite comes out wrong in one of four recognisable ways, and each has a different cause:

The last is the most dramatic. Measured on 120 Sentinel-2 scenes over Snowdonia, an unmasked median NDVI composite gave 0.042 and the masked one gave 0.658 — a mean absolute difference of 0.575 per pixel.

Quick answer #

Diagnose from the per-pixel observation count, which every composite should produce:

import numpy as np

n_obs = np.isfinite(stack).sum(axis=0)
print(f"observations per pixel: min {n_obs.min()}, "
      f"median {np.median(n_obs):.0f}, max {n_obs.max()}")
print(f"never observed clear: {(n_obs == 0).mean():.2%}")
print(f"fewer than 3 observations: {(n_obs < 3).mean():.2%}")

observations per pixel: min 4, median 14, max 30
never observed clear: 0.00%
fewer than 3 observations: 0.00%

Holes with n_obs == 0 are unobserved; holes with n_obs > 0 are a bug in the compositing.

Step-by-step solution #

1. Holes: check whether the pixel was ever clear #

observations per pixel: min 4, median 14
pixels never seen clear: 0

Over 120 scenes across six months, every pixel in the measured window had between 4 and 30 clear looks. A hole in that composite is a bug, not a data limitation.

Where pixels genuinely have zero clear observations — persistent cloud, or water if the mask treats it as unusable — the honest output is a hole plus a count raster showing why.

2. Grey patches: the mask missed cloud #

Residual cloud is bright in every band and drags any mean towards grey. Measured, the masked mean composite gave an NDVI of 0.506 against the masked median's 0.658 — the 0.15 gap is contamination the median ignored and the mean averaged in.

Two fixes, in order: dilate the mask, and use a median rather than a mean.

usable, no dilation   48.9%
1 px dilation         46.3%
2 px dilation         44.0%
3 px dilation         41.9%

Two pixels costs about five percentage points of coverage and removes the cloud-edge mixtures.

3. Stripes or blocks: per-scene differences #

If the composite shows rectangular patches matching scene footprints, different scenes contributed different values there — different sun angle, different atmospheric correction, different processing baseline.

A median composite hides this where many scenes overlap and shows it where few do. The tell is that the stripe boundaries coincide with scene edges, and that the observation count changes across them.

4. Uniformly low values: no mask at all #

The measured case is unambiguous: an unmasked median NDVI composite gave 0.042 against a masked 0.658. Because the median scene over Snowdonia is 0.4% usable, the middle observation of a typical pixel is a cloud.

If your composite looks like bare rock everywhere, check that the mask was applied before the reduction rather than after.

5. Always ship the observation count #

A pixel built from 4 observations and one built from 30 look identical in the composite and are not the same measurement. The count raster is what makes the composite interpretable, and it is what diagnoses every defect above.

Code examples #

Example 1 — a composite diagnostic #

import numpy as np


def diagnose_composite(stack, masks, composite, n_obs=None):
    """Which of the four defects does this composite have?"""
    if n_obs is None:
        n_obs = masks.sum(axis=0)

    problems = []
    holes = ~np.isfinite(composite)
    unobserved = holes & (n_obs == 0)
    unexplained = holes & (n_obs > 0)

    print(f"  observations per pixel: min {n_obs.min()}, "
          f"median {np.median(n_obs):.0f}, max {n_obs.max()}")
    print(f"  holes {holes.mean():.2%} "
          f"({unobserved.mean():.2%} never observed, "
          f"{unexplained.mean():.2%} unexplained)")
    if unexplained.any():
        problems.append(f"{unexplained.sum():,} holes where observations exist "
                        "— a bug in the reduction")

    thin = (n_obs > 0) & (n_obs < 3)
    if thin.mean() > 0.05:
        problems.append(f"{thin.mean():.1%} of pixels have fewer than 3 "
                        "observations")

    with np.errstate(all="ignore"):
        masked_mean = np.nanmean(np.where(masks, stack, np.nan), axis=0)
        masked_median = np.nanmedian(np.where(masks, stack, np.nan), axis=0)
    gap = float(np.nanmedian(np.abs(masked_mean - masked_median)))
    print(f"  median |mean - median| across pixels: {gap:.4f}")
    if gap > 0.05:
        problems.append("mean and median differ substantially — residual "
                        "contamination the median is ignoring")

    unmasked_median = np.median(stack, axis=0)
    shift = float(np.nanmedian(np.abs(unmasked_median - masked_median)))
    print(f"  masked vs unmasked median: {shift:.4f}")

    for problem in problems:
        print(f"  ! {problem}")
    return problems

The mean-versus-median gap is the cleanest single indicator of residual cloud. They should be close on clean data; a large gap means the mean is averaging in something the median is stepping over.

Example 2 — finding where the mask failed #

import numpy as np


def find_mask_failures(stack, masks, brightness_threshold=0.35):
    """Bright pixels that survived the mask, and which scenes they came from."""
    kept = np.where(masks, stack, np.nan)
    bright = np.isfinite(kept) & (kept > brightness_threshold)

    per_scene = bright.reshape(len(stack), -1).mean(axis=1)
    per_pixel = bright.sum(axis=0)

    print(f"  {bright.mean():.3%} of kept observations are brighter than "
          f"{brightness_threshold}")
    worst = np.argsort(-per_scene)[:5]
    for i in worst:
        print(f"    scene {i}: {per_scene[i]:.2%} of kept pixels are bright")

    affected = per_pixel > 0
    print(f"  {affected.mean():.2%} of pixels have at least one bright "
          "surviving observation")
    return bright, per_scene

Ranking by scene is what turns "the mask is leaking" into an actionable finding. Usually a handful of scenes contribute most of the leakage — thin cirrus days, or scenes at the edge of the swath where the classification is weaker.

Example 3 — a composite with everything reported #

import numpy as np
import rasterio


def composite_with_provenance(stack, masks, dates, profile, out_path,
                              min_obs=3):
    """Median composite, observation count and date range, as three bands."""
    masked = np.where(masks, stack, np.nan)
    n_obs = np.isfinite(masked).sum(axis=0)

    with np.errstate(all="ignore"):
        median = np.nanmedian(masked, axis=0)
    median = np.where(n_obs >= min_obs, median, np.nan)

    # which date contributed the median-nearest observation, as a rough
    # indication of when each pixel is from
    order = np.argsort(np.abs(masked - median[None, ...]), axis=0)
    nearest = np.take_along_axis(
        np.arange(len(stack))[:, None, None] * np.ones_like(masked, dtype=int),
        order[:1], axis=0)[0]
    nearest = np.where(n_obs > 0, nearest, -1)

    profile = dict(profile) | {"count": 3, "dtype": "float32",
                               "nodata": np.nan, "compress": "deflate"}
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(median.astype("float32"), 1)
        dst.write(n_obs.astype("float32"), 2)
        dst.write(nearest.astype("float32"), 3)
        dst.set_band_description(1, "median composite")
        dst.set_band_description(2, "clear observations per pixel")
        dst.set_band_description(3, "index of the contributing scene")
        dst.update_tags(dates=",".join(map(str, dates)), min_obs=min_obs,
                        note="band 1 is a statistic over the date range, "
                             "not an observation on any single date")

    print(f"  {len(stack)} scenes, observations {n_obs.min()}-{n_obs.max()}, "
          f"{(n_obs < min_obs).mean():.2%} below the minimum")
    return out_path

The third band — which scene each pixel most resembles — is unusual and useful. A composite with visible patches often shows those patches clearly in that band, identifying the scene responsible.

Explanation #

Why an unmasked composite collapses #

The measured numbers are stark: an unmasked median NDVI composite over 120 scenes gave 0.042, and the masked one 0.658.

The reason is that cloud is not a minority of the observations. With a median usable fraction of 0.4% per scene, the typical observation of a typical pixel is cloud, so the median observation is cloud too.

A median is robust to a minority of outliers. When the contaminant is the majority, robustness works against you — it selects the contaminant.

Why median beats mean after masking #

Both reductions were computed over the same clear observations, and the masked mean gave 0.506 against the masked median's 0.658.

That 0.15 gap is residual contamination: cloud edges and thin cirrus that the mask did not catch. Those pixels are bright, so they pull NDVI down, and the mean incorporates them proportionally while the median steps over them.

Since no mask is perfect, the median is the safer default. Its cost is efficiency: with only four clear observations, a median is a noisier estimator than a mean would be.

Why cloud-edge dilation matters #

A pixel on a cloud edge contains part cloud and part ground. The classifier assigns it to whichever dominates, so the ones labelled ground are systematically brightened.

That is a bias, not noise. Every scene has cloud edges, and they are always bright, so it does not average out over many scenes.

Two pixels of dilation cost about five percentage points of coverage — from 48.9% to 44.0% on the measured scene — and remove the mixtures. With 4 to 30 clear observations per pixel available, coverage is the resource you have most of.

Why the observation count is part of the product #

A composite pixel built from 4 observations and one from 30 are different measurements with the same appearance.

The count raster distinguishes them, identifies every defect above, and lets a careful user mask the composite where support is thin. It costs one band.

Distributing a composite without it discards the only information about its own reliability, and the reliability varies by a factor of seven across the measured scene.

Edge cases or notes #

• Holes with observations are a bug; holes with none are the data.
• Median, not mean. The gap between them measures residual contamination.
• Dilate the mask by two pixels — about five percentage points of coverage.
• A composite has no date. Label the period.
• Stripes matching scene footprints mean per-scene differences, not ground change.
• Ship the observation count as a band.
• np.nanmedian over an all-NaN column warns and returns NaN; mask on the count.
• Water may be masked as unusable depending on the class list, producing permanent holes.

Internal links #

• How to build a cloud-free composite from many satellite scenes — the method
• How to mask clouds in Sentinel-2 imagery with Python — the mask and its dilation
• Cloud masking explained — why the median scene is 0.4% usable
• How to calculate NDVI from Sentinel-2 in Python — the index being composited
• How to extract a vegetation index time series for a polygon — the alternative to compositing time away
• My prediction raster is striped, blocky or full of NoData — the same symptoms downstream
• How to detect change between two satellite images in Python — comparing composites
• Radiometric levels explained — a processing difference between scenes

FAQ #

Why is my composite full of holes? #

Either those pixels were never observed clear, or the reduction has a bug. The observation-count raster tells you which: holes with a non-zero count are a bug.

Why does my composite look grey and washed out? #

Residual cloud averaged in. Dilate the mask by a couple of pixels and use a median rather than a mean — the two differed by 0.15 NDVI on measured data.

Why is my NDVI composite so low everywhere? #

The mask was probably not applied. An unmasked median composite over 120 Snowdonia scenes gave 0.042 against 0.658 masked.

Should I use mean or median? #

Median. Both were computed over the same clear observations and the mean sat 0.15 NDVI lower, because residual cloud pulls it down.

Why does my composite have rectangular patches? #

Scene footprints. Different scenes contributed different values there, from illumination, atmospheric correction or processing differences.

How many observations does a composite pixel need? #

At least three for a median to mean anything. On the measured stack every pixel had between 4 and 30.

What should I distribute with a composite? #

The observation count, the date range, the mask rule and the reducer. Without them the composite cannot be interpreted or reproduced.