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