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 127 Sentinel-2 scenes around Snowdon, an unmasked median NDVI composite gave 0.044 and the masked one gave 0.667 — a mean absolute difference of 0.587 per pixel (the masking guide measures it).

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 0, median 14, max 31
never observed clear: 0.09%
fewer than 3 observations: 1.15%

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

Holes, grey patches and stripes, each diagnosed

1. Holes: check whether the pixel was ever clear

Over 127 scenes across seven months, the median pixel in the measured window had 14 clear looks and 904 pixels — 0.09% — had none. Those 904 are holes in the data, not in the code; a hole anywhere else is a bug.

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 on the same stack, the masked mean composite gave an NDVI of 0.543 against the masked median's 0.667 — the 0.12 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.

On the 22 September 2025 scene, dilating the SCL mask on its 20 m grid left 42.7% usable with none, 40.0% with one pixel, 37.7% with two and 35.7% with three (the masking guide). Two pixels cost five percentage points of coverage and remove 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.044 against a masked 0.667. Because the median scene over the Snowdon window is 0.2% 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 3 observations and one built from 31 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.

A composite diagnostic

The data used on this page

Every Sentinel-2 L2A scene in Earth Search over a 10 km window around Snowdon (EPSG:32630, 422,880–432,880 E, 5,875,440–5,885,440 N) from 1 June to 31 December 2025 — 127 scenes. The red, near-infrared and SCL windows are cached per scene, shared with the composite guide. stack is per-scene NDVI with unusable pixels set to NaN, ndvi_all the same without the mask, masks the usable pixels and brightness the mean of red and near-infrared.

import os

import numpy as np
import pystac_client
import rasterio
from pyproj import Transformer
from rasterio.transform import from_origin
from rasterio.windows import from_bounds

BOUNDS = (422_880, 5_875_440, 432_880, 5_885_440)
to_wgs84 = Transformer.from_crs(32630, 4326, always_xy=True)
aoi = [*to_wgs84.transform(*BOUNDS[:2]), *to_wgs84.transform(*BOUNDS[2:])]
items = sorted(pystac_client.Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=aoi, datetime="2025-06-01/2025-12-31").items(), key=lambda i: i.id)


def window(href):
    with rasterio.open(href) as src:
        return src.read(1, window=from_bounds(*BOUNDS, transform=src.transform))


os.makedirs("scenes", exist_ok=True)
ndvi_all, masks, brightness, dates = [], [], [], []
for item in items:
    path = f"scenes/{item.id}.npz"
    if not os.path.exists(path):
        scl = window(item.assets["scl"].href).repeat(2, axis=0).repeat(2, axis=1)
        np.savez_compressed(path, red=window(item.assets["red"].href), nir=window(item.assets["nir"].href),
                            scl=scl, date=str(item.datetime.date()),
                            cloud_cover=item.properties["eo:cloud_cover"])
    z = np.load(path)
    red, nir = z["red"].astype("float32") / 10_000, z["nir"].astype("float32") / 10_000
    with np.errstate(all="ignore"):
        ndvi_all.append((nir - red) / (nir + red))
    masks.append(np.isin(z["scl"], [4, 5, 6, 7, 11]))
    brightness.append((red + nir) / 2)
    dates.append(str(z["date"]))
ndvi_all, masks, brightness = np.stack(ndvi_all), np.stack(masks), np.stack(brightness)
stack = np.where(masks, ndvi_all, np.nan)
profile = {"driver": "GTiff", "width": 1000, "height": 1000, "crs": "EPSG:32630",
           "transform": from_origin(BOUNDS[0], BOUNDS[3], 10, 10)}
print(stack.shape, dates[0], "to", dates[-1])

(127, 1000, 1000) 2025-06-04 to 2025-12-29

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


with np.errstate(all="ignore"):
    composite = np.nanmedian(stack, axis=0)
problems = diagnose_composite(ndvi_all, masks, composite)

  observations per pixel: min 0, median 14, max 31
  holes 0.09% (0.09% never observed, 0.00% unexplained)
  median |mean - median| across pixels: 0.0653
  masked vs unmasked median: 0.5929
  ! mean and median differ substantially — residual contamination the median is ignoring

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. On this stack the typical pixel's gap was 0.065 NDVI — over the 0.05 line, so the diagnostic flags it — while every one of the 0.09% of holes was a pixel never observed clear, so none is a bug.

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 its pixels kept and bright")

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


bright, per_scene = find_mask_failures(brightness, masks)
for i in np.argsort(-per_scene)[:5]:
    print(f"    scene {i} is {items[i].id}")

  0.899% of kept observations are brighter than 0.35
    scene 27: 14.25% of its pixels kept and bright
    scene 69: 13.82% of its pixels kept and bright
    scene 119: 12.38% of its pixels kept and bright
    scene 35: 9.10% of its pixels kept and bright
    scene 77: 8.18% of its pixels kept and bright
  56.77% of pixels have at least one bright surviving observation
    scene 27 is S2A_30UVD_20251015_1_L2A
    scene 69 is S2B_30UVD_20251015_0_L2A
    scene 119 is S2C_30UVD_20251119_0_L2A
    scene 35 is S2A_30UVD_20251124_0_L2A
    scene 77 is S2B_30UVD_20251124_1_L2A

Ranking by scene is what turns "the mask is leaking" into an actionable finding. Here five scenes on three dates stood out. On 15 October — both satellites — 28% of the window was labelled bare soil with a median red reflectance of 0.306, against 0.125 for bare soil on the clear 12 July: cloud the classifier called ground, on a scene that advertised 88% cloud. On 19 and 24 November the bright survivors were snow, class 11, which the keep list includes on purpose. The two need opposite fixes — drop the October scenes, keep or separate the snow — and only the per-scene ranking tells them apart.

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


composite_with_provenance(stack, masks, dates, profile, "composite.tif")
with rasterio.open("composite.tif") as ds:
    print([ds.descriptions[b] for b in range(3)], ds.tags()["min_obs"])

  127 scenes, observations 0-31, 1.15% below the minimum
['median composite', 'clear observations per pixel', 'index of the contributing scene'] 3

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.

Why unmasked composites collapse

Why an unmasked composite collapses

The measured numbers are stark: an unmasked median NDVI composite over 127 scenes gave 0.044, and the masked one 0.667.

The reason is that cloud is not a minority of the observations. With a median usable fraction of 0.2% 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.543 against the masked median's 0.667.

That 0.12 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 three or 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 five percentage points of coverage — from 42.7% to 37.7% on the measured scene — and remove the mixtures. With a median of 14 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 3 observations and one from 31 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 across the measured window the support ran from no clear looks to 31.

Holes with observations and dilation

• 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.2% 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.12 NDVI on measured data.

Why is my NDVI composite so low everywhere?

The mask was probably not applied. An unmasked median composite over 127 Snowdon scenes gave 0.044 against 0.667 masked.

Should I use mean or median?

Median. Both were computed over the same clear observations and the mean sat 0.12 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 98.9% of pixels had three or more, the median pixel 14, and 904 pixels none.

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.