⚙️ How-To Satellite imagery Raster & rasterio Satellite Imagery · step 14 of 16 Updated

How to Extract a Vegetation Index Time Series for a Polygon

Problem statement

Extracting NDVI for a field over a season looks like a loop: for each scene, mask the clouds, clip to the polygon, take the median. The output is a tidy table of dates and values that plots beautifully.

Almost all of the variation in that plot can be cloud.

Measured on a 1 km² vegetated polygon in Snowdonia across 120 Sentinel-2 scenes in 2025:

all 29 observations with any clear pixels    NDVI  -0.038 to 0.889,  sd 0.273
observations with the polygon >=90% clear    NDVI   0.759 to 0.877,  sd 0.043
observations with under 20% clear            NDVI  -0.012 to 0.805,  sd 0.248

The same polygon, the same year. Restricting to well-observed dates cuts the standard deviation by a factor of six. Everything else was residual cloud and shadow in a partly visible polygon.

Quick answer

Record the clear fraction with every observation, and filter on it:

import numpy as np


def polygon_observation(ndvi, usable, mask, min_fraction=0.5, min_pixels=30):
    """One time-series point, with the evidence behind it."""
    inside = mask & usable
    n = int(inside.sum())
    total = int(mask.sum())
    fraction = n / total if total else 0.0

    if n < min_pixels or fraction < min_fraction:
        return {"n": n, "fraction": round(fraction, 4), "median": None,
                "reason": "insufficient clear coverage"}

    values = ndvi[inside]
    return {
        "n": n,
        "fraction": round(fraction, 4),
        "median": round(float(np.median(values)), 4),
        "p25": round(float(np.percentile(values, 25)), 4),
        "p75": round(float(np.percentile(values, 75)), 4),
    }

A time-series point without its clear fraction is not a measurement. It is a number that might be the field or might be the cloud over it.

An NDVI time series scattering between minus 0.04 and 0.89 when unfiltered and collapsing to a narrow 0.76 to 0.88 band when restricted to observations that are 90 percent clear.
The spread is not phenology. It is the fraction of the polygon that was visible.

Step-by-step solution

1. Expect far fewer usable dates than scenes 

120 scenes over the area
 29 gave at least 2% of the polygon clear
  5 gave at least 90% clear

Five usable observations from 120 overpasses is what a 1 km² polygon in a maritime climate actually yields. Any analysis assuming a regular time step needs to confront that first.

2. Compute the clear fraction, not just the clear pixel count

The pixel count depends on polygon size; the fraction does not, and the fraction is what tells you whether the median is representative.

A polygon 10% visible gives a median of whatever corner happened to be clear. On this polygon, one 10.6%-clear observation reported 0.321 while a fully clear observation two weeks later reported 0.864.

3. Set a minimum clear fraction and state it

Fifty percent is a reasonable default; 90% is defensible when you have enough dates. Whatever you choose, put it in the output table so the filter travels with the numbers.

Look at what the threshold buys here:

filter observations NDVI range sd
any clear pixels 29 −0.038 to 0.889 0.273
at least 50% clear 11 0.486 to 0.889 0.130
at least 80% clear 6 0.759 to 0.877 0.039
at least 90% clear 5 0.759 to 0.877 0.043

Note that 20% is not a useful threshold here — it keeps 17 observations and barely moves the spread, because a 39%-clear date produced the series minimum of −0.038. The gain arrives between 50% and 80%.

4. Use the median, not the mean, within the polygon

Residual cloud edge is bright and pulls a mean down; a median over a polygon with a few contaminated pixels barely moves. Report the interquartile range too — a wide spread within a single date is itself evidence of contamination or genuine heterogeneity.

5. Watch for two observations on the same day

Sentinel-2 tiles overlap, so a polygon near a tile boundary can appear in two items with the same date. On this polygon, 2025-11-24 produced two observations: −0.038 at 39% clear and −0.012 at 8.2% clear.

Both are the same ground on the same day. Deduplicate by date, keeping the one with the higher clear fraction, or the analysis double-weights that day.

6. Interpret before you smooth

Winter NDVI here stays high — 0.862 on 27 December at 90% clear — because these are upland grasslands and conifer plantations that do not senesce. A gap-filling routine that assumes a temperate deciduous curve would "correct" that to something lower and wrong.

120 scenes narrowing to 29 with any clear pixels over the polygon and 5 with at least 90 percent clear.
The funnel from overpasses to usable observations is the real constraint on any time series.

Code examples

Example 1 — a time series with its evidence attached 

import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.features import geometry_mask

UNUSABLE = (0, 1, 3, 8, 9, 10)


def ndvi_series(items, geometry, aoi, min_fraction=0.5, min_pixels=30):
    """One row per scene: date, clear fraction, NDVI statistics."""
    rows = []
    for item in items:
        assets = {k: v.href for k, v in item.assets.items()}
        try:
            with rasterio.open(assets["scl"]) as ds:
                window = window_for(ds, aoi)
                scl = ds.read(1, window=window)
                transform = ds.window_transform(window)
            with rasterio.open(assets["red"]) as ds:
                red = ds.read(1, window=window_for(ds, aoi),
                              out_shape=scl.shape,
                              resampling=Resampling.bilinear).astype("float32")
            with rasterio.open(assets["nir"]) as ds:
                nir = ds.read(1, window=window_for(ds, aoi),
                              out_shape=scl.shape,
                              resampling=Resampling.bilinear).astype("float32")
        except rasterio.errors.RasterioIOError:
            continue

        inside = ~geometry_mask([geometry], out_shape=scl.shape,
                                transform=transform, invert=False)
        usable = ~np.isin(scl, UNUSABLE) & inside

        total = int(inside.sum())
        n = int(usable.sum())
        fraction = n / total if total else 0.0

        row = {"date": item.properties["datetime"][:10], "id": item.id,
               "clear_pixels": n, "clear_fraction": round(fraction, 4)}

        if n >= min_pixels and fraction >= min_fraction:
            with np.errstate(all="ignore"):
                total_band = nir + red
                ndvi = np.where(np.abs(total_band) < 1e-6, np.nan,
                                (nir - red) / total_band)
            values = ndvi[usable]
            values = values[np.isfinite(values)]
            row.update(ndvi_median=round(float(np.median(values)), 4),
                       ndvi_p25=round(float(np.percentile(values, 25)), 4),
                       ndvi_p75=round(float(np.percentile(values, 75)), 4))
        else:
            row["ndvi_median"] = None

        rows.append(row)

    usable_rows = [r for r in rows if r["ndvi_median"] is not None]
    print(f"  {len(rows)} scenes, {len(usable_rows)} usable at "
          f"{min_fraction:.0%} clear")
    return rows

Returning the rejected rows as well as the accepted ones is deliberate. A gap in a plot should be visibly a gap, not a missing entry nobody noticed.

Example 2 — deduplicating same-day observations 

from collections import defaultdict


def deduplicate_by_date(rows):
    """Overlapping tiles produce two items for one date. Keep the better."""
    best = {}
    for row in rows:
        current = best.get(row["date"])
        if current is None or row["clear_fraction"] > current["clear_fraction"]:
            best[row["date"]] = row

    duplicates = len(rows) - len(best)
    if duplicates:
        print(f"  removed {duplicates} same-day duplicate observation(s)")
    return sorted(best.values(), key=lambda r: r["date"])

Example 3 — how much your threshold changes the answer 

import numpy as np


def threshold_sensitivity(rows, thresholds=(0.0, 0.2, 0.5, 0.8, 0.9)):
    """Report the series statistics at several clear-fraction thresholds."""
    for threshold in thresholds:
        kept = [r for r in rows
                if r.get("ndvi_median") is not None
                and r["clear_fraction"] >= threshold]
        if not kept:
            print(f"  >= {threshold:4.0%} clear:  no observations")
            continue
        values = np.array([r["ndvi_median"] for r in kept])
        print(f"  >= {threshold:4.0%} clear:  {len(kept):3d} obs   "
              f"{values.min():6.3f} to {values.max():6.3f}   "
              f"sd {values.std():.3f}")
  >=   0% clear:   29 obs   -0.038 to  0.889   sd 0.273
  >=  20% clear:   17 obs   -0.038 to  0.889   sd 0.232
  >=  50% clear:   11 obs    0.486 to  0.889   sd 0.130
  >=  80% clear:    6 obs    0.759 to  0.877   sd 0.039
  >=  90% clear:    5 obs    0.759 to  0.877   sd 0.043

Run this before choosing a threshold. If the series statistics change a lot across it — as they do here, by a factor of six in the standard deviation — then the threshold is the most important parameter in your analysis and it belongs in the caption of every plot.

Explanation

Why partial visibility biases rather than adds noise

If cloud removed a random subset of pixels, a median over the rest would be an unbiased estimate. Cloud does not do that.

Cloud shadow is the largest unusable class, and shadowed pixels that escape the mask are dark in both bands, giving a plausible mid-range NDVI. Cloud edges are bright and give low NDVI. Neither is centred on the true value, so a partly visible polygon is biased in a direction that depends on where the cloud happened to be.

That is why the sub-20%-clear observations spanned −0.012 to 0.805 while the fully clear ones spanned 0.759 to 0.877. The former are not noisy measurements of the field; many of them are measurements of something else.

Why a fixed time step is the wrong model

Five usable observations across a year, unevenly spaced and clustered in the drier months, is not a time series in the sense most statistical tools assume. Applying a moving average, a seasonal decomposition or a Fourier fit to it imports an assumption of regular sampling that the data does not meet.

Better tools for this shape of data: interpolate to a regular grid explicitly and carry the interpolation flag, or fit a smooth function of day-of-year (a harmonic regression or a spline) directly to the irregular observations, weighting each by its clear fraction.

Why per-pixel and per-polygon series behave differently

A per-pixel series over the same period had a median of 14 clear observations across the study window, because different parts of the window are clear on different dates. A whole 1 km polygon needs all of it clear at once, which is far rarer — five dates against fourteen.

So the choice between "NDVI per pixel, then aggregate" and "aggregate first, then NDVI" is not cosmetic. Per-pixel series have many more observations and no requirement that the polygon be simultaneously visible; the cost is that a polygon summary built from them mixes dates.

Why local knowledge beats gap-filling

The December value of 0.862 is real. Upland Wales is grass and conifer, and neither drops its leaves.

Any gap-filling or smoothing routine encodes an expectation about the seasonal shape, and that expectation is a model of a particular vegetation type. Applied to the wrong one it produces smooth, plausible, wrong curves — and because they are smooth and plausible, they survive review.

Three time-series symptoms and their causes: duplicate dates, partial visibility and gap-filling.
All three produce a plausible curve. Only the clear-fraction column reveals the second.

Edge cases or notes

  • Record the clear fraction with every observation. Without it a value is uninterpretable.
  • Deduplicate by date. Overlapping tiles produce two items for one day.
  • Use the median within the polygon, and report the interquartile range.
  • A polygon smaller than a few pixels cannot support a median at all; at 10 m, a 0.5 ha field is 50 pixels.
  • Mixed 10 m and 20 m bands must be resampled onto one grid before the index.
  • geometry_mask needs the window transform, not the dataset transform.
  • Sentinel-2A, 2B and 2C differ slightly in spectral response; a step between platforms is instrumental, not seasonal.
  • Do not smooth before you have looked at the raw points. The smoothing will hide the cloud problem rather than solve it.

FAQ

How do I get an NDVI time series for a field in Python?

For each scene, mask clouds, clip to the polygon, and take the median NDVI over the clear pixels — recording the fraction of the polygon that was clear with every value.

Why is my NDVI time series so noisy?

Usually because partly cloudy dates are included. Filtering to observations with the polygon at least 90% clear cut the standard deviation from 0.273 to 0.043 on the series measured here.

What clear fraction should I require?

At least 50%, and 90% if you have enough dates. Run the sensitivity table first — if the series statistics move a lot across thresholds, the threshold is your most important parameter.

How many usable observations should I expect?

Far fewer than the number of scenes. Over a 1 km² polygon in Snowdonia, 120 scenes in 2025 gave 29 with any clear pixels and 5 that were at least 90% clear.

Why do I have two observations on the same date?

Sentinel-2 tiles overlap, so a polygon near a boundary appears in two items. Keep the one with the higher clear fraction.

Should I smooth or gap-fill the series?

Only after looking at the raw points and knowing the vegetation. Smoothing encodes an assumed seasonal shape, and applying a deciduous curve to evergreen upland grassland produces plausible, wrong output.

Mean or median within the polygon?

Median. Residual cloud edge is bright and drags a mean; a median tolerates a few contaminated pixels.