--- title: "NDVI time series in Python with Sentinel-2, odc-stac and xarray" description: "Sentinel-2 NDVI time series in Python with odc-stac, dask and xarray: SCL cloud masking, monthly composites, per-field curves and crop phenology." url: https://www.technicalgis.com/tutorials/ndvi-time-series-odc-stac/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-10-17 updated: 2026-10-17 language: en-GB --- # NDVI time series in Python with Sentinel-2, odc-stac and xarray > Sentinel-2 NDVI time series in Python with odc-stac, dask and xarray: SCL cloud masking, monthly composites, per-field curves and crop phenology. - Level: Intermediate - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/ndvi-time-series-odc-stac/) ## Key takeaways 1. odc.stac.load turns a STAC search into a lazy xarray datacube: 396 Sentinel-2 scenes from 2021 to 2025 became a 0.77 GB dask-backed cube in seconds, and nothing was downloaded until .compute(). 2. Read the small SCL band first, drop mostly cloudy dates, then compute NDVI only where it is needed. Here 245 of 396 dates had at least 20% clear pixels, and 67% of the remaining pixel-dates were clear. 3. Since processing baseline 04.00 (January 2022), Sentinel-2 L2A digital numbers include a +1000 offset. Forgetting to remove it biases NDVI low across most of the archive. 4. Smoothed per-field curves separate the crops clearly: in 2023 winter wheat peaked on 16 May, potato and sugar beet on 15 July and maize on 26 August. Sugar beet stayed green until about 28 October, five weeks after potato. ## Before you start - Comfortable with NumPy, pandas and the basics of xarray (the Sentinel-2 Lake Mead tutorial introduces STAC searches and Cloud Optimised GeoTIFFs) - Python 3.10+ with odc-stac, dask, xarray, pystac-client, planetary-computer, rasterio, geopandas, scipy, matplotlib and imageio - An internet connection: imagery is streamed from the Microsoft Planetary Computer and needs no account or API key **Short answer:** to build an NDVI time series in Python, search the Sentinel-2 L2A archive with `pystac-client`, then pass the items to **`odc.stac.load(..., chunks={"time": 1})`** to get a lazy, dask-backed **xarray** datacube with dimensions `time × y × x`. Mask cloud and shadow with the **SCL** band (keep classes 4–7), remove the +1000 offset from scenes processed with baseline 04.00 or later, and compute NDVI = (B08 − B04) / (B08 + B04). Then use xarray to aggregate: `ndvi.resample(time="1MS").median()` for monthly composites, a zonal mean per field for field curves, and a Savitzky–Golay filter to smooth them before reading off phenology dates. This tutorial follows **57 fields in Methwold Fen, Norfolk**, an area of deep, flat Fenland soils in East Anglia that grows winter wheat, sugar beet, potatoes and maize, through **five growing seasons (2021–2025)**. The field boundaries and 2023 crop labels come from the Rural Payments Agency’s **Crop Map of England (CROME)**. You will learn how to: 1. download crop labels for an area from the CROME 2023 web service and turn them into field polygons; 2. search **every** Sentinel-2 L2A scene of one tile for five years with STAC; 3. load them as a **lazy datacube with odc-stac and dask**, and see how much would be read; 4. build an **SCL cloud mask** and keep only usable dates; 5. compute **NDVI**, handling the processing-baseline offset; 6. make **monthly median composites** and small-multiple maps; 7. extract a **per-field NDVI time series** quickly with `np.bincount`; 8. smooth the curves with a **Savitzky–Golay** filter and estimate **green-up, peak and senescence** for each crop. All outputs below are real results from [`ndvi_time_series.py`](https://www.technicalgis.com/code/ndvi-time-series-odc-stac/), run on 17 October 2026 with Python 3.12, odc-stac 0.5.3, xarray 2026.9.0 and dask 2026.8.0. ## Why use odc-stac for a Sentinel-2 time series? A five-year time series of one Sentinel-2 tile is hundreds of scenes, each with a dozen bands stored as separate Cloud Optimised GeoTIFFs. You could open each file with rasterio and stack the arrays yourself, as in our [Lake Mead tutorial](https://www.technicalgis.com/tutorials/sentinel-2-stac-lake-mead/), but odc-stac does the hard parts for you: - it reprojects and resamples every item onto **one common pixel grid**, so scenes from different orbits and tiles line up; - it **mosaics** items from the same day (`groupby="solar_day"`); - with `chunks=`, it returns **dask arrays**, so building the cube costs almost nothing and only the pixels you actually compute are read, in parallel; - the result is a normal **xarray Dataset**, so `resample`, `groupby`, `where` and `median` work along the time axis. `stackstac` is a popular alternative with a similar idea. The code below would need only small changes to use it. ## Step 1: field boundaries and crop labels from CROME CROME classifies England into hexagons of about 0.4 ha each year, using Sentinel-1 and Sentinel-2 imagery and a machine-learning model trained on field inspections. The 2023 map is published as an ArcGIS feature service. The script downloads the hexagons inside the bounding box in pages of 2000, caches them in a GeoPackage, and then dissolves touching hexagons of the same crop into blocks: `ndvi_time_series.py` ```python BBOX = (0.40, 52.50, 0.46, 52.54) # Methwold Fen, Norfolk (lon/lat): about 4 x 4.5 km of Fenland fields CRS = "EPSG:32631" # UTM zone 31N, the native grid of Sentinel-2 tile 31UCU CROME = ("https://services-eu1.arcgis.com/KB6uNVj5ZcJr7jUP/ArcGIS/rest/services/" "Crop_Map_of_England_2023/FeatureServer/0/query") CROPS = {"AC66": "Winter wheat", "AC03": "Sugar beet", "AC44": "Potato", "AC17": "Maize", "PG01": "Grass"} MIN_FIELD_HA = 3 while True: # the service returns at most 2000 hexagons per request: page through them r = requests.get(CROME, timeout=180, params={ "where": "1=1", "geometry": geom, "geometryType": "esriGeometryEnvelope", "inSR": 4326, "spatialRel": "esriSpatialRelIntersects", "outFields": "LUCODE,PROB,CROMEID", "outSR": 4326, "orderByFields": "OBJECTID", "resultOffset": offset, "resultRecordCount": 2000, "f": "geojson"}) page = r.json() parts.append(gpd.GeoDataFrame.from_features(page["features"], crs=4326)) if len(page["features"]) < 2000: break offset += 2000 blocks = hexes[hexes["LUCODE"].isin(CROPS)].dissolve("LUCODE").explode(index_parts=False).reset_index() blocks["crop"] = blocks["LUCODE"].map(CROPS) blocks["geometry"] = blocks.buffer(-10) # shrink 10 m so edge pixels (hedges, tracks) drop out blocks = blocks[~blocks.is_empty].explode(index_parts=False) blocks["area_ha"] = blocks.area / 1e4 blocks = blocks[blocks["area_ha"] >= MIN_FIELD_HA] ``` `Output — CROME 2023 hexagons and crop blocks` ```text 4859 hexagons of 0.42 ha; median classification probability 0.61 LUCODE AC66 Winter wheat 1750 AC17 Maize 897 AC03 Sugar beet 542 AC44 Potato 387 PG01 Grass 286 FA01 229 AC01 121 AC67 117 WO12 106 LG20 104 median probability by crop: {'Sugar beet': 0.68, 'Maize': 0.37, 'Potato': 0.74, 'Winter wheat': 0.95, 'Grass': 0.26} === crop blocks (fields) used for time series === count sum median crop Grass 4 38.5 10.8 Maize 14 236.2 5.7 Potato 13 85.5 4.7 Sugar beet 14 148.4 7.5 Winter wheat 12 637.2 23.6 ``` That leaves 57 blocks of at least 3 ha. Note the `PROB` column: CROME is itself a model output, and its confidence varies a lot by crop. It is very sure about winter wheat (median 0.95) but much less sure about maize (0.37) and grass (0.26), so treat those labels with care. ![Two maps of the same 4.2 by 4.6 km area. Left: a Sentinel-2 true-colour image from 10 June 2023 showing a patchwork of rectangular green and brown fields. Right: the same area in grey with the crop blocks coloured by crop, mostly large winter wheat blocks in gold, with sugar beet, potato, maize and a few grass blocks.](https://www.technicalgis.com/_astro/01-study-area.Cef4wyQz_1DPpxv.webp) _The study area in Methwold Fen on 10 June 2023 (left) and the 57 CROME 2023 crop blocks of at least 3 ha used for the time series (right). Blocks are dissolved hexagons, so one block can span several neighbouring fields with the same crop, as the large wheat blocks do._ ## Step 2: search five years of Sentinel-2 with STAC The Microsoft Planetary Computer hosts the full Sentinel-2 L2A archive and its STAC API needs no key; `planetary_computer.sign_inplace` adds a short-lived token to each asset URL. Filtering on `s2:mgrs_tile` keeps one tile, so every scene has the same footprint. A few dates exist twice, an original and a reprocessed version, so the script keeps the newest processing baseline: `ndvi_time_series.py` ```python catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace) search = catalog.search(collections=["sentinel-2-l2a"], bbox=BBOX, datetime="2021-01-01/2025-12-31", query={"eo:cloud_cover": {"lt": 80}, "s2:mgrs_tile": {"eq": "31UCU"}}) items = list(search.items()) best = {} # a few dates exist twice (original and reprocessed): keep the newest processing baseline for it in items: key = it.datetime.date() if key not in best or it.properties["s2:processing_baseline"] > best[key].properties["s2:processing_baseline"]: best[key] = it items = sorted(best.values(), key=lambda i: i.datetime) ``` `Output — STAC search` ```text 396 Sentinel-2 L2A scenes of tile 31UCU with < 80% cloud, 03 Jan 2021 to 31 Dec 2025 {2021: 75, 2022: 81, 2023: 73, 2024: 63, 2025: 104} ``` The tile-level `eo:cloud_cover` filter is loose on purpose. It describes the whole 110 km tile, not our 4 km area, so a scene that is 70% cloudy overall can still be clear over the fields. The real cloud filtering happens per pixel in Step 4. The jump in 2025 is mostly Sentinel-2C, which began routine operations in January 2025 while Sentinel-2A continued to take some images. ## Step 3: load a lazy datacube with odc-stac This is the core of the tutorial. `odc.stac.load` builds the datacube from the STAC items: `ndvi_time_series.py` ```python odc.stac.configure_rio(cloud_defaults=True) # sensible GDAL settings for reading COGs over HTTP ds = odc.stac.load( items, bands=["B02", "B03", "B04", "B08", "SCL"], bbox=BBOX, crs=CRS, resolution=10, chunks={"time": 1, "x": 1024, "y": 1024}, # one chunk per scene and band: nothing is read yet groupby="solar_day", dtype="uint16", nodata=0, ) ds = ds.assign_coords(time=ds.time.dt.floor("D")) # one scene per day: keep just the date print(ds) print(f"{ds.nbytes / 1e9:.2f} GB if loaded in full, in {ds['B04'].data.npartitions} chunks per band") ``` `Output — lazy datacube` ```text Size: 772MB Dimensions: (y: 460, x: 424, time: 396) Coordinates: * y (y) float64 4kB 5.824e+06 5.824e+06 ... 5.82e+06 5.82e+06 * x (x) float64 3kB 3.235e+05 3.235e+05 ... 3.277e+05 3.277e+05 * time (time) datetime64[us] 3kB 2021-01-03 2021-01-06 ... 2025-12-31 spatial_ref int32 4B 32631 Data variables: B02 (time, y, x) uint16 154MB dask.array B03 (time, y, x) uint16 154MB dask.array B04 (time, y, x) uint16 154MB dask.array B08 (time, y, x) uint16 154MB dask.array SCL (time, y, x) uint16 154MB dask.array 0.77 GB if loaded in full, in 396 chunks per band ``` The call returns almost immediately because no pixels have been read: each band is a dask array of 396 chunks, one per date, and each chunk knows which windowed reads of which COGs it needs. A few points about the arguments: - **`crs` and `resolution`** fix the output grid. Using the tile’s own UTM zone at 10 m avoids resampling the 10 m bands. The 20 m SCL band is upsampled with nearest neighbour. - **`chunks={"time": 1, ...}`** gives one chunk per scene. That suits per-date operations such as cloud masking. For per-pixel operations along time, such as a median over a whole season, rechunk later to `{"time": -1}`. - **`dtype="uint16", nodata=0`** keep the raw digital numbers. Converting to float32 would double the memory before you need it. - **`groupby="solar_day"`** mosaics items acquired on the same day, which matters when your area spans two tiles or two adjacent orbit segments. ## Step 4: SCL cloud masking The scene classification layer labels each pixel. The script keeps vegetation (4), bare soil (5), water (6) and unclassified (7), and treats cloud shadow (3), medium and high cloud (8, 9) and cirrus (10) as cloudy. SCL usually misses the soft edges of clouds, so the cloudy pixels are grown by three pixels (30 m). Reading only the SCL band first is cheap, and it tells us which dates are worth reading in full: `ndvi_time_series.py` ```python KEEP_SCL = [4, 5, 6, 7] # vegetation, bare soil, water, unclassified MIN_CLEAR = 0.2 # skip scenes where less than 20% of the area is clear dask.config.set(scheduler="threads", num_workers=16) scl = ds["SCL"].compute() cloudy = np.isin(scl.values, [3, 8, 9, 10]) # cloud shadow, cloud (medium, high), cirrus grown = ndimage.binary_dilation(cloudy, iterations=3, axes=(1, 2)) # grow by 30 m: edges are often missed clear = scl.isin(KEEP_SCL) & ~xr.DataArray(grown, coords=scl.coords, dims=scl.dims) clear_frac = clear.mean(("y", "x")).to_series() ``` `Output — clear-sky fraction per scene` ```text read SCL for 396 dates in 4 s === clear-sky fraction per scene === 245 of 396 dates have at least 20% clear pixels share of usable scenes by month (%): {1: 66, 2: 58, 3: 61, 4: 59, 5: 50, 6: 63, 7: 59, 8: 68, 9: 55, 10: 61, 11: 74, 12: 65} ``` `axes=(1, 2)` dilates in space only, not across dates. It needs SciPy 1.15 or later. ![Top: a scatter plot of the clear-sky percentage of each of 396 Sentinel-2 dates from 2021 to 2025. Many dates are near 0% or 100%, with a dashed red line at 20%. Bottom: a bar chart of usable dates per month, averaging between about three and five and a half per year in each month.](https://www.technicalgis.com/_astro/02-clear-observations.BfD_TP85_dTgtL.webp) _Clear-sky share of the study area on each Sentinel-2 date (top) and the mean number of usable dates per calendar month (bottom). Even in East Anglia, one of the driest parts of the UK, about 40% of scenes are unusable, and February, May, September and December average about three usable dates a year or fewer._ ## Step 5: compute NDVI, with the processing-baseline offset Sentinel-2 L2A digital numbers are reflectance × 10,000. From **processing baseline 04.00**, introduced on 25 January 2022, ESA adds an offset of **+1000** so that negative reflectances can be stored. On the Planetary Computer, most older scenes have not been reprocessed, so a five-year cube mixes both conventions. The script builds a per-date offset from the items’ `s2:processing_baseline` property: `ndvi_time_series.py` ```python baseline = pd.Series({pd.Timestamp(i.datetime.date()): float(i.properties["s2:processing_baseline"]) for i in items}) offset = xr.DataArray(np.where(baseline.reindex(ds.time.dt.floor("D").values).values >= 4.0, 1000, 0), dims="time", coords={"time": ds.time}) def ndvi_of(sub: xr.Dataset) -> xr.DataArray: """NDVI from L2A digital numbers. Baseline 04.00+ (from 25 Jan 2022) stores reflectance x 10000 + 1000.""" off = offset.sel(time=sub.time).astype("float32") red = (sub["B04"].astype("float32") - off) / 10000 nir = (sub["B08"].astype("float32") - off) / 10000 return ((nir - red) / (nir + red)).where((sub["B04"] > 0) & (sub["B08"] > 0)) keep_times = clear_frac.index[clear_frac >= MIN_CLEAR] lazy = ndvi_of(ds.sel(time=keep_times)).where(clear.sel(time=keep_times)) ndvi = lazy.compute().rename("ndvi") ``` `Output — NDVI cube` ```text lazy NDVI: float32, {'time': 245, 'y': 460, 'x': 424}, 245 chunks computed NDVI for 245 dates in 14 s === NDVI cube === {'time': 245, 'y': 460, 'x': 424}, 67.4% of pixel-dates clear NDVI percentiles 5/50/95: [0.189 0.559 0.898] ``` `ds.sel(time=keep_times)` drops the cloudy dates **before** anything is computed, so dask reads the red and NIR bands for 245 dates instead of 396. The whole step took 14 seconds. The script caches the result in a NetCDF file, so later runs skip the download. ## Step 6: monthly median composites With the time series in xarray, monthly composites are one line. A median ignores NaN (cloud) and is robust to a single bad observation, such as cloud shadow that SCL missed: `ndvi_time_series.py` ```python monthly = ndvi.resample(time="1MS").median() n_obs = ndvi.notnull().resample(time="1MS").sum() cover = monthly.notnull().mean(("y", "x")) monthly.astype("float32").to_netcdf(OUT / "ndvi_monthly_median_2021_2025.nc") ``` `Output — monthly median composites (2023)` ```text 60 months; 47 have NDVI for at least 95% of pixels, 1 have none Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec median NDVI 0.44 0.42 0.49 0.54 0.69 0.73 0.74 0.63 0.55 0.43 0.51 0.41 clear obs/pixel 3 3 2 1 3 4 2 3 2 3 4 1 cover % 100 100 100 87 100 100 100 99 100 100 100 91 ``` Small multiples, one map per month, are the clearest way to show a datacube on a static page. The script also writes the same twelve frames as an animated GIF (`output/ndvi_2023.gif`) with `imageio`, which is handy for slides and social media. ![Twelve small NDVI maps of the study area, January to December 2023, coloured from red (low) through yellow to blue (high). Some fields turn blue in March to May, others only in June to August, and a few in September to November. April and December have grey cloud gaps.](https://www.technicalgis.com/_astro/03-monthly-ndvi-2023.-ISBsakf_FBtEY.webp) _Monthly median NDVI in 2023 with the Spectral colour scale. Winter wheat fields are already blue (green and leafy) by March and turn orange after harvest in August, while sugar beet, potato and maize fields only turn blue in June and July. Grey areas had no clear observation in that month._ ## Step 7: a per-field NDVI time series To get one value per field and date, rasterise the field polygons onto the same grid as the cube, then sum the NDVI of each field’s pixels with `np.bincount`. That is much faster than a zonal-statistics loop over polygons, and it lets you require that most of a field is clear before using its value: `ndvi_time_series.py` ```python gbox = ds.odc.geobox # the cube's grid: shape, CRS and affine transform labels = rasterio.features.rasterize(zip(blocks.geometry, blocks["field_id"]), out_shape=gbox.shape, transform=gbox.affine, fill=0, dtype="int32") n_pix = np.bincount(labels.ravel(), minlength=len(blocks) + 1) MIN_FIELD_CLEAR = 0.8 # a field counts as observed on a date only if 80% of its pixels are clear rows = [] for t in range(ndvi.sizes["time"]): a = ndvi.isel(time=t).values.ravel() ok = np.isfinite(a) s = np.bincount(labels.ravel()[ok], weights=a[ok], minlength=len(blocks) + 1) c = np.bincount(labels.ravel()[ok], minlength=len(blocks) + 1) frac = c / np.maximum(n_pix, 1) for f in np.nonzero(frac[1:] >= MIN_FIELD_CLEAR)[0] + 1: rows.append((pd.Timestamp(ndvi.time.values[t]).normalize(), f, s[f] / c[f], frac[f])) ts = pd.DataFrame(rows, columns=["date", "field_id", "ndvi", "clear_frac"]) ``` `Output — per-field time series (field-dates per crop and year)` ```text 8159 field-dates for 57 fields; median 144 clear observations per field over five years date 2021 2022 2023 2024 2025 crop Grass 93 104 105 91 196 Maize 300 348 364 317 664 Potato 299 336 360 303 632 Sugar beet 280 359 359 322 666 Winter wheat 261 274 314 250 562 ``` `import odc.geo.xr` registers the `.odc` accessor on xarray objects, which gives the cube’s geobox. The result is a tidy table, also written to `field_ndvi_timeseries.csv`, with about 29 clear observations per field per year. ## Step 8: smooth the curves and estimate phenology Raw field means are still noisy: thin cloud and shadow that SCL missed pull NDVI **down**, almost never up. A standard fix (Chen et al., 2004) is an **upper-envelope Savitzky–Golay filter**: interpolate onto a regular 5-day grid, fit the filter, replace observations far below the curve with the fitted value, and repeat: `ndvi_time_series.py` ```python STEP = "5D" # resample to a regular 5-day grid before filtering WINDOW, ORDER = 9, 2 # 9 x 5 days = 45-day window, quadratic fit def smooth(series: pd.Series, start: str = "2021-01-01", end: str = "2025-12-31") -> pd.Series: """Upper-envelope Savitzky-Golay smoothing (after Chen et al., 2004) on a regular 5-day grid.""" grid = pd.date_range(start, end, freq=STEP) s = series.groupby(level=0).mean() for _ in range(3): # observations far below the curve are probably cloud or shadow: pull them up filled = s.reindex(s.index.union(grid)).interpolate("time", limit_direction="both").reindex(grid) fit = pd.Series(savgol_filter(filled.values, WINDOW, ORDER), index=grid) on_fit = fit.reindex(s.index, method="nearest") s = s.where(s >= on_fit - 0.05, on_fit) return fit.clip(-0.2, 1) ``` ![Four line charts from 2021 to 2025, one each for a winter wheat, sugar beet, potato and maize field. Grey dots are observations and a coloured line is the smoothed curve. Each field shows one main peak per year, but the timing and shape of the peak change from year to year.](https://www.technicalgis.com/_astro/04-field-time-series.Cc53cmBx_Z1EO5cx.webp) _Five years of NDVI for the most typical field of each 2023 crop (the field whose 2023 curve is closest to its crop's median). The 2023 peaks follow the crop: May for wheat, July for beet and potato, August for maize. In other years the curves differ (field 47 peaks in June or July rather than May, for example) because crops are rotated, and the CROME labels describe 2023 only._ For phenology, the script finds the 2023 peak and the dates where the smoothed curve crosses **50% of the seasonal amplitude** on the way up (green-up) and down (senescence or harvest). Thresholds of 20–50% are all common; 50% is the most robust to noise. `ndvi_time_series.py` ```python def phenology(curve: pd.Series, year: int) -> dict: """Green-up, peak and senescence dates where the curve crosses 50% of the seasonal amplitude.""" c = curve[f"{year - 1}-10-01":f"{year}-12-31"] peak = c[f"{year}-03-01":f"{year}-10-31"].idxmax() lo_before, lo_after = c[:peak].min(), c[peak:].min() up = c[:peak][c[:peak] < lo_before + 0.5 * (c[peak] - lo_before)] down = c[peak:][c[peak:] < lo_after + 0.5 * (c[peak] - lo_after)] return {"green_up": up.index[-1] if len(up) else pd.NaT, "peak": peak, "peak_ndvi": c[peak], "senescence": down.index[0] if len(down) else pd.NaT, "amplitude": c[peak] - min(lo_before, lo_after)} curves = {f: smooth(g.set_index("date")["ndvi"]) for f, g in ts.groupby("field_id") if len(g) >= 40} pheno = pd.DataFrame([dict(field_id=f, **phenology(c, 2023)) for f, c in curves.items()]) ``` `Output — 2023 phenology by crop (median of fields)` ```text fields green_up peak peak_ndvi senescence crop Winter wheat 12 25 Feb 16 May 0.87 20 Jul Grass 4 01 Apr 15 Jun 0.80 16 Aug Potato 13 10 Jun 15 Jul 0.91 23 Sep Sugar beet 14 05 Jun 15 Jul 0.91 28 Oct Maize 14 25 Jun 26 Aug 0.83 15 Oct ``` The results match the Fenland crop calendar well: - **Winter wheat**, sown in autumn, is already green over winter, greens up strongly in late February and peaks in mid-May. NDVI falls as the crop ripens and is harvested in late July and August. - **Potato and sugar beet** are planted in spring, cover the ground in June and both peak in mid-July at an NDVI of 0.91. Potato haulm is destroyed before lifting, so potato NDVI falls in September; sugar beet stays green until it is lifted, here from late October. - **Maize** is sown in May, greens up last (late June) and peaks in late August. ![A line chart from October 2022 to December 2023 with one curve per crop and a shaded interquartile band. Winter wheat peaks first in May, potato and sugar beet rise steeply in June to peak in July, maize peaks in late August, and grass stays between 0.5 and 0.8 all year.](https://www.technicalgis.com/_astro/05-phenology-by-crop.tl7axacA_ZtqaNH.webp) _The 2023 phenology curves: median smoothed NDVI of the fields of each crop, with the interquartile range shaded and a dot at the median peak date. The four arable crops are separated mainly by timing, not by peak height._ Finally, a heatmap of the monthly maximum NDVI of every field shows crop rotation. Each row is a field and each column a month. In 2023 almost every potato and sugar beet row is dark blue from late June to September, and the winter wheat rows peak earlier, in spring. In the other years the same rows are far less consistent, because each field carried a different crop. ![A heatmap with one row per field and one column per month from 2021 to 2025, grouped into maize, potato, sugar beet and winter wheat by their 2023 crop. Coloured from red (low NDVI) to blue (high). In 2023 each group has its blue period at a different time; in other years the patterns are mixed.](https://www.technicalgis.com/_astro/06-crop-rotation-heatmap.HMIU469F_Z1zBKa0.webp) _Monthly maximum NDVI of every field, grouped by the 2023 crop. The groups are consistent in 2023 but not in other years, because Fenland farms rotate their crops. CROME labels are for one year only._ ## Pitfalls and limitations - **The labels are a model too.** CROME is produced by a machine-learning classifier, and its own confidence is low for maize (median 0.37) and grass (0.26) here. Some “maize” blocks may be another spring crop. For research, use field-inspection data or the farm’s own records. - **Blocks are not exactly fields.** Dissolving hexagons merges neighbouring fields with the same crop, and the hexagon edges cut across field boundaries. The 10 m inward buffer helps, but some blocks still include tracks and margins. Real field boundaries (for example from the Rural Payments Agency’s land parcels, where you have access) are better. - **One year of labels.** The crop labels apply to 2023 only. The 2021, 2022, 2024 and 2025 curves belong to whatever was grown then. - **SCL misses some cloud.** Thin cirrus, haze and shadow from small clouds get through; the upper-envelope filter and the 80% field-clear rule reduce, but do not remove, the effect. Low outliers in the field plots are mostly missed cloud. For better masks, try s2cloudless or the OmniCloudMask model. - **Gaps change the dates.** Phenology dates are only as good as the observations around them. A cloudy fortnight at green-up can shift the 50% crossing by a week or more, and the 45-day filter window smooths very quick changes such as harvest. - **One tile.** The script uses only tile 31UCU. Neighbouring tiles and orbits add more dates, but they need care with the processing baseline and different viewing angles. - **NDVI saturates.** Above an NDVI of about 0.85, dense canopies look alike. For biomass or yield work, try the red-edge index NDRE or EVI. ## Frequently asked questions ### How do I create an NDVI time series in Python? Search Sentinel-2 L2A scenes with `pystac-client`, load them into an xarray datacube with `odc.stac.load` (or `stackstac`), mask cloud with the SCL band, compute NDVI from bands B08 and B04 after removing the +1000 offset in newer scenes, and reduce over space: a zonal mean per field or a median over an area. Then smooth with a Savitzky–Golay filter. The full, runnable script is linked above. ### What is the difference between odc-stac and stackstac? Both turn STAC items into a dask-backed xarray object. odc-stac returns an `xarray.Dataset` with one variable per band, groups items by day, and supports per-band resampling and data types. stackstac returns one `DataArray` with a `band` dimension and is a little more minimal. Both are maintained and fast; odc-stac is used by Digital Earth Australia and Digital Earth Africa. ### How do I mask clouds in Sentinel-2 with Python? For L2A data, the simplest method is the SCL band: keep vegetation, bare soil, water and unclassified (classes 4–7) and mask cloud shadow, cloud and cirrus (3, 8, 9 and 10). Grow the cloud mask by a few pixels, because SCL often misses cloud edges. For better accuracy, use a machine-learning mask such as s2cloudless (the `S2_CLOUD_PROBABILITY` product) or OmniCloudMask. ### Why do my Sentinel-2 NDVI values drop in 2022? Probably because of the processing-baseline offset. Since baseline 04.00 (25 January 2022), L2A digital numbers have 1000 added. If you divide by 10,000 without first subtracting 1000, red and NIR are both 0.1 too high, which lowers NDVI. Read `s2:processing_baseline` from each STAC item and subtract the offset where it is 04.00 or later, as in Step 5. ### How many clear Sentinel-2 images will I get per year? It depends on the climate. In this East Anglian example, 245 of 396 scenes (62%) over five years had at least 20% of the area clear, about 49 dates per year, and the median field had 144 clear observations over five years. Western Britain and Ireland get fewer. Sentinel-2C, which began routine operations in 2025, increased the number of scenes noticeably here. ### Can I use Landsat as well? Yes. Harmonised Landsat Sentinel-2 (HLS) data provides Landsat 8 and 9 and Sentinel-2 on a common 30 m grid with matching bands, which roughly doubles the number of observations. On the Planetary Computer, `landsat-c2-l2` can be loaded with the same `odc.stac.load` call, but you need to apply its own scale and offset and its QA cloud mask. ## Next steps - Use the phenology curves as features to **classify crops**, as in our companion tutorial [land cover classification with a random forest](https://www.technicalgis.com/tutorials/land-cover-classification-random-forest/). - Swap the 4 km area for a whole farm or county and run the same code on a `dask.distributed` cluster. - Learn the STAC basics on a single date in [Measuring Lake Mead’s shrinking shoreline with Sentinel-2 and STAC](https://www.technicalgis.com/tutorials/sentinel-2-stac-lake-mead/), or apply time-series thinking to water in [reservoir water quality with Sentinel-2](https://www.technicalgis.com/tutorials/reservoir-water-quality-sentinel-2/). - When cloud blocks optical data, track crops and flooding with radar instead: [Sentinel-1 SAR flood mapping in Python](https://www.technicalgis.com/tutorials/sentinel-1-sar-flood-mapping/). - New to vector data in Python? Start with [GeoPandas: getting started](https://www.technicalgis.com/tutorials/geopandas-getting-started/). - When cloud leaves gaps in the NDVI record, follow wheat and oilseed rape with radar instead in [Sentinel-1 crop monitoring in Python](https://www.technicalgis.com/tutorials/sentinel-1-crop-monitoring/). ## Data and licences - **Sentinel-2 L2A**: contains modified Copernicus Sentinel data 2021–2025, processed by ESA. Accessed as Cloud Optimised GeoTIFFs from the [Microsoft Planetary Computer](https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a) on 17 October 2026. - **Crop Map of England (CROME) 2023**: © Rural Payments Agency, licensed under the Open Government Licence v3.0. Accessed through the [CROME 2023 feature service](https://services-eu1.arcgis.com/KB6uNVj5ZcJr7jUP/ArcGIS/rest/services/Crop_Map_of_England_2023/FeatureServer) on 17 October 2026. - References: Chen, J. et al. (2004), “A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter”, _Remote Sensing of Environment_ 91(3–4); Rouse, J. W. et al. (1974), “Monitoring vegetation systems in the Great Plains with ERTS”, NASA SP-351. - Code: [`ndvi_time_series.py`](https://www.technicalgis.com/code/ndvi-time-series-odc-stac/), tested with Python 3.12.0, odc-stac 0.5.3, xarray 2026.9.0, dask 2026.8.0, pystac-client 0.9.0, planetary-computer 1.0.0 and scipy 1.18.1. ## Key terms - **NDVI**: Normalised Difference Vegetation Index, (NIR − red) / (NIR + red), using Sentinel-2 bands B08 and B04. Dense green vegetation scores about 0.8–0.9, bare soil 0.1–0.3 and water below 0. - **odc-stac**: A Python library from the Open Data Cube project that loads STAC items into an xarray Dataset on a common pixel grid, optionally as lazy dask arrays. - **Datacube**: A stack of co-registered images with dimensions such as time, y and x, so that a pixel's whole history can be read as one array. - **SCL**: The Sentinel-2 L2A scene classification layer: a 20 m band that labels each pixel as cloud, cloud shadow, vegetation, bare soil, water, snow and so on. - **Savitzky–Golay filter**: A smoothing filter that fits a low-order polynomial in a moving window. It keeps the height and timing of peaks better than a moving average. - **Phenology**: The timing of the stages in a plant's life cycle, such as green-up, peak growth and senescence, here estimated from the NDVI curve. ## Data for this tutorial - [Microsoft Planetary Computer](https://www.technicalgis.com/data-sources/#src-microsoft-planetary-computer) · Sentinel-2 L2A as Cloud Optimised GeoTIFFs via STAC - [UK data sources](https://www.technicalgis.com/data-sources/united-kingdom/#src-data-gov-uk) · Crop Map of England (CROME) 2023, Rural Payments Agency Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/ndvi-time-series-odc-stac/