Choosing a Pixel-Selection Rule for Composites

The rule is a statement of what the composite is for. Four common choices in one place:

median = stack.median("time")                                  # robust, mixes dates
ndvi = (stack.nir - stack.red) / (stack.nir + stack.red)
maxndvi = stack.isel(time=ndvi.fillna(-1).argmax("time"))     # peak greenness
first = stack.bfill("time").isel(time=0)                       # one real observation
p25 = stack.quantile(0.25, "time")                             # darker, suppresses haze

Every compositing rule answers the same question — which value represents this pixel for this period? — and they answer it differently. This page belongs to seamless mosaicking and edge blending in Satellite Processing Workflows & Index Pipelines.

What Each Rule Keeps and Loses

Five rules, three properties Median is robust to residual cloud but mixes bands from different dates, so band ratios can be inconsistent. Medoid picks a real observation close to the median, keeping bands coherent. Maximum NDVI captures peak greenness but favours cloud shadow edges and wet soil can be missed. First-valid keeps one real observation per pixel but inherits any unmasked cloud. Low percentiles suppress haze but darken the result. Compositing rules compared rule cloud robust bands coherent one date median high partial no medoid high yes yes max NDVI high yes yes first-valid low yes yes 25th percentile high partial no Max NDVI is coherent but biased toward the greenest date, which is a feature only if you want it.

Band coherence is the property most often overlooked. A per-band median takes the red value from one date and the near-infrared from another, so indices computed on a median composite are not the index of any real observation. For display this rarely matters; for classification using band ratios it can.

Environment & Setup

Package Version pin Used for
xarray >=2023.1 Reductions along time
numpy >=1.23 Medoid distance computation
dask >=2023.1 Lazy reductions over large stacks 
pip install "xarray>=2023.1" "numpy>=1.23" "dask>=2023.1"

Complete Working Example

import numpy as np
import xarray as xr


def select_by_index(stack: xr.Dataset, idx: xr.DataArray) -> xr.Dataset:
    return stack.isel(time=idx)


def composite(stack: xr.Dataset, rule: str = "medoid") -> xr.Dataset:
    """stack: masked reflectance with dims (time, y, x); cloud already set to NaN."""
    bands = [b for b in stack.data_vars]
    if rule == "median":
        return stack.median("time", skipna=True)
    if rule == "p25":
        return stack.quantile(0.25, "time", skipna=True).drop_vars("quantile")
    if rule == "first":
        return stack.bfill("time").isel(time=0)
    if rule == "maxndvi":
        ndvi = (stack.nir - stack.red) / (stack.nir + stack.red)
        return select_by_index(stack, ndvi.fillna(-2).argmax("time"))
    if rule == "medoid":
        arr = stack[bands].to_array("band")                     # band, time, y, x
        med = arr.median("time", skipna=True)
        dist = ((arr - med) ** 2).sum("band", skipna=False)     # NaN where any band missing
        best = dist.fillna(np.inf).argmin("time")
        return select_by_index(stack, best)
    raise ValueError(rule)

The medoid is the most useful rule that most pipelines do not use: it picks, for each pixel, the real observation closest to the multi-band median. The result is as robust to outliers as the median but keeps every band from the same acquisition, so indices remain physically meaningful. Its cost is one extra pass over the stack to compute distances.

Matching the Rule to the Purpose

Pick the rule from the use For a visual base map, median or medoid give clean results. For land-cover classification using band ratios, medoid keeps bands coherent. For crop or vegetation peak mapping, maximum NDVI captures the growing-season peak. For analysis that needs a single known date, first-valid in a sensible order is required. Purpose → rule base map median classification medoid vegetation peak max NDVI known date first-valid Store the rule in the output's metadata; two composites built differently are not comparable.

Maximum NDVI deserves a caution. It picks the greenest observation, which in croplands is the growing-season peak — often exactly what is wanted — but over water, snow and bare soil the “greenest” date can be a slightly contaminated one where cloud edges raise NIR. Combining max-NDVI with a good mask, as in dilating cloud masks to catch thin cirrus, keeps it honest.

Carrying the Date Along

For any rule that selects a real observation — medoid, max-NDVI, first-valid — keep the index of the selected time slice as an extra band. A “date of observation” layer lets analysts see where the composite jumps between seasons, explains visible seams, and allows date-dependent corrections such as BRDF normalisation to be applied after the fact. It costs one small integer band and saves a great deal of confusion when two neighbouring fields look different for reasons that turn out to be a six-week gap between their source images.

Cost of Each Rule at Scale

The rules differ widely in compute. Median and percentiles need the whole time series of each pixel in memory at once, which is why the time dimension should be a single chunk; on a large stack they dominate runtime. First-valid is the cheapest, since it can stop at the first clear observation and never touches later ones. Max-NDVI needs two bands to compute the index and then one gather per band, which is cheap. The medoid is the most expensive — a median, then a distance for every observation, then a gather — typically two to three times the cost of a plain median.

For continental composites that cost difference is real money, and it is reasonable to use different rules for different products: a median for the quick-look base map that is rebuilt often, and a medoid for the annual classification input that is built once and used many times.

Combining Rules

Rules can also be chained. A common robust recipe is to discard the brightest and darkest quarter of observations per pixel — removing most residual cloud and shadow — and then take the medoid of what remains. Another is to use max-NDVI within the growing season and median outside it. Whatever the recipe, apply it identically across every tile and year so the outputs stay comparable, and store the recipe as metadata rather than in a notebook.

Verification

Coherent composites reproduce their own indices For a coherent rule such as medoid, NDVI computed from the composite bands equals NDVI of the selected source observation, so points fall on the one-to-one line. For a per-band median, the points scatter around the line because bands come from different dates. Band coherence check medoid — on the line median — scattered 
import numpy as np

c = composite(stack, "medoid")
ndvi_c = (c.nir - c.red) / (c.nir + c.red)
ndvi_stack = (stack.nir - stack.red) / (stack.nir + stack.red)
# medoid output must equal an actual observation's NDVI at every pixel
diff = np.abs(ndvi_stack - ndvi_c).min("time")
assert float(diff.max()) < 1e-6

Common Errors

argmax fails on all-NaN pixels

Fully clouded pixels have no valid time. Fill NaN with a sentinel before argmax and mask those pixels afterwards.

Median composite looks hazy

Thin cloud survives masking and pulls the median up. Use a lower percentile or improve the mask.

Max-NDVI composite has speckle over water

Water NDVI is near zero on every date, so small noise decides the winner. Mask water or use medoid there.

Composites from different years disagree

They were built with different rules or date windows. Record the rule and window with each output.

Frequently Asked Questions

Q: Which compositing rule is best? None universally. Median for clean display, medoid for classification with coherent bands, max-NDVI for vegetation peaks, first-valid when you need a known acquisition date per pixel.

Q: What is a medoid composite? For each pixel, the actual observation whose band values are closest to the multi-band median. It is robust like a median but every band comes from the same date, so band ratios stay meaningful.

Q: Why do indices differ between median composites and single scenes? Because a per-band median takes each band from a potentially different date, so the ratio of medians is not the ratio of any real observation.

Q: Should I keep the source date? Yes, for any rule that selects real observations. A date band explains seams and enables later corrections.