Computing Phenology Metrics from NDVI Time Series

Smooth a regular series, then find where it crosses a fraction of its seasonal amplitude:

import numpy as np
from scipy.signal import savgol_filter

s = savgol_filter(ndvi_5day, window_length=7, polyorder=2)   # regular 5-day series
base, peak = np.percentile(s, 5), s.max()
thr = base + 0.5 * (peak - base)
above = np.where(s >= thr)[0]
sos, pos, eos = doy[above[0]], doy[s.argmax()], doy[above[-1]]

Phenology metrics turn a noisy year of observations into a handful of numbers — when growth started, when it peaked, when it ended, and how strong it was. This page belongs to temporal aggregation and time series analysis in Satellite Processing Workflows & Index Pipelines.

The Metrics on One Curve

A relative threshold — a fraction of each pixel’s own amplitude — is the standard approach because it adapts to land cover. A fixed NDVI threshold of 0.4 would mark the start of season for a lush crop in early spring and never trigger for a sparse rangeland, while half-amplitude marks an equivalent point on both curves.

Environment & Setup

pip install "xarray>=2023.1" "scipy>=1.11" "numpy>=1.23" "dask>=2023.1"

Complete Working Example

import numpy as np
import xarray as xr
from scipy.signal import savgol_filter


def regularise(ndvi: xr.DataArray, step: str = "5D", max_gap: str = "30D") -> xr.DataArray:
    reg = ndvi.resample(time=step).max()                     # max favours clear over hazy
    return reg.interpolate_na("time", max_gap=max_gap)


def _metrics_1d(y: np.ndarray, doy: np.ndarray, frac: float) -> np.ndarray:
    out = np.full(6, np.nan, dtype="float32")
    if np.isnan(y).mean() > 0.3:
        return out
    y = np.interp(np.arange(len(y)), np.flatnonzero(~np.isnan(y)), y[~np.isnan(y)])
    s = savgol_filter(y, window_length=7, polyorder=2)
    base, ipk = np.percentile(s, 5), int(np.argmax(s))
    amp = s[ipk] - base
    if amp < 0.1:                                            # no real season
        return out
    thr = base + frac * amp
    up = np.flatnonzero(s[:ipk + 1] < thr)
    down = np.flatnonzero(s[ipk:] < thr)
    if len(up) == 0 or len(down) == 0:
        return out
    i_sos, i_eos = up[-1] + 1, ipk + down[0] - 1
    integ = float(np.trapz(np.clip(s[i_sos:i_eos + 1] - base, 0, None), doy[i_sos:i_eos + 1]))
    out[:] = [doy[i_sos], doy[ipk], doy[i_eos], doy[i_eos] - doy[i_sos], amp, integ]
    return out


def phenology(ndvi: xr.DataArray, frac: float = 0.5) -> xr.Dataset:
    reg = regularise(ndvi)
    doy = reg.time.dt.dayofyear.values.astype("float32")
    res = xr.apply_ufunc(
        _metrics_1d, reg, kwargs={"doy": doy, "frac": frac},
        input_core_dims=[["time"]], output_core_dims=[["metric"]],
        vectorize=True, dask="parallelized", output_dtypes=["float32"],
        dask_gufunc_kwargs={"output_sizes": {"metric": 6}})
    names = ["sos", "pos", "eos", "length", "amplitude", "integrated"]
    return res.assign_coords(metric=names).to_dataset("metric")

The search for the start of season walks backwards from the peak rather than forwards from January. Walking forwards finds the first crossing of the year, which is often a noise spike or a winter cover crop; walking back from the peak finds the crossing that belongs to the main season. The same logic, mirrored, finds the end. Gap filling before smoothing follows filling gaps in NDVI time series with interpolation.

Smoothing Choices and Their Bias

Residual cloud pulls NDVI down, never up, so noise in these series is one-sided. Taking the maximum within each regularisation step, as the example does, favours the clearest observation; a further refinement is an upper-envelope fit that iteratively replaces values below the smoothed curve with the curve itself, which is how established tools such as TIMESAT treat cloud-contaminated series. Whichever smoother is used, keep it identical across years, because a change of smoother changes every metric.

Two Seasons in One Year

Double-cropped fields, and savannas with two rainy seasons, have two peaks. The single-peak method reports only the larger and computes a season length that may span both. Detecting a second season means finding local maxima separated by a minimum that falls below the threshold, then applying the same crossing logic around each peak. Report a count of seasons as its own layer so users know which pixels carry a second set of metrics. For most temperate applications one season is enough, but in South and Southeast Asia assuming one season silently discards half the agricultural calendar.

Storing and Comparing Years

Phenology metrics become most useful as a multi-year record. Write each year’s metrics as a small multi-band COG — start, peak, end, length, amplitude, integrated NDVI and a season count — with day-of-year stored as int16 and a nodata value outside the valid range. Comparing a year against the average of previous years then reveals early or late springs, shortened seasons in drought years, and fields whose management changed. Because every metric depends on the smoother, threshold fraction and regularisation step, record those settings in the file’s tags; a two-day shift in start of season is meaningless if the method changed between years, so treat the settings as part of the data itself.

Verification

import numpy as np

ph = phenology(ndvi_2025).compute()
ok = ph.sos.notnull()
assert bool(((ph.sos < ph.pos) & (ph.pos < ph.eos)).where(ok, True).all())
print("median length", float(ph.length.median()), "days;",
      "pixels with a season", f"{float(ok.mean()):.1%}")

Common Errors

Start of season in January for every pixel

The search runs forward from the start of the year and finds noise. Search backwards from the peak.

Amplitude is lower than expected

A moving-average smoother flattened the peak. Use Savitzky–Golay or an upper-envelope fit.

Many pixels have no season

The amplitude floor is too high for sparse vegetation, or too many observations are missing. Lower the floor or lengthen the gap-fill limit.

Seasons that cross New Year are cut in half

Southern-hemisphere and winter crops span two calendar years. Run on a July-to-June window for those regions.

Frequently Asked Questions

Q: What threshold should define start of season? Half the seasonal amplitude is the most common choice and is robust. Lower fractions such as 20% find earlier green-up but are more sensitive to noise.

Q: Why smooth before extracting metrics? Raw series contain cloud dips and noise that create false crossings. Smoothing gives a curve with one clear rise and fall per season.

Q: How do I handle seasons that cross the new year? Shift the analysis window — for example July to June — so the whole season falls inside one window, and convert day-of-year accordingly.

Q: Is integrated NDVI a measure of productivity? It is a widely used proxy for seasonal productivity, since it combines how green and how long. It is relative, not an absolute biomass measurement, unless calibrated against field data.

Related

• Temporal Aggregation and Time Series Analysis — the parent topic.
• Filling Gaps in NDVI Time Series with Interpolation — preparing the series.
• Detecting Seasonal Trends with Rolling Windows — trends across years.
• Building Temporal Feature Stacks from Image Time Series — phenology metrics as model features.