How to Extract a Vegetation Index Time Series for a Polygon

On this page

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.

Much of the variation in that plot is not the vegetation. Some is cloud in a partly visible polygon; in an upland winter, some is snow and low sun, which a cloud mask was never going to remove.

Measured on a 1 km² square near Snowdon across the 127 Sentinel-2 scenes from June to December 2025 (Examples 3 and 4):

all 29 observations with any clear pixels    NDVI  -0.092 to 0.804,  sd 0.295
observations with the square >=90% clear     NDVI   0.170 to 0.762,  sd 0.200
observations with under 20% clear            NDVI  -0.067 to 0.804,  sd 0.246
also dropping snow, dark area, unclassified,
  and >=50% clear                            NDVI   0.513 to 0.778,  sd 0.092

Filtering on the cloud-mask clear fraction alone cut the standard deviation by a third. The 90%-clear set still held 19 November at NDVI 0.170, when 58.6% of the square was snow. Dropping the three classes that are not cloud but are not vegetation either brought the spread down by a factor of three.

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.

The spread is partly the fraction of the square that was visible, and partly what the visible part was.

What you get — from The data used on this page

127 scenes from 2025-06-02 to 2025-12-31

Few usable dates, clear fractions, per-date medians

1. Expect far fewer usable dates than scenes

127 scenes over the area
 29 left at least 30 pixels (1.2% of the square) outside the cloud classes
  6 left at least 90%
  3 were at least 90% clear of cloud, snow, dark area and unclassified

Three to six usable observations from 127 overpasses is what a 1 km² square in a maritime upland climate yielded over seven months. 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 square the twelve observations under 20% clear ranged from −0.067 to 0.804 — nearly the whole range of the series, from pixels that were mostly cloud edge.

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	cloud classes dropped		+ snow, dark area, unclassified
	observations	sd	observations	sd
any clear pixels	29	0.295	22	0.263
at least 20% clear	17	0.317	15	0.238
at least 50% clear	11	0.293	7	0.092
at least 80% clear	8	0.303	3	0.038
at least 90% clear	6	0.200	3	0.038

With only the cloud classes dropped, no threshold below 90% helped: 29 November was 87% "clear" and gave the series minimum, −0.092, with 66% of the square labelled unclassified. A clear-fraction threshold only means something once the class list behind it does.

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

In 2025, Sentinel-2A imaged tile 30UVD on the same days as Sentinel-2B, every tenth day, so 21 dates in this stack have two items. Overlapping tiles do the same thing to a polygon near a tile edge.

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

The square's NDVI fell from 0.762 on 12 July to 0.436 on 29 December, with 0.170 and −0.092 in between. Only part of that is the vegetation. On 19 November the square was 58.6% snow. On 27 December, at a 13.5° sun, 36.2% of it was classed as dark area — at that sun elevation, most likely terrain shadow. A smoother fitted through those points would present snow and shadow as a phenology curve.

The funnel from overpasses to usable observations is the real constraint on any time series.

A time series with its evidence attached

The data used on this page

The polygon is a 1 km square in the Llanberis community, 3 km north-west of Snowdon's summit — EPSG:32630, 425,880–426,880 E, 5,882,440–5,883,440 N — which the clear 12 July 2025 scene classifies entirely as vegetation. The scenes are every Sentinel-2 L2A item in Earth Search that covers it from 1 June to 31 December 2025; Example 1 reads the red, near-infrared and SCL windows over the square from each one.

import pystac_client
from pyproj import Transformer
from shapely.geometry import box, mapping

SQUARE = (425_880, 5_882_440, 426_880, 5_883_440)            # EPSG:32630
geometry = mapping(box(*SQUARE))
to_wgs84 = Transformer.from_crs(32630, 4326, always_xy=True)
bbox = [*to_wgs84.transform(*SQUARE[:2]), *to_wgs84.transform(*SQUARE[2:])]
items = sorted(pystac_client.Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=bbox, datetime="2025-06-01/2025-12-31").items(),
    key=lambda i: i.datetime)
print(len(items), "scenes from", items[0].datetime.date(), "to", items[-1].datetime.date())

127 scenes from 2025-06-02 to 2025-12-31

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 window_for(ds, aoi):
    """The window of ds covering aoi, a bounds tuple in the raster's CRS."""
    return rasterio.windows.from_bounds(*aoi, transform=ds.transform).round_offsets().round_lengths()


_scenes = {}


def read_scene(item, aoi):
    """SCL, red and NIR over aoi on the 20 m SCL grid; read once per item."""
    if item.id not in _scenes:
        with rasterio.open(item.assets["scl"].href) as ds:
            window = window_for(ds, aoi)
            scl = ds.read(1, window=window)
            transform = ds.window_transform(window)
        bands = []
        for name in ("red", "nir"):
            with rasterio.open(item.assets[name].href) as ds:
                bands.append(ds.read(1, window=window_for(ds, aoi), out_shape=scl.shape,
                                     resampling=Resampling.bilinear).astype("float32"))
        _scenes[item.id] = (scl, *bands, transform)
    return _scenes[item.id]


def ndvi_series(items, geometry, aoi, min_fraction=0.5, min_pixels=30, unusable=UNUSABLE):
    """One row per scene: date, clear fraction, NDVI statistics."""
    rows = []
    for item in items:
        try:
            scl, red, nir, transform = read_scene(item, aoi)
        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


rows = ndvi_series(items, geometry, SQUARE, min_fraction=0.0)

  127 scenes, 29 usable at 0% clear

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"])


daily = deduplicate_by_date(rows)
for r in daily:
    if r["ndvi_median"] is not None and r["clear_fraction"] >= 0.5:
        print(f"  {r['date']}  {r['clear_fraction']:6.1%} clear  NDVI {r['ndvi_median']:.3f}")

  removed 21 same-day duplicate observation(s)
  2025-07-12  100.0% clear  NDVI 0.762
  2025-08-16   75.2% clear  NDVI 0.749
  2025-09-22   95.2% clear  NDVI 0.705
  2025-10-13   99.8% clear  NDVI 0.670
  2025-11-17   74.8% clear  NDVI 0.540
  2025-11-19  100.0% clear  NDVI 0.170
  2025-11-29   86.8% clear  NDVI -0.092
  2025-12-19   81.8% clear  NDVI 0.053
  2025-12-27  100.0% clear  NDVI 0.546
  2025-12-29  100.0% clear  NDVI 0.436

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}")


threshold_sensitivity(rows)
low = [r for r in rows if r["ndvi_median"] is not None and r["clear_fraction"] < 0.2]
low_values = [r["ndvi_median"] for r in low]
print(f"  under 20% clear:  {len(low):3d} obs   {min(low_values):6.3f} to {max(low_values):6.3f}   "
      f"sd {np.std(low_values):.3f}")

  >=   0% clear:   29 obs   -0.092 to  0.804   sd 0.295
  >=  20% clear:   17 obs   -0.092 to  0.778   sd 0.317
  >=  50% clear:   11 obs   -0.092 to  0.778   sd 0.293
  >=  80% clear:    8 obs   -0.092 to  0.762   sd 0.303
  >=  90% clear:    6 obs    0.170 to  0.762   sd 0.200
  under 20% clear:   12 obs   -0.067 to  0.804   sd 0.246

Run this before choosing a threshold. If the series statistics change a lot across it, then the threshold is the most important parameter in your analysis and it belongs in the caption of every plot.

Example 4 — what "clear" leaves in

The drop list above keeps SCL classes 2 (dark area), 7 (unclassified) and 11 (snow). Over an upland square in November and December, those are not small classes:

from collections import Counter

for item in items:
    scl, *_ = read_scene(item, SQUARE)
    kept = ~np.isin(scl, UNUSABLE)
    if kept.mean() >= 0.8 and item.datetime.month >= 11:
        shares = Counter(scl[kept].ravel().tolist())
        print(f"  {item.datetime.date()}  " + "  ".join(
            f"{code}: {count / scl.size:5.1%}" for code, count in sorted(shares.items())))

STRICT = UNUSABLE + (2, 7, 11)          # also drop dark area, unclassified and snow
strict_rows = ndvi_series(items, geometry, SQUARE, min_fraction=0.0, unusable=STRICT)
threshold_sensitivity(strict_rows)

  2025-11-19  2:  1.4%  4:  6.0%  5: 34.0%  11: 58.6%
  2025-11-29  2: 11.8%  5:  7.0%  6:  2.0%  7: 66.0%
  2025-12-19  2: 35.1%  5: 24.9%  6:  0.1%  7: 21.7%
  2025-12-27  2: 36.2%  4: 39.9%  5: 23.9%
  2025-12-29  2: 42.3%  4:  5.4%  5: 52.3%
  127 scenes, 22 usable at 0% clear
  >=   0% clear:   22 obs   -0.103 to  0.810   sd 0.263
  >=  20% clear:   15 obs    0.022 to  0.778   sd 0.238
  >=  50% clear:    7 obs    0.513 to  0.778   sd 0.092
  >=  80% clear:    3 obs    0.670 to  0.762   sd 0.038
  >=  90% clear:    3 obs    0.670 to  0.762   sd 0.038

Five of the late dates were mostly something other than vegetation or cloud: snow on 19 November, unclassified pixels on 29 November, and a third or more dark area in all three December scenes. Dropped as well, they take the 50%-clear series from 11 observations with a standard deviation of 0.293 to 7 with 0.092 — and the drop list, not the threshold, did most of that.

Example 4's class shares beside Example 2's NDVI. The low values are snow, shadow and pixels the classifier would not name.

Partial visibility biases the result

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.067 to 0.804 while the three dates at least 90% clear of every doubtful class spanned 0.670 to 0.762. 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

Three to seven usable observations across seven months, unevenly spaced, 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 — three to six 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 late-season values are real measurements, but of snow, of shadow and of whatever the classifier would not name, as much as of vegetation. Only knowing the place, and reading the class shares behind each value, tells you which.

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.

All three produce a plausible curve. The clear fraction reveals the second; only the class shares reveal the third.

Clear fractions, duplicates and medians

Record the clear fraction with every observation. Without it a value is uninterpretable.
Deduplicate by date. Two satellites on one day, or overlapping tiles, produce two items for one date.
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.
Snow is class 11, dark area class 2. A cloud-only drop list keeps both; in an upland winter they were the largest classes in the "clear" pixels.

Internal links

How to calculate NDVI from Sentinel-2 in Python — the per-scene calculation
Cloud masking explained — why partial visibility is biased, not noisy
How to mask clouds in Sentinel-2 imagery with Python — the mask this depends on
How to build a cloud-free composite from many satellite scenes — the alternative when dates do not matter
Spectral indices explained — choosing an index other than NDVI
How to calculate zonal statistics in Python — the polygon summary in general
Radiometric levels explained — a processing change looks like an event
How to detect change between two satellite images in Python — two dates rather than many
How to extract phenology metrics from an NDVI time series — turning the series into numbers

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?

Partly cloudy dates, and — in winter — snow and shadow the cloud mask keeps. On the series measured here, requiring 50% clear of cloud alone left a standard deviation of 0.293; also dropping snow, dark area and unclassified pixels brought it to 0.092.

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² square near Snowdon, 127 scenes from June to December 2025 gave 29 with any clear pixels, 6 at least 90% clear of cloud, and 3 at least 90% clear of cloud, snow and shadow.

Why do I have two observations on the same date?

In 2025 Sentinel-2A and 2B imaged the same tile on the same day every ten days, and overlapping tiles do the same near a tile edge. Keep the one with the higher clear fraction.

Should I smooth or gap-fill the series?

Only after looking at the raw points and the class shares behind them. Smoothing encodes an assumed seasonal shape, and fitted through snow and shadow dates it presents them as phenology.

Mean or median within the polygon?

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