How to extract phenology metrics from an NDVI time series

Problem statement

The metrics are simple to define and sensitive to every choice around them. Start of season, peak, end of season, length and integral are all read off a curve — and the curve is an interpolation through irregular, cloud-limited observations, smoothed by a window you chose, thresholded at a level you chose.

On a real Sentinel-2 series over a Flevoland spring-barley field in 2025, 39 usable dates ran from 7 March to 30 September, and the field-median NDVI went 0.163 → 0.959 → 0.163 before the stubble greened again. The metrics from that series are defensible. The same metrics with a 50% threshold instead of 20% start 9 days later and end 10 days earlier, and neither number is wrong.

Quick answer

import numpy as np, pandas as pd

def phenology(dates, values, threshold=0.2, smooth_days=15):
    s = (pd.Series(np.asarray(values, dtype=float), index=pd.DatetimeIndex(dates))
         .reindex(pd.date_range(min(dates), max(dates), freq="D"))
         .interpolate("time")
         .rolling(smooth_days, center=True, min_periods=smooth_days // 3).mean()
         .dropna())

    base, peak = float(s.min()), float(s.max())
    amp = peak - base
    peak_date = s.idxmax()
    level = base + threshold * amp
    rising = s[(s >= level) & (s.index <= peak_date)]
    falling = s[(s >= level) & (s.index >= peak_date)]

    return {"base": round(base, 3), "peak": round(peak, 3),
            "amplitude": round(amp, 3), "peak_date": str(peak_date.date()),
            "sos": str(rising.index.min().date()) if len(rising) else None,
            "eos": str(falling.index.max().date()) if len(falling) else None,
            "integral": round(float((s - base).clip(lower=0).sum()), 1),
            "threshold": threshold, "smooth_days": smooth_days,
            "observations": int(np.isfinite(values).sum())}

Returning the threshold, the smoothing window and the observation count inside the result is the point: without them the metrics cannot be compared with anyone else's, including your own from last year.

Gerst, zomer- parcel, 11.4 ha; 47 dates with a clear view of it, 2025-03-07 to 2025-09-30

Step-by-step solution

1. Build the series per field

Field medians are far less noisy than pixel series and the field is the management unit. Use pixel series only where within-field variation is the question.

2. Screen the observations before interpolating

Cloud, shadow and snow all produce valid-looking index values. A scene classification mask and a minimum clear fraction per date are the cheapest screen; on the reference series a 60% threshold kept 39 of the 47 dates that had any clear view of the field.

3. Interpolate onto a daily grid

Linear interpolation in time is adequate for most work. The important part is that every metric is computed from a regular grid, so a gap does not silently become a feature.

4. Smooth, and report the window

A 15-day centred rolling mean removes observation noise. It also shifts the peak of an asymmetric curve by roughly half the window, which is a systematic bias that belongs in the documentation.

5. Choose a threshold definition

Relative to amplitude (20% is common), absolute (NDVI 0.3), or the maximum of the derivative. Relative thresholds transfer between crops; absolute ones are easier to explain; the derivative maximum is the most physically meaningful and the noisiest.

6. Handle the double peak

Cover crops, cut grassland and catch crops give two peaks. A single-peak extractor reports a season running from the first green-up to the last senescence, which is not a season.

7. Report the evidence behind each metric

How many clear observations fall within ten days of the start of season, the peak and the end. A metric interpolated across a 24-day gap deserves less confidence than one bracketed by observations a day or two apart.

Code examples

The data used on this page

One real field: the largest spring-barley parcel (Gerst, zomer-) declared in 2025 in the Dutch parcel register BRP Gewaspercelen within 5.55–5.75° E, 52.48–52.56° N, east of Dronten in Flevoland. Its series comes from every Sentinel-2 L2A scene in Earth Search over the parcel from 1 March to 30 September 2025. For each date the setup records the share of the parcel that the scene classification calls vegetation or bare soil (clear_frac), and the median NDVI (B8, B4) and NDRE (B8, B5) over those pixels. Where two scenes share a date, the clearer one is kept. Dates with no clear pixel are dropped. Reading about a hundred scenes takes several minutes the first time; after that season.csv is reused.

import os

import geopandas as gpd
import numpy as np
import pandas as pd
import rasterio
import requests
from rasterio.features import geometry_mask
from rasterio.windows import from_bounds

os.environ.update(AWS_NO_SIGN_REQUEST="YES", GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR",
                  CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif")

if not os.path.exists("field.gpkg"):
    url = "https://api.pdok.nl/rvo/gewaspercelen/ogc/v1/collections/brpgewas/items"
    feats, params = [], {"f": "json", "bbox": "5.55,52.48,5.75,52.56", "limit": 1000}
    while url:
        page = requests.get(url, params=params, timeout=120).json()
        feats += page["features"]
        url, params = next((l["href"] for l in page["links"] if l["rel"] == "next"), None), None
    brp = gpd.GeoDataFrame.from_features(feats, crs=4326)
    brp["ha"] = brp.to_crs(28992).area / 1e4
    brp.to_file("brp.gpkg")
    brp[brp.gewas.str.startswith("Gerst, zomer")].nlargest(1, "ha").to_file("field.gpkg")
field = gpd.read_file("field.gpkg")


def window(href, bounds, shape=None):
    with rasterio.open(href) as src:
        win = from_bounds(*bounds, transform=src.transform)
        extra = {"out_shape": shape} if shape else {}
        return src.read(1, window=win, boundless=True, **extra).astype("float32"), src.window_transform(win)


if not os.path.exists("season.csv"):
    body = {"collections": ["sentinel-2-l2a"], "intersects": field.geometry.iloc[0].__geo_interface__,
            "datetime": "2025-03-01T00:00:00Z/2025-09-30T23:59:59Z", "limit": 100}
    items = []
    while body:
        page = requests.post("https://earth-search.aws.element84.com/v1/search", json=body, timeout=120).json()
        items += page["features"]
        body = next((l["body"] for l in page["links"] if l["rel"] == "next"), None)
    rows = []
    for day in sorted({i["properties"]["datetime"][:10] for i in items}):
        best = None
        for item in (i for i in items if i["properties"]["datetime"][:10] == day):
            a = item["assets"]
            utm = field.to_crs(item["properties"]["proj:epsg"])
            b = utm.total_bounds + np.array([-20, -20, 20, 20])
            red, transform = window(a["red"]["href"], b)
            inside = ~geometry_mask(utm.geometry, red.shape, transform)
            scl, _ = window(a["scl"]["href"], b, red.shape)
            clear = inside & np.isin(scl, [4, 5])          # vegetation, bare soil
            if best is None or clear.sum() > best[1].sum():
                best = (item, clear, inside, b, red)
        item, clear, inside, b, red = best
        if not clear.any():
            continue
        # DN / 10,000 is reflectance: no offset to remove, whatever the item's
        # boa_offset_applied flag says (water reads 1.5–2.6% on 7 and 10 March either way)
        nir, _ = window(item["assets"]["nir"]["href"], b, red.shape)
        re1, _ = window(item["assets"]["rededge1"]["href"], b, red.shape)
        r, n, e = red[clear], nir[clear], re1[clear]
        rows.append({"date": day, "scene": item["id"],
                     "clear_frac": round(float(clear.sum() / inside.sum()), 3),
                     "ndvi_med": round(float(np.median((n - r) / (n + r))), 3),
                     "ndre_med": round(float(np.median((n - e) / (n + e))), 3)})
    pd.DataFrame(rows).to_csv("season.csv", index=False)

season = pd.read_csv("season.csv", parse_dates=["date"])
print(f"{field.gewas.iloc[0]} parcel, {field.ha.iloc[0]:.1f} ha; "
      f"{len(season)} dates with a clear view of it, {season.date.min().date()} to {season.date.max().date()}")

Gerst, zomer- parcel, 11.4 ha; 47 dates with a clear view of it, 2025-03-07 to 2025-09-30

Example 1 — the real series, and the metrics from it

usable = season[season.clear_frac >= 0.6]
crop = usable[usable.date <= "2025-09-15"]        # before the stubble regrows
print(f"{len(usable)} usable dates of {len(season)}; the barley season uses {len(crop)}")
print(crop.iloc[::4][["date", "ndvi_med", "ndre_med"]].to_string(index=False))

print(phenology(crop.date, crop.ndvi_med))
print(phenology(crop.date, crop.ndre_med))
for t in (0.1, 0.2, 0.3, 0.5):
    m = phenology(crop.date, crop.ndvi_med, threshold=t)
    days = (pd.Timestamp(m["eos"]) - pd.Timestamp(m["sos"])).days
    print(f"threshold {t:.0%}: {m['sos']} to {m['eos']}, {days} days")
plateau = usable[usable.ndvi_med > 0.9]
print(f"NDVI above 0.9 from {plateau.date.min():%d %b} to {plateau.date.max():%d %b}: "
      f"{plateau.ndvi_med.min():.3f}–{plateau.ndvi_med.max():.3f}")

39 usable dates of 47; the barley season uses 37
      date  ndvi_med  ndre_med
2025-03-07     0.163     0.103
2025-04-01     0.400     0.233
2025-04-09     0.558     0.383
2025-04-29     0.938     0.733
2025-05-13     0.938     0.770
2025-06-12     0.899     0.698
2025-06-30     0.510     0.268
2025-07-20     0.163     0.091
2025-08-14     0.200     0.103
2025-09-13     0.281     0.153
{'base': 0.169, 'peak': 0.95, 'amplitude': 0.781, 'peak_date': '2025-05-18', 'sos': '2025-03-28', 'eos': '2025-07-07', 'integral': 67.1, 'threshold': 0.2, 'smooth_days': 15, 'observations': 37}
{'base': 0.105, 'peak': 0.792, 'amplitude': 0.687, 'peak_date': '2025-05-22', 'sos': '2025-03-30', 'eos': '2025-07-02', 'integral': 51.3, 'threshold': 0.2, 'smooth_days': 15, 'observations': 37}
threshold 10%: 2025-03-23 to 2025-07-14, 113 days
threshold 20%: 2025-03-28 to 2025-07-07, 101 days
threshold 30%: 2025-03-30 to 2025-07-03, 95 days
threshold 50%: 2025-04-06 to 2025-06-27, 82 days
NDVI above 0.9 from 26 Apr to 13 Jun: 0.923–0.959

Running the same function on NDVI and on NDRE is worth doing. NDVI sits above 0.92 from late April to mid-June, so its peak date is poorly determined; the NDRE curve still has a peak in it, on 22 May against NDVI's 18 May. The threshold alone moves the season length from 82 to 113 days.

Example 2 — the evidence behind each metric

import numpy as np, pandas as pd

def metric_support(dates, metric_date, window_days=10):
    d = pd.DatetimeIndex(dates)
    m = pd.Timestamp(metric_date)
    within = np.abs((d - m).days) <= window_days
    gaps = np.diff(np.sort(d)).astype("timedelta64[D]").astype(int)
    before = d[d <= m].max() if (d <= m).any() else None
    after = d[d >= m].min() if (d >= m).any() else None
    return {"observations_within": int(within.sum()),
            "days_to_previous": (m - before).days if before is not None else None,
            "days_to_next": (after - m).days if after is not None else None,
            "largest_gap_days": int(gaps.max()) if len(gaps) else None}

m = phenology(crop.date, crop.ndvi_med)
for key in ("sos", "peak_date", "eos"):
    print(key, m[key], metric_support(crop.date, m[key]))

sos 2025-03-28 {'observations_within': 5, 'days_to_previous': 8, 'days_to_next': 4, 'largest_gap_days': 24}
peak_date 2025-05-18 {'observations_within': 6, 'days_to_previous': 2, 'days_to_next': 1, 'largest_gap_days': 24}
eos 2025-07-07 {'observations_within': 3, 'days_to_previous': 5, 'days_to_next': 8, 'largest_gap_days': 24}

All three metrics here are well bracketed — the nearest observations are one to eight days away. The 24-day gap from 19 May to 12 June is in the record, though, and in a cloudier year the same check is what tells you a metric was interpolated rather than seen.

Example 3 — a double-peak-aware extractor

import numpy as np, pandas as pd
from scipy.signal import find_peaks

def seasons(dates, values, threshold=0.2, smooth_days=15, min_prominence=0.1,
            min_separation_days=45):
    s = (pd.Series(np.asarray(values, dtype=float), index=pd.DatetimeIndex(dates))
         .reindex(pd.date_range(min(dates), max(dates), freq="D"))
         .interpolate("time")
         .rolling(smooth_days, center=True, min_periods=smooth_days // 3).mean()
         .dropna())

    base, amp = float(s.min()), float(s.max() - s.min())
    idx, props = find_peaks(s.values, prominence=min_prominence * amp,
                            distance=min_separation_days)
    if len(idx) == 0:
        return []

    out = []
    level = base + threshold * amp
    for i, p in enumerate(idx):
        left = s.index[:p][s.values[:p] < level]
        right = s.index[p:][s.values[p:] < level]
        out.append({
            "peak_date": str(s.index[p].date()),
            "peak": round(float(s.values[p]), 3),
            "prominence": round(float(props["prominences"][i]), 3),
            "sos": str(left.max().date()) if len(left) else None,
            "eos": str(right.min().date()) if len(right) else None,
        })
    return out


for label, part in (("March to mid-September", crop), ("all of March to September", usable)):
    print(label)
    for found in seasons(part.date, part.ndvi_med):
        print("  ", found)
    print("   single-peak extractor:",
          {k: phenology(part.date, part.ndvi_med)[k] for k in ("sos", "eos")})

March to mid-September
   {'peak_date': '2025-05-18', 'peak': 0.95, 'prominence': 0.778, 'sos': '2025-03-27', 'eos': '2025-07-08'}
   single-peak extractor: {'sos': '2025-03-28', 'eos': '2025-07-07'}
all of March to September
   {'peak_date': '2025-05-18', 'peak': 0.95, 'prominence': 0.778, 'sos': '2025-03-27', 'eos': '2025-07-08'}
   single-peak extractor: {'sos': '2025-03-28', 'eos': '2025-09-30'}

prominence is the parameter that separates a genuine second crop from a wobble in the senescence curve. Setting it as a fraction of the seasonal amplitude makes it transfer between fields and indices.

The September regrowth shows the other half. Given all of March to September, the single-peak extractor ran the season on to 30 September, because the stubble had climbed back over its threshold. seasons() still ended it on 8 July: it looks for the first fall below the threshold after the peak. The regrowth is not reported as a second season either, because the series ends while it is still rising and find_peaks does not count a rise at the edge.

Explanation

Why interpolation is not optional

Metrics computed from the observations alone are metrics about when the sky was clear. Across a 24-day gap over the peak, the "highest observed value" is whichever cloud-free date happened to fall either side of it. Interpolating onto a daily grid makes the metric a property of the curve, which is at least a consistent thing to compare between fields.

Why smoothing shifts the peak

A centred rolling mean is symmetric, so on a symmetric curve it leaves the peak in place. A crop curve is not symmetric — green-up is faster than senescence — so the mean pulls the peak towards the slower side by roughly half the window for a strongly asymmetric curve. The bias is systematic and therefore comparable between fields, provided the window is the same.

Why relative thresholds transfer and absolute ones do not

An absolute NDVI threshold of 0.3 is early emergence for maize and mid-season for a thin grassland. A threshold at 20% of the field's own amplitude adapts to the crop and the soil background, which is what makes it usable across a mixed farm. The cost is that a field with a small amplitude gets a threshold close to its noise floor.

Why the peak date is unreliable with NDVI

From 26 April to 13 June the reference series sat between 0.923 and 0.959. Any peak-finding on that plateau is picking between values separated by little more than the measurement noise, so the peak date is effectively arbitrary within a seven-week window. A red-edge index has a real maximum in the same period.

Edge cases or notes

• Winter crops span the year end. Define the season window explicitly.
• Perennials never reach a low base. Amplitude-relative thresholds fail on orchards.
• Snow gives negative values. Mask rather than interpolating across.
• Series must be sorted. Duplicated dates break the interpolation.
• A partial series has no peak. Require coverage before and after.
• Savitzky–Golay and double-logistic fits impose a shape; the shape is an assumption.
• The integral is the most robust metric and the least interpretable.
• Report the observation count alongside every metric.

Internal links

• Crop phenology and growing seasons explained — what the metrics mean
• Vegetation indices explained: NDVI, EVI, NDRE and when each fails — why the NDVI peak is unreliable
• How to classify crop types from a satellite time series — the main consumer
• Growing degree days explained — expressing the metrics in thermal time
• How to build an NDVI time series for a polygon in Python — assembling the series
• How to resample a time series in xarray — the gridded equivalent
• How to mask clouds in Sentinel-2 imagery in Python — screening the dates
• NDVI stops responding in a dense canopy — the flat peak

FAQ

What phenology metrics should I extract?

Start of season, peak value and date, end of season, length, amplitude and the integral. The integral is the most robust and the least interpretable.

Do I have to interpolate the series?

Yes. Metrics computed from raw observations are metrics about when the sky was clear — even a fairly clear real season had a 24-day gap across the peak.

What threshold defines the start of season?

Whatever you publish. Twenty per cent of the seasonal amplitude is common and transfers between crops better than an absolute value.

Does smoothing bias the metrics?

It shifts the peak of an asymmetric curve by roughly half the window. The bias is systematic, so it is comparable as long as the window is fixed.

Why is my peak date unstable?

NDVI saturation. A real barley field sat between 0.923 and 0.959 for seven weeks, so the peak date within that window is arbitrary.

How do I handle two crops in one season?

Use a peak finder with a prominence threshold expressed as a fraction of the seasonal amplitude, and return one set of metrics per peak.