Creating Monthly NDVI Composites with xarray

Given a cloud-masked Sentinel-2 NDVI stack with an irregular time axis, collapse it into clean monthly median composites in one call: group the time dimension by calendar month with resample(time="1MS") and reduce with a NaN-aware .median():

import xarray as xr

# ndvi: DataArray with dims (time, y, x); cloud pixels already NaN
monthly = ndvi.resample(time="1MS").median()   # skipna=True by default

Each output pixel is the median of that month’s valid observations, and masked NaNs are skipped automatically. This page belongs to Temporal Aggregation and Time-Series Analysis and assumes you already have per-scene NDVI — produced with the arithmetic in Calculating NDVI directly from xarray DataArrays — and a clean cloud mask.

Why This Arises in Remote Sensing Workflows

Sentinel-2 revisits a location every five days at best, but clouds, sensor gaps, and orbit geometry mean the usable NDVI observations for a pixel arrive on an irregular calendar: three clear dates one month, none the next, seven the following. Downstream analyses — phenology curves, anomaly detection, dashboards — want a regular grid: one value per pixel per month. Temporal composites bridge that gap by reducing all the observations that fall inside a period to a single representative value.

xarray.resample is the tool because it understands the datetime coordinate. It bins the irregular time axis into fixed calendar periods and applies a reducer to each bin, mirroring the groupby semantics of pandas but generalised to N-dimensional arrays. You describe the period with a frequency string — 1MS for month starts — and xarray handles the arithmetic of assigning each acquisition to its bin, however uneven the spacing.

Two details make it correct for masked satellite data. First, cloud pixels must already be NaN, not zero or a fill value, so the reducer ignores them — xarray’s .median() skips NaNs by default. This is the single most common source of a wrong composite: a fill value of 0 stored where a cloud was masked reads as a real, very low NDVI observation and pulls the monthly median down, most visibly along coastlines and cloud edges. Second, the median rather than the mean resists any residual haze or thin cloud that survived masking, because one bright outlier cannot drag the median the way it drags an average. On a series that has already been carefully masked the two agree closely; the median simply fails more gracefully when the mask let something through.

This topic lives in Temporal Aggregation and Time-Series Analysis under Satellite Processing Workflows & Index Pipelines. The masking that produces the NaNs comes from Cloud and Shadow Masking Strategies, and when the stack is too large for one machine you run the same call on a Dask cluster per Computing NDVI with xarray on a Dask Cluster.

How resample Bins an Irregular Time Axis

The diagram shows irregular acquisition dates, some masked to NaN, being grouped into month-start bins and reduced by a NaN-aware median.

Irregular NDVI observations resampled to monthly medians Irregular acquisition dates, some cloud-masked to NaN, are grouped into month-start bins. Each bin is reduced by a NaN-aware median to one composite value per month. June July August valid NaN Jun medianof 2 valid Jul medianof 2 valid Aug medianof 2 valid resample(time="1MS").median() — NaN observations skipped

A month whose observations are all NaN returns NaN rather than a fabricated value, which is the correct signal that the month has no usable data — never a zero. This distinction matters as soon as the composites feed anything quantitative: a phenology fit or an anomaly baseline can be told to skip NaN months, but it will happily and wrongly consume a fabricated zero. Preserving genuine gaps as NaN keeps the honesty of the series intact all the way through the pipeline, and pairing the composite with the per-month observation count makes the density of support explicit for anyone reading the result later.

Environment & Setup

Package Minimum version Why required
xarray 2023.6 resample on a datetime coordinate and NaN-aware reducers
numpy 1.24 Backing arrays and finite-count logic
pandas 2.0 Datetime index parsing behind xarray’s time coordinate 
pip install "xarray>=2023.6" "numpy>=1.24" "pandas>=2.0"

Complete Working Example

The function takes a list of per-scene NDVI arrays with their acquisition dates, builds a time-indexed DataArray, resamples to monthly medians, and returns both the composite and a per-month valid-observation count.

import numpy as np
import pandas as pd
import xarray as xr

def monthly_ndvi_composite(
    ndvi_scenes: list[np.ndarray],
    dates: list[str],
    y: np.ndarray,
    x: np.ndarray,
    min_valid: int = 1,
) -> xr.Dataset:
    """
    Build monthly median NDVI composites from an irregular, cloud-masked series.

    ndvi_scenes : per-scene 2-D NDVI arrays, cloud pixels already set to NaN.
    dates       : acquisition dates, one per scene (ISO strings).
    min_valid   : mask a monthly pixel that had fewer than this many observations.
    """
    time = pd.to_datetime(dates)
    stack = xr.DataArray(
        np.stack(ndvi_scenes).astype("float32"),
        dims=("time", "y", "x"),
        coords={"time": time, "y": y, "x": x},
        name="ndvi",
    )

    grouped = stack.resample(time="1MS")           # calendar-month bins, month-start label
    composite = grouped.median()                    # skipna=True -> NaNs ignored
    valid_count = grouped.count()                   # finite observations per month

    # Suppress months with too few observations to be trustworthy
    composite = composite.where(valid_count >= min_valid)

    return xr.Dataset(
        {"ndvi_median": composite, "n_obs": valid_count}
    )


if __name__ == "__main__":
    rng = np.random.default_rng(0)
    ny, nx = 64, 64
    y = np.arange(ny, dtype="float32")
    x = np.arange(nx, dtype="float32")

    # Simulate 10 irregular acquisitions across three months, some cloudy (NaN)
    dates = [
        "2024-06-03", "2024-06-11", "2024-06-28",
        "2024-07-06", "2024-07-16", "2024-07-31",
        "2024-08-05", "2024-08-15", "2024-08-20", "2024-08-30",
    ]
    scenes = []
    for _ in dates:
        arr = rng.uniform(0.1, 0.9, size=(ny, nx)).astype("float32")
        cloud = rng.random((ny, nx)) < 0.3         # ~30% cloudy pixels -> NaN
        arr[cloud] = np.nan
        scenes.append(arr)

    ds = monthly_ndvi_composite(scenes, dates, y, x, min_valid=2)
    print(ds["ndvi_median"].time.values)           # 3 month-start timestamps
    print("Mean obs per month:", float(ds["n_obs"].mean()))

The valid_count layer is what turns a composite into something you can trust: a month-start pixel backed by one hazy observation is very different from one backed by six, and .where(valid_count >= min_valid) lets you drop the thin ones instead of plotting them as if they were solid. Carrying the count alongside the median also lets you weight or flag composites later — for a phenology fit, for instance, you might down-weight months with sparse coverage rather than discard them outright.

The whole aggregation is a single vectorised pass. resample never materialises a Python loop over dates; it computes the bin membership once and dispatches the median across the stack, so the cost scales with the array size rather than the number of scenes. That is also why the pattern transfers unchanged to a lazy, chunked stack: the same call builds a Dask task graph instead of computing eagerly, which is the escape hatch when a full time series will not fit in memory.

Variant Patterns

1. Seasonal or ten-day composites

The frequency string is the only change. Use QS-DEC for meteorological seasons or 10D for dekads; the reducer and NaN handling are identical.

seasonal = ndvi.resample(time="QS-DEC").median()   # DJF, MAM, JJA, SON
dekads = ndvi.resample(time="10D").median()        # 10-day composites

2. Maximum-value compositing

Classic MVC keeps the greenest observation per period, which strongly rejects residual cloud and off-nadir darkening. Swap the reducer for .max().

mvc = ndvi.resample(time="1MS").max()   # maximum-value composite, still skipna

3. Lazy compositing over a Dask-backed stack

If the NDVI stack is chunked with Dask, the same resample(...).median() call builds a task graph and runs out-of-core. Chunk the spatial dims and keep time whole so each period reduces within a partition.

lazy = ndvi.chunk({"time": -1, "y": 1024, "x": 1024})
monthly = lazy.resample(time="1MS").median().compute()

See Computing NDVI with xarray on a Dask Cluster for the Dask cluster setup this pattern targets.

Common Errors

The composite is dragged toward zero or a fill value

Cloud pixels are stored as 0 or a nodata sentinel rather than NaN, so the reducer treats them as real low-NDVI observations. Set masked pixels to np.nan before resampling so skipna excludes them.

ValueError: Could not convert time to a datetime index

The time coordinate is strings or plain integers. Convert it with pandas.to_datetime when building the DataArray so resample can bin it into calendar periods.

Empty months appear as 0 instead of NaN

A custom reducer or an explicit skipna=False filled all-NaN months with a computed value. Use the default .median() (which returns NaN for all-NaN groups) and pair it with a count() layer to flag empty periods.