Building Temporal Feature Stacks from Image Time Series

Collapse the masked series into percentiles per season, and keep the count of what went into them:

import xarray as xr

masked = ndvi.where(cloud_free)                      # never reduce unmasked data
season = masked.groupby("time.season")

features = xr.Dataset({
    "ndvi_p10": season.quantile(0.10, dim="time"),
    "ndvi_p50": season.quantile(0.50, dim="time"),
    "ndvi_p90": season.quantile(0.90, dim="time"),
    "n_obs": masked.groupby("time.season").count(dim="time"),
})

Percentiles survive residual cloud; means do not. This page belongs to feature engineering for pixel-based models in Raster Machine Learning & Model Inference.

Phenology Is the Strongest Signal in the Data

A single date tells you what a surface looked like on one morning. A year tells you what it does — when it greens up, how long it stays green, whether it is cut twice or grazed continuously. For anything involving vegetation, the temporal shape separates classes that are spectrally identical on any given date.

Four classes, one crossing point Winter wheat peaks early and is harvested by midsummer, maize peaks late, grassland stays moderately green all year and forest is high and flat. All four curves pass through nearly the same value in late spring, so a single image from that date cannot separate them, while the full-year shape separates them cleanly. Annual NDVI trajectories 0.9 0.1 a single date here separates nothing Jan Jun Dec winter wheat maize grassland forest

The engineering problem is that the series is irregular and full of holes. Different pixels have different valid dates, so any feature that assumes an even sampling has to reckon with that first.

Environment & Setup

Package Version pin Used for
xarray >=2023.12 Time-indexed reductions and grouping
rioxarray >=0.15 Loading the series and writing the feature stack
dask >=2024.1 Chunked reduction when the cube exceeds memory
numpy >=1.23 Harmonic design matrix and least squares 
pip install "xarray>=2023.12" "rioxarray>=0.15" "dask>=2024.1" "numpy>=1.23"

Complete Working Example

import numpy as np
import xarray as xr


def temporal_features(series: xr.DataArray, *, min_obs: int = 6) -> xr.Dataset:
    """Robust seasonal features from an irregular, masked NDVI series.

    series: dims (time, y, x), already masked so invalid observations are NaN.
    """
    out = {}
    for name, group in series.groupby("time.season"):
        q = group.quantile([0.1, 0.5, 0.9], dim="time")
        n = group.count(dim="time")
        enough = n >= min_obs

        out[f"{name}_p10"] = q.sel(quantile=0.1).where(enough).astype("float32")
        out[f"{name}_p50"] = q.sel(quantile=0.5).where(enough).astype("float32")
        out[f"{name}_p90"] = q.sel(quantile=0.9).where(enough).astype("float32")
        out[f"{name}_range"] = (out[f"{name}_p90"] - out[f"{name}_p10"]).astype("float32")
        out[f"{name}_n"] = n.astype("float32")

    ds = xr.Dataset(out)
    # Amplitude and timing across the whole year, independent of season bins
    ds["annual_amplitude"] = (series.max("time") - series.min("time")).astype("float32")
    ds["peak_doy_sin"], ds["peak_doy_cos"] = _peak_timing(series)
    return ds.drop_vars("quantile", errors="ignore")


def _peak_timing(series: xr.DataArray) -> tuple[xr.DataArray, xr.DataArray]:
    """Day-of-year of the maximum, encoded circularly."""
    doy = xr.DataArray(series["time"].dt.dayofyear, dims="time")
    peak = doy.isel(time=series.fillna(-9).argmax("time"))
    angle = 2 * np.pi * peak / 365.25
    return np.sin(angle).astype("float32"), np.cos(angle).astype("float32")

Peak timing is encoded circularly for the same reason aspect is: day 360 and day 5 are five days apart on the ground and 355 apart as integers. The fillna(-9) before argmax is a guard — argmax on an all-NaN pixel raises, and those pixels are masked out downstream by the min_obs condition anyway.

Choosing Between Reduction Strategies

Four ways to collapse a time series Annual mean is one feature, cheap, and destroyed by residual cloud. Seasonal percentiles are a dozen features, robust, and tolerant of gaps. A harmonic fit is six features and describes shape well but needs an evenly sampled series. Raw stacked dates preserve everything and explode the feature count while being unusable where dates differ per pixel. Reduction strategy trade-offs annual mean 1 feature cloud-sensitive gap tolerant too coarse for crops seasonal percentiles 12-20 features robust to cloud gap tolerant the default choice harmonic fit 4-6 features moderately robust needs even sampling elegant where data allows stacked dates 30-70 features cloud-sensitive breaks on gaps only for dense archives Percentiles win because they are the only column that is robust and gap tolerant at once. Add a harmonic fit on top only where the archive is dense enough to support it.

The harmonic option deserves a note because it is genuinely better where it applies. Fitting a + b·sin(2πt) + c·cos(2πt) per pixel gives an amplitude and a phase that describe the season in four numbers, and the residual describes how unusual the year was. It needs a reasonably even series, which in practice means a dense archive or an interpolation step first — the techniques in filling gaps in NDVI time series with interpolation apply directly.

Interpolating for percentile features, on the other hand, is actively harmful: it replaces missing observations with a smooth guess, which narrows the distribution and biases the tails inward. The percentiles then describe the interpolation rather than the surface.

Verification

import numpy as np

for name in ds.data_vars:
    if name.endswith("_n"):
        continue
    finite = float(np.isfinite(ds[name]).mean())
    print(f"{name:<18} {finite:6.1%} finite")

assert float(ds["DJF_n"].max()) <= len(series.time), "count exceeds available dates"
assert float((ds["JJA_p90"] >= ds["JJA_p10"]).mean(skipna=True)) > 0.999
Valid observations are not evenly distributed Across one region the dry-season feature is built from fifteen to twenty valid dates almost everywhere, while the wet-season feature falls to three or four dates over the uplands where cloud persists. Carrying the count as its own feature lets the model discount the pixels where the summary rests on almost nothing. dry season: 15-20 dates wet season: 3-12 dates uniform coverage; percentiles well determined pale uplands have too few dates to trust Darker means more valid observations. The count belongs in the feature stack, not just the log.

The ordering assertion is the useful one: a p90 below a p10 means the quantile dimension was selected wrongly, which is easy to do and produces features that look plausible. The finite fractions expose the seasons where cloud left too few observations — in monsoon regions the wet-season features are routinely 40% missing, and the model must be given the count feature so it can discount them rather than treating a noisy percentile as gospel. Where the count falls below the threshold across a whole region rather than scattered pixels, the honest response is to drop that season’s features from the model entirely for that region, because a feature that is absent in half the study area teaches the model to use geography as a proxy for data availability.

Common Errors

ValueError: All-NaN slice encountered

A pixel had no valid observation in a season and a reduction that cannot handle it was used. Guard with the min_obs mask before computing, or use the nan-aware variants throughout.

Features vary with the number of dates in the request

The series was reduced without masking, so cloudy scenes shifted the statistics. Always mask first, using the approach in masking clouds with the Sentinel-2 SCL band.

Memory blows up on a full tile

Percentiles require sorting along time, which Dask cannot do lazily in one pass. Rechunk so each chunk holds the whole time axis for a small spatial block — the guidance in tuning Dask chunk sizes for raster cubes applies directly.

Frequently Asked Questions

Q: Why percentiles rather than means? Because residual cloud survives every mask, and a single bright contaminated observation moves a mean far more than it moves a median or a tenth percentile. Percentiles also describe the shape of the season, which is what separates crop types.

Q: Should I interpolate gaps before computing features? Not for percentile features — they handle gaps natively by ignoring missing dates. Interpolate only when the feature genuinely requires an evenly spaced series, such as a harmonic fit or a rolling window.

Q: How many observations does a pixel need? Enough that the percentile is not the extreme value itself: at least six to eight valid dates per season for a tenth percentile to mean anything. Below that, mask the pixel out rather than letting it contribute a noisy feature.

Q: Should the seasons be calendar quarters or locally defined? Locally defined, wherever the growing season is known. Calendar quarters split a southern-hemisphere growing season across two bins and merge a monsoon onset with the dry period before it. Defining seasons from the local agricultural calendar — or simply from the long-term mean NDVI curve for the region — produces features that line up with what the vegetation is actually doing, and it costs nothing more than a different grouping key.