Source code: NDVI time series in Python with Sentinel-2, odc-stac and xarray

"""NDVI time series with Sentinel-2, STAC and odc-stac: five growing seasons of fields in Methwold Fen, Norfolk.

Searches the Microsoft Planetary Computer STAC API for every Sentinel-2 L2A scene of one tile from 2021
to 2025, loads them lazily as a dask-backed xarray datacube with odc-stac, masks cloud and shadow with
the scene classification (SCL) band, computes NDVI, builds monthly median composites, extracts per-field
time series for crop blocks from the Rural Payments Agency's Crop Map of England (CROME) 2023, smooths
them with a Savitzky-Golay filter and derives simple phenology dates for each crop.
Run with:  python ndvi_time_series.py
"""
import json
import platform
import time
import warnings
from importlib.metadata import version
from pathlib import Path

import dask
import geopandas as gpd
import imageio.v3 as iio
import matplotlib

matplotlib.use("Agg")
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import odc.geo.xr  # noqa: F401  (adds the .odc accessor to xarray objects)
import odc.stac
import pandas as pd
import planetary_computer
import pystac_client
import rasterio.features
import requests
import xarray as xr
from matplotlib.patches import Patch
from scipy import ndimage
from scipy.signal import savgol_filter

warnings.filterwarnings("ignore", category=RuntimeWarning)   # all-NaN medians in fully cloudy months
HERE = Path(__file__).parent
OUT, DATA = HERE / "output", HERE / "data"
for d in (OUT, DATA):
    d.mkdir(exist_ok=True)
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
pd.set_option("display.width", 140)

print("=== versions ===")
print(f"python {platform.python_version()} | odc-stac {odc.stac.__version__} | xarray {xr.__version__} | "
      f"dask {dask.__version__} | pystac-client {pystac_client.__version__} | "
      f"planetary-computer {version('planetary-computer')} | numpy {np.__version__} | pandas {pd.__version__}")


def section(title: str) -> None:
    print(f"\n=== {title} ===")


# --- Step 1: crop blocks from the Crop Map of England (CROME) 2023 ---------------------------------
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"}
COLOURS = {"Winter wheat": "#d4a017", "Sugar beet": "#c21f3a", "Potato": "#7b3fa0", "Maize": "#2f9e44",
           "Grass": "#1d6fa3"}
MIN_FIELD_HA = 3

crome_file = DATA / "crome_2023.gpkg"
if not crome_file.exists():
    geom = json.dumps({"xmin": BBOX[0], "ymin": BBOX[1], "xmax": BBOX[2], "ymax": BBOX[3],
                       "spatialReference": {"wkid": 4326}})
    parts, offset = [], 0
    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"})
        r.raise_for_status()
        page = r.json()
        parts.append(gpd.GeoDataFrame.from_features(page["features"], crs=4326))
        if len(page["features"]) < 2000:
            break
        offset += 2000
    pd.concat(parts, ignore_index=True).to_file(crome_file, driver="GPKG")
hexes = gpd.read_file(crome_file).to_crs(CRS)
section("CROME 2023 hexagons")
print(f"{len(hexes)} hexagons of {hexes.area.median() / 1e4:.2f} ha; "
      f"median classification probability {hexes['PROB'].median():.2f}")
print(hexes["LUCODE"].value_counts().head(10).rename(lambda c: f"{c} {CROPS.get(c, '')}").to_string())

# Dissolve touching hexagons of the same crop into blocks: usually one field, sometimes neighbours with the same crop.
print("median probability by crop:", hexes[hexes["LUCODE"].isin(CROPS)].groupby("LUCODE")["PROB"].median()
      .rename(CROPS).round(2).to_dict())
blocks = hexes[hexes["LUCODE"].isin(CROPS)].dissolve("LUCODE").explode(index_parts=False).reset_index()[["LUCODE", "geometry"]]
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].sort_values(["crop", "area_ha"], ascending=[True, False])
blocks = blocks.reset_index(drop=True)
blocks["field_id"] = np.arange(1, len(blocks) + 1)
section("crop blocks (fields) used for time series")
print(blocks.groupby("crop")["area_ha"].agg(["count", "sum", "median"]).round(1).to_string())

# --- Step 2: search the Sentinel-2 L2A archive ------------------------------------------------------
STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
PERIOD = "2021-01-01/2025-12-31"
catalog = pystac_client.Client.open(STAC_URL, modifier=planetary_computer.sign_inplace)
search = catalog.search(collections=["sentinel-2-l2a"], bbox=BBOX, datetime=PERIOD,
                        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)
section("STAC search")
print(f"{len(items)} Sentinel-2 L2A scenes of tile 31UCU with < 80% cloud, "
      f"{items[0].datetime:%d %b %Y} to {items[-1].datetime:%d %b %Y}")
print(pd.Series([i.datetime.year for i in items]).value_counts().sort_index().to_dict())

# --- Step 3: a lazy datacube with odc-stac -----------------------------------------------------------
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
section("lazy datacube")
print(ds)
print(f"{ds.nbytes / 1e9:.2f} GB if loaded in full, in {ds['B04'].data.npartitions} chunks per band")

# --- Step 4: SCL cloud mask, computed first because it is small -----------------------------------
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_file = DATA / "scl.nc"
if scl_file.exists():
    scl = xr.load_dataset(scl_file)["SCL"]
else:
    t0 = time.time()
    scl = ds["SCL"].compute()
    scl.to_netcdf(scl_file)
    print(f"read SCL for {scl.sizes['time']} dates in {time.time() - t0:.0f} s")
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()
section("clear-sky fraction per scene")
print(f"{(clear_frac >= MIN_CLEAR).sum()} of {len(clear_frac)} dates have at least {MIN_CLEAR:.0%} clear pixels")
by_month = clear_frac.groupby(clear_frac.index.month).apply(lambda s: (s >= MIN_CLEAR).mean())
print("share of usable scenes by month (%):", (by_month * 100).round().astype(int).to_dict())

# --- Step 5: NDVI, with the processing-baseline offset ----------------------------------------------
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]
ndvi_file = DATA / "ndvi.nc"
if ndvi_file.exists():
    ndvi = xr.load_dataset(ndvi_file)["ndvi"]
else:
    t0 = time.time()
    lazy = ndvi_of(ds.sel(time=keep_times)).where(clear.sel(time=keep_times))
    print(f"lazy NDVI: {lazy.dtype}, {dict(lazy.sizes)}, {lazy.data.npartitions} chunks")
    ndvi = lazy.compute().rename("ndvi")
    ndvi.to_netcdf(ndvi_file)
    print(f"computed NDVI for {ndvi.sizes['time']} dates in {time.time() - t0:.0f} s")
section("NDVI cube")
print(f"{dict(ndvi.sizes)}, {float(ndvi.notnull().mean()):.1%} of pixel-dates clear")
print(f"NDVI percentiles 5/50/95: {np.nanpercentile(ndvi, [5, 50, 95]).round(3)}")

# --- Step 6: monthly median composites ------------------------------------------------------------------
monthly = ndvi.resample(time="1MS").median()
n_obs = ndvi.notnull().resample(time="1MS").sum()
cover = monthly.notnull().mean(("y", "x"))
section("monthly median composites")
print(f"{monthly.sizes['time']} months; {int((cover >= 0.95).sum())} have NDVI for at least 95% of pixels, "
      f"{int((cover == 0).sum())} have none")
m23 = monthly.sel(time="2023")
table = pd.DataFrame({"median NDVI": [f"{v:.2f}" for v in m23.median(("y", "x")).values],
                      "clear obs/pixel": n_obs.sel(time="2023").median(("y", "x")).values.astype(int),
                      "cover %": (cover.sel(time="2023").values * 100).round().astype(int)},
                     index=m23.time.dt.strftime("%b").values)
print(table.T.to_string())
monthly.astype("float32").to_netcdf(OUT / "ndvi_monthly_median_2021_2025.nc")

# --- Step 7: per-field time series ----------------------------------------------------------------------
gbox = ds.odc.geobox
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"])
ts = ts.merge(blocks[["field_id", "crop", "area_ha"]], on="field_id")
ts.round(4).to_csv(OUT / "field_ndvi_timeseries.csv", index=False, date_format="%Y-%m-%d")
section("per-field time series")
print(f"{len(ts)} field-dates for {ts['field_id'].nunique()} fields; "
      f"median {ts.groupby('field_id').size().median():.0f} clear observations per field over five years")
print(ts.groupby(["crop", ts["date"].dt.year])["ndvi"].count().unstack().to_string())

# --- Step 8: Savitzky-Golay smoothing and phenology ---------------------------------------------------------
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)


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()])
pheno = pheno.merge(blocks[["field_id", "crop", "area_ha"]], on="field_id")
section("2023 phenology by crop (median of fields)")
summary = pheno.groupby("crop").agg(
    fields=("field_id", "size"),
    green_up=("green_up", lambda d: d.dropna().median().strftime("%d %b")),
    peak=("peak", lambda d: d.median().strftime("%d %b")),
    peak_ndvi=("peak_ndvi", "median"),
    senescence=("senescence", lambda d: d.dropna().median().strftime("%d %b")),
).sort_values("peak", key=lambda s: pd.to_datetime(s + " 2023", format="%d %b %Y"))
print(summary.round(2).to_string())
pheno.round(3).to_csv(OUT / "phenology_2023.csv", index=False, date_format="%Y-%m-%d")

# --- Figures ---------------------------------------------------------------------------------------------------
bb = gbox.boundingbox
ext = [bb.left, bb.right, bb.bottom, bb.top]


def km_axes(ax) -> None:
    """Label map axes in km from the lower-left corner of the study area."""
    xt, yt = np.arange(0, ext[1] - ext[0], 1000), np.arange(0, ext[3] - ext[2], 1000)
    ax.set_xticks(ext[0] + xt, [f"{v / 1000:.0f}" for v in xt])
    ax.set_yticks(ext[2] + yt, [f"{v / 1000:.0f}" for v in yt])
    ax.set_xlabel("km")


cmap = plt.get_cmap("Spectral").copy()
cmap.set_bad("#d9dcdb")

# 1. study area: true colour on a clear summer day, and the CROME crop blocks
summer = clear_frac[(clear_frac.index.month.isin([6, 7])) & (clear_frac.index.year == 2023)]
rgb_time = summer.idxmax()
rgb_ds = ds[["B02", "B03", "B04"]].sel(time=rgb_time).compute()
off = int(offset.sel(time=rgb_time))
rgb = np.dstack([(rgb_ds[b].values.astype("float32") - off) / 10000 for b in ("B04", "B03", "B02")])
rgb = np.clip(rgb / np.nanpercentile(rgb, 99), 0, 1) ** 0.8
section("figures")
print(f"true-colour date {pd.Timestamp(rgb_time):%d %b %Y}, {summer.max():.0%} clear")
fig, axes = plt.subplots(1, 2, figsize=(11, 5.6), sharey=True)
axes[0].imshow(rgb, extent=ext)
axes[0].set_title(f"Sentinel-2 true colour, {pd.Timestamp(rgb_time):%d %B %Y}")
axes[1].imshow(rgb.mean(axis=2), cmap="gray", extent=ext, vmin=0, vmax=1.4)
for crop, g in blocks.groupby("crop"):
    g.plot(ax=axes[1], color=COLOURS[crop], edgecolor="white", linewidth=0.3)
axes[1].set_title(f"CROME 2023 crop blocks of at least {MIN_FIELD_HA} ha")
axes[1].legend(handles=[Patch(color=COLOURS[c], label=f"{c} ({(blocks['crop'] == c).sum()})")
                        for c in COLOURS], loc="lower left", fontsize=8, framealpha=0.9)
for ax in axes:
    ax.set_xlim(ext[:2]); ax.set_ylim(ext[2:])
    km_axes(ax)
axes[0].set_ylabel("km")
axes[1].set_ylabel("")
fig.savefig(OUT / "01-study-area.png")
plt.close(fig)

# 2. how many clear observations: per date and per month
fig, axes = plt.subplots(2, 1, figsize=(10, 5.6), gridspec_kw={"height_ratios": [2, 1.2]})
cf = clear_frac
axes[0].scatter(cf.index, cf * 100, s=7, c=np.where(cf >= MIN_CLEAR, "#1d6fa3", "#b0b7b5"))
axes[0].axhline(MIN_CLEAR * 100, color="#c21f3a", lw=0.8, ls="--")
axes[0].set_ylabel("clear pixels (%)")
axes[0].set_title(f"Clear-sky share of the study area on each of {len(cf)} Sentinel-2 dates, 2021–2025 "
                  f"(blue = used, {int((cf >= MIN_CLEAR).sum())} dates)")
axes[0].xaxis.set_major_locator(mdates.YearLocator())
axes[0].xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
usable = cf[cf >= MIN_CLEAR]
counts = pd.Series(usable.index.month).value_counts().reindex(range(1, 13), fill_value=0) / 5
axes[1].bar(range(1, 13), counts.values, color="#2f9e44")
axes[1].set_xticks(range(1, 13), ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"])
axes[1].set_ylabel("usable dates\nper year (mean)")
for ax in axes:
    ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(OUT / "02-clear-observations.png")
plt.close(fig)

# 3. small multiples: monthly median NDVI in 2023
fig, axes = plt.subplots(3, 4, figsize=(11, 9.4))
for ax, t in zip(axes.ravel(), m23.time.values):
    im = ax.imshow(m23.sel(time=t), cmap=cmap, vmin=-0.1, vmax=0.9, extent=ext)
    nobs = int(n_obs.sel(time=t).median())
    ax.set_title(f"{pd.Timestamp(t):%B} ({nobs} clear obs.)", fontsize=8.5)
    ax.set_xticks([]); ax.set_yticks([])
fig.subplots_adjust(wspace=0.04, hspace=0.16, bottom=0.08)
cax = fig.add_axes([0.25, 0.04, 0.5, 0.015])
fig.colorbar(im, cax=cax, orientation="horizontal", label="NDVI, monthly median, 2023 (grey = no clear observation)")
fig.savefig(OUT / "03-monthly-ndvi-2023.png")
plt.close(fig)

# 4. one field per crop, raw and smoothed, over five years
def most_typical(crop: str) -> int:
    """The field whose 2023 curve correlates best with the median 2023 curve of its crop."""
    fields = pheno.loc[pheno["crop"] == crop, "field_id"]
    c23 = pd.DataFrame({f: curves[f]["2023"] for f in fields})
    return int(c23.corrwith(c23.median(axis=1)).idxmax())


area = blocks.set_index("field_id")["area_ha"]
example = {crop: most_typical(crop) for crop in ["Winter wheat", "Sugar beet", "Potato", "Maize"]}
print("most typical field per crop:", example)
fig, axes = plt.subplots(4, 1, figsize=(10, 8.6), sharex=True)
for ax, (crop, f) in zip(axes, example.items()):
    g = ts[ts["field_id"] == f].set_index("date")["ndvi"]
    ax.scatter(g.index, g.values, s=8, color="#8a9491", label="clear observations (field mean)")
    ax.plot(curves[f].index, curves[f].values, color=COLOURS[crop], lw=1.6, label="Savitzky–Golay upper envelope")
    for y in range(2021, 2026):
        ax.axvspan(pd.Timestamp(f"{y}-06-01"), pd.Timestamp(f"{y}-08-31"), color="#fff1b8", lw=0, zorder=0)
    ax.text(0.01, 0.95, f"Field {f}: {crop} in 2023 ({area[f]:.0f} ha, {len(g)} observations)",
            transform=ax.transAxes, fontweight="bold", color=COLOURS[crop], va="top")
    ax.set_ylim(-0.15, 1.05)
    ax.grid(alpha=0.3)
axes[0].legend(loc="lower right", fontsize=7.5, ncol=2, framealpha=0.9)
axes[2].set_ylabel("NDVI")
axes[-1].xaxis.set_major_locator(mdates.YearLocator())
axes[-1].xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
fig.suptitle("Five years of NDVI for four fields, labelled by their 2023 crop (yellow = June–August)", y=0.995)
fig.tight_layout()
fig.savefig(OUT / "04-field-time-series.png")
plt.close(fig)

# 5. phenology curves by crop in 2023
fig, ax = plt.subplots(figsize=(10, 5.2))
for crop in summary.index:
    ids = pheno.loc[pheno["crop"] == crop, "field_id"]
    stack = pd.concat([curves[f]["2022-10-01":"2023-12-31"] for f in ids], axis=1)
    q25, q50, q75 = stack.quantile(0.25, axis=1), stack.median(axis=1), stack.quantile(0.75, axis=1)
    ax.fill_between(q50.index, q25, q75, color=COLOURS[crop], alpha=0.15, lw=0)
    ax.plot(q50.index, q50, color=COLOURS[crop], lw=2, label=f"{crop} ({len(ids)} fields)")
    pk = pd.Timestamp(f"{summary.loc[crop, 'peak']} 2023")
    ax.plot(pk, q50.asof(pk), "o", color=COLOURS[crop], mec="white", ms=7)
ax.set_xlim(pd.Timestamp("2022-10-01"), pd.Timestamp("2023-12-31"))
ax.xaxis.set_major_locator(mdates.MonthLocator(bymonth=[1, 3, 5, 7, 9, 11]))
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b %Y"))
ax.set_ylabel("NDVI (smoothed)")
ax.set_title("2023 phenology curves: median and interquartile range of fields per crop; dots mark the median peak")
ax.legend(fontsize=8, loc="upper left", ncol=2)
ax.grid(alpha=0.3)
fig.savefig(OUT / "05-phenology-by-crop.png")
plt.close(fig)

# 6. crop rotation: monthly maximum NDVI per field, 2021-2025, grouped by 2023 crop
mon = ts.groupby(["field_id", pd.Grouper(key="date", freq="1MS")])["ndvi"].max().unstack()
mon = mon.reindex(columns=pd.date_range("2021-01-01", "2025-12-01", freq="MS"))
order = pheno[pheno["crop"].isin(example)].sort_values(["crop", "area_ha"], ascending=[True, False])
mat = mon.loc[order["field_id"]]
fig, ax = plt.subplots(figsize=(10, 6.4))
im = ax.imshow(np.ma.masked_invalid(mat.values), aspect="auto", cmap=cmap, vmin=-0.1, vmax=0.9,
               interpolation="nearest")
pos = 0
for crop, g in order.groupby("crop", sort=False):
    ax.axhline(pos - 0.5, color="black", lw=0.8)
    ax.text(-1, pos + len(g) / 2, crop, ha="right", va="center", color=COLOURS[crop], fontweight="bold")
    pos += len(g)
ax.set_yticks([])
ax.set_xticks([i * 12 - 0.5 for i in range(5)], [str(y) for y in range(2021, 2026)], ha="left")
ax.set_title("Monthly maximum NDVI per field, 2021–2025, grouped by the 2023 crop (grey = no clear observation)")
fig.colorbar(im, ax=ax, fraction=0.03, pad=0.01, label="NDVI")
fig.savefig(OUT / "06-crop-rotation-heatmap.png")
plt.close(fig)

# Animated GIF of the 2023 monthly composites (written to output/, not shown in the article)
frames = []
for t in m23.time.values:
    f, a = plt.subplots(figsize=(5, 5.4), dpi=100)
    a.imshow(m23.sel(time=t), cmap=cmap, vmin=-0.1, vmax=0.9, extent=ext)
    a.set_title(f"Methwold Fen NDVI, {pd.Timestamp(t):%B %Y}")
    a.set_xticks([]); a.set_yticks([])
    f.canvas.draw()
    frames.append(np.asarray(f.canvas.buffer_rgba())[..., :3].copy())
    plt.close(f)
iio.imwrite(OUT / "ndvi_2023.gif", np.stack(frames), duration=700, loop=0)
print(f"wrote {len(list(OUT.glob('*.png')))} figures, ndvi_2023.gif ({len(frames)} frames) and the CSV and NetCDF files")