Source code: Mann-Kendall trend test on rasters in Python: Theil-Sen NDVI trends with xarray

"""Pixel-wise NDVI trend analysis in Python: Mann-Kendall and Theil-Sen with xarray.apply_ufunc on
25 years of MODIS Terra NDVI (MOD13A1 v6.1) over the Loess Plateau around Yan'an, China, with the
seasonal Mann-Kendall test, serial-correlation corrections and Benjamini-Hochberg FDR control.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/raster-trend-analysis-mann-kendall-python/
Data: MODIS/Terra Vegetation Indices 16-Day L3 Global 500 m (MOD13A1 v061), NASA LP DAAC, served by the
Microsoft Planetary Computer. MODIS data are freely reusable without restriction.
"""
import sys
import time
import warnings
from concurrent.futures import ThreadPoolExecutor
from importlib.metadata import version
from pathlib import Path

import numpy as np
import pandas as pd
import xarray as xr
import rioxarray  # noqa: F401  (registers the .rio accessor)
import rasterio
from rasterio.warp import transform_bounds
from rasterio.windows import from_bounds
import pystac_client
import planetary_computer
import pymannkendall as pmk
from scipy.stats import norm
from statsmodels.stats.multitest import multipletests
import matplotlib
matplotlib.use("Agg")   # render to files only; no GUI backend
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm, TwoSlopeNorm
from matplotlib.patches import Patch

warnings.filterwarnings("ignore")
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9,
                     "figure.facecolor": "white", "axes.facecolor": "white"})

DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)
SEED = 42
rng = np.random.default_rng(SEED)
BBOX = (108.5, 36.0, 110.5, 38.5)        # lon/lat: Yan'an (south) to Yulin (north), Shaanxi, China
YEARS = range(2001, 2025)                # 24 complete years (2025 is incomplete in the catalogue)
SEASON = [5, 6, 7, 8, 9]                 # growing season, May-September
ALPHA = 0.05
CITIES = {"Yan'an": (109.49, 36.59), "Yulin": (109.73, 38.29), "Suide": (110.25, 37.50)}


def section(title):
    print(f"\n=== {title} ===")


section("versions")
print(f"python {sys.version.split()[0]} | numpy {version('numpy')} | xarray {version('xarray')} | "
      f"rasterio {version('rasterio')} | pymannkendall {version('pymannkendall')} | scipy {version('scipy')}")

# ---------------------------------------------------------------- Step 1: find and read the MODIS composites
section("MOD13A1 v6.1 items on the Planetary Computer")
cube_file = DATA / "mod13a1_ndvi_cube.npz"
if not cube_file.exists():
    catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",
                                        modifier=planetary_computer.sign_inplace)
    items = catalog.search(collections=["modis-13A1-061"], bbox=BBOX,
                           datetime=f"{YEARS[0]}-01-01/{YEARS[-1] + 1}-12-31").item_collection()
    items = sorted([it for it in items if it.id.startswith("MOD13A1")], key=lambda it: it.properties["start_datetime"])
    start = [pd.Timestamp(it.properties["start_datetime"]).tz_localize(None) for it in items]
    tiles = sorted({it.id.split(".")[2] for it in items})
    print(f"{len(items)} Terra composites, tiles {tiles}, {start[0]:%Y-%m-%d} .. {start[-1]:%Y-%m-%d}")

    with rasterio.open(items[0].assets["500m_16_days_NDVI"].href) as src:
        crs = src.crs
        win = from_bounds(*transform_bounds("EPSG:4326", crs, *BBOX, densify_pts=51), src.transform)
        win = win.round_offsets().round_lengths()
        transform = rasterio.windows.transform(win, src.transform)

    def read(item):
        out = []
        for key in ("500m_16_days_NDVI", "500m_16_days_pixel_reliability"):
            for attempt in range(5):
                try:
                    with rasterio.open(item.assets[key].href) as src:
                        out.append(src.read(1, window=win))
                    break
                except rasterio.errors.RasterioIOError:
                    time.sleep(5 * (attempt + 1))
                    planetary_computer.sign_inplace(item)
            else:
                raise RuntimeError(f"could not read {item.id} {key}")
        return out

    t0 = time.time()
    with rasterio.Env(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", GDAL_HTTP_MAX_RETRY=5, GDAL_HTTP_RETRY_DELAY=2):
        with ThreadPoolExecutor(16) as pool:
            arrays = list(pool.map(read, items))
    print(f"read {len(items)} x 2 windows of {win.height} x {win.width} pixels in {time.time() - t0:.0f} s")
    np.savez_compressed(cube_file, ndvi=np.stack([a[0] for a in arrays]), reliability=np.stack([a[1] for a in arrays]),
                        time=np.array(start, dtype="datetime64[s]"),
                        x=transform.c + transform.a * (np.arange(win.width) + 0.5),
                        y=transform.f + transform.e * (np.arange(win.height) + 0.5), crs=crs.to_wkt())
z = np.load(cube_file)
cube = xr.Dataset({"ndvi": (("time", "y", "x"), z["ndvi"]), "reliability": (("time", "y", "x"), z["reliability"])},
                  coords={"time": pd.DatetimeIndex(z["time"]), "y": z["y"], "x": z["x"]})
crs = rasterio.crs.CRS.from_wkt(str(z["crs"]))
dup = cube.indexes["time"].duplicated()
print(f"cube: {cube.sizes['time']} composites x {cube.sizes['y']} x {cube.sizes['x']} pixels (500 m, MODIS sinusoidal); "
      f"{dup.sum()} duplicate composites (reprocessed copies) dropped")
cube = cube.isel(time=~dup)
expected = pd.DatetimeIndex([pd.Timestamp(f"{y}-01-01") + pd.Timedelta(days=16 * k) for y in range(2001, 2026) for k in range(23)])
missing = expected.difference(cube.indexes["time"])
print(f"composites missing from the catalogue: {', '.join(missing.strftime('%Y-%m-%d'))}")
print(f"-> 2025 lacks July-November, so the analysis uses {YEARS[0]}-{YEARS[-1]} ({len(YEARS)} growing seasons)")
cube = cube.sel(time=cube.time.dt.year.isin(list(YEARS)))

# ---------------------------------------------------------------- Step 2: QA masking
section("QA masking with the pixel reliability layer")
rel = cube.reliability
for k, name in [(-1, "fill / no data"), (0, "good"), (1, "marginal"), (2, "snow or ice"), (3, "cloudy")]:
    print(f"  reliability {k:>2} ({name}): {float((rel == k).mean()):.1%} of pixel-composites in the window")
ndvi = cube.ndvi.where((rel >= 0) & (rel <= 1) & (cube.ndvi > -2000)) * 1e-4
ndvi = ndvi.assign_coords(time=cube.time + pd.Timedelta(days=8))       # date each composite by its mid-point

# ---------------------------------------------------------------- Step 3: aggregate to remove seasonality
section("growing-season aggregation")
gs = ndvi.where(ndvi.time.dt.month.isin(SEASON), drop=True)
n_valid = gs.notnull().groupby("time.year").sum()
annual = gs.groupby("time.year").mean().where(n_valid >= 5)              # need 5 of ~10 composites
monthly = gs.resample(time="MS").mean()
monthly = monthly.where(monthly.time.dt.month.isin(SEASON), drop=True)
monthly = monthly.assign_coords(year=("time", monthly.time.dt.year.data), month=("time", monthly.time.dt.month.data))
monthly = monthly.set_index(time=["year", "month"]).unstack("time").transpose("year", "month", "y", "x")
mean_ndvi = annual.mean("year")
years_ok = annual.notnull().sum("year")
# the sinusoidal window is a parallelogram in lon/lat: keep only pixels whose centre is inside BBOX
lon, lat = rasterio.warp.transform(crs, "EPSG:4326", *[a.ravel() for a in np.meshgrid(cube.x, cube.y)])
lon = np.reshape(lon, mean_ndvi.shape); lat = np.reshape(lat, mean_ndvi.shape)
inside = xr.DataArray((lon >= BBOX[0]) & (lon <= BBOX[2]) & (lat >= BBOX[1]) & (lat <= BBOX[3]), dims=("y", "x"))
mean_ndvi = mean_ndvi.where(inside)
land = inside & (mean_ndvi > 0.1) & (years_ok == len(YEARS))
print(f"growing-season composites per year: {gs.time.size / len(YEARS):.1f} on average; median valid (QA-passed) "
      f"per pixel-year inside the box {float(n_valid.where(inside).median()):.0f}")
print(f"pixels analysed: {int(land.sum()):,} of {int(inside.sum()):,} inside the box (dropped: mean growing-season "
      f"NDVI <= 0.1 or any year with < 5 valid composites)")
annual = annual.where(land)
monthly = monthly.where(land)


# ---------------------------------------------------------------- Step 4: vectorised Mann-Kendall + Theil-Sen
def mk_theil_sen(x, t):
    """Mann-Kendall test and Theil-Sen slope along the last axis of x (time), for any leading shape.
    Missing values are skipped pairwise. Variance ignores ties (fine for continuous NDVI means)."""
    i, j = np.triu_indices(x.shape[-1], 1)
    d = x[..., j] - x[..., i]
    s = np.nansum(np.sign(d), axis=-1)
    n = np.sum(~np.isnan(x), axis=-1)
    var = n * (n - 1) * (2 * n + 5) / 18.0
    z = np.where(s > 0, (s - 1) / np.sqrt(var), np.where(s < 0, (s + 1) / np.sqrt(var), 0.0))
    p = 2 * norm.sf(np.abs(z))
    slope = np.nanmedian(d / (t[j] - t[i]), axis=-1)
    return np.stack([s, var, z, p, slope], axis=-1)


def seasonal_mk(x, t):
    """Seasonal Mann-Kendall (Hirsch et al. 1982): x has dims (..., year, season). S and Var(S) are summed
    over seasons; the slope is the median of all within-season pairwise slopes."""
    i, j = np.triu_indices(x.shape[-2], 1)
    s = var = 0
    slopes = []
    for m in range(x.shape[-1]):
        xm = x[..., m]
        d = xm[..., j] - xm[..., i]
        s = s + np.nansum(np.sign(d), axis=-1)
        n = np.sum(~np.isnan(xm), axis=-1)
        var = var + n * (n - 1) * (2 * n + 5) / 18.0
        slopes.append(d / (t[j] - t[i]))
    z = np.where(s > 0, (s - 1) / np.sqrt(var), np.where(s < 0, (s + 1) / np.sqrt(var), 0.0))
    p = 2 * norm.sf(np.abs(z))
    slope = np.nanmedian(np.concatenate(slopes, axis=-1), axis=-1)
    return np.stack([s, var, z, p, slope], axis=-1)


section("vectorised Mann-Kendall and Theil-Sen (annual growing-season means)")
t_years = np.asarray(annual.year, dtype=float)
t0 = time.time()
stats = xr.apply_ufunc(mk_theil_sen, annual.chunk({"y": 64, "x": 64}), kwargs={"t": t_years},
                       input_core_dims=[["year"]], output_core_dims=[["stat"]], dask="parallelized",
                       output_dtypes=[float], dask_gufunc_kwargs={"output_sizes": {"stat": 5}})
stats = stats.assign_coords(stat=["S", "varS", "z", "p", "slope"]).compute()
print(f"{int(land.sum()):,} pixels in {time.time() - t0:.1f} s")
p_mk = stats.sel(stat="p").where(land)
slope = stats.sel(stat="slope").where(land) * 10                        # NDVI per decade

# cross-check against pymannkendall on random pixels
yy, xx = np.nonzero(land.values)
pick = rng.choice(len(yy), 6, replace=False)
print("cross-check with pymannkendall.original_test (pixel: S ours/pmk, p ours/pmk, slope ours/pmk):")
for k in pick:
    series = annual.values[:, yy[k], xx[k]]
    ref = pmk.original_test(series)
    ours = stats.values[yy[k], xx[k]]
    print(f"  ({yy[k]:>3},{xx[k]:>3}): S {ours[0]:+5.0f}/{ref.s:+5.0f}  p {ours[3]:.4f}/{ref.p:.4f}  "
          f"slope {ours[4]:+.5f}/{ref.slope:+.5f}")

section("seasonal Mann-Kendall (monthly means, May-September as seasons)")
sstats = xr.apply_ufunc(seasonal_mk, monthly.chunk({"y": 64, "x": 64}), kwargs={"t": t_years},
                        input_core_dims=[["year", "month"]], output_core_dims=[["stat"]], dask="parallelized",
                        output_dtypes=[float], dask_gufunc_kwargs={"output_sizes": {"stat": 5}})
sstats = sstats.assign_coords(stat=["S", "varS", "z", "p", "slope"]).compute()
print("cross-check with pymannkendall.seasonal_test(period=5):")
for k in pick[:3]:
    mser = monthly.values[:, :, yy[k], xx[k]]           # (year, month)
    ref = pmk.seasonal_test(mser.reshape(-1), period=len(SEASON))
    ours = sstats.values[yy[k], xx[k]]
    print(f"  ({yy[k]:>3},{xx[k]:>3}): S {ours[0]:+5.0f}/{ref.s:+5.0f}  p {ours[3]:.4f}/{ref.p:.4f}  "
          f"slope {ours[4]:+.5f}/{ref.slope:+.5f}")
p_smk = sstats.sel(stat="p").where(land)
both = land.values
agree = np.mean(np.sign(stats.sel(stat="z").values[both]) == np.sign(sstats.sel(stat="z").values[both]))
print(f"annual MK and seasonal MK agree on the sign of the trend for {agree:.1%} of pixels; "
      f"p < 0.05 (unadjusted): annual {float((p_mk < ALPHA).sum() / land.sum()):.1%}, "
      f"seasonal {float((p_smk < ALPHA).sum() / land.sum()):.1%}")

# ---------------------------------------------------------------- Step 5: serial correlation
section("serial correlation")


def lag1(x):
    """Lag-1 autocorrelation along the last axis after removing a least-squares linear trend."""
    t = np.arange(x.shape[-1], dtype=float)
    xd = x - np.nanmean(x, -1, keepdims=True)
    tt = t - t.mean()
    b = np.nansum(xd * tt, -1, keepdims=True) / np.sum(tt ** 2)
    e = xd - b * tt
    return np.nansum(e[..., 1:] * e[..., :-1], -1) / np.nansum(e ** 2, -1)


r1_annual = lag1(np.moveaxis(annual.values, 0, -1))[both]
anom = ndvi.groupby("time.dayofyear") - ndvi.groupby("time.dayofyear").mean()   # 16-day deseasonalised anomalies
anom_v = np.moveaxis(anom.fillna(0).values, 0, -1)[both]
r1_16day = lag1(anom_v)
crit = 1.96 / np.sqrt(len(YEARS))
print(f"lag-1 autocorrelation, annual growing-season means: median {np.median(r1_annual):+.2f}; "
      f"|r1| > {crit:.2f} (approx. 5% test) for {np.mean(np.abs(r1_annual) > crit):.1%} of pixels")
print(f"lag-1 autocorrelation, 16-day deseasonalised anomalies: median {np.median(r1_16day):+.2f}")

sub = rng.choice(len(yy), 2000, replace=False)
tests = {"Original MK": pmk.original_test, "Hamed-Rao modified MK": pmk.hamed_rao_modification_test,
         "Yue-Wang modified MK": pmk.yue_wang_modification_test, "Pre-whitening (PW)": pmk.pre_whitening_modification_test,
         "Trend-free pre-whitening (TFPW)": pmk.trend_free_pre_whitening_modification_test}
sig_share = {}
for name, f in tests.items():
    ps = np.array([f(annual.values[:, yy[k], xx[k]]).p for k in sub])
    sig_share[name] = np.mean(ps < ALPHA)
    print(f"  {name:<32} p < 0.05 for {sig_share[name]:.1%} of 2,000 random pixels")
naive16 = np.array([pmk.original_test(anom_v[k]).p for k in sub[:500]])
hr16 = np.array([pmk.hamed_rao_modification_test(anom_v[k]).p for k in sub[:500]])
print(f"  for contrast, MK on the 16-day anomaly series ({anom_v.shape[1]} values) of 500 pixels: p < 0.05 for "
      f"{np.mean(naive16 < ALPHA):.1%} (original) vs {np.mean(hr16 < ALPHA):.1%} (Hamed-Rao)")

# ---------------------------------------------------------------- Step 6: multiple testing
section("false discovery rate across pixels")
pv = p_mk.values[both]
m = len(pv)
order = np.argsort(pv)
ranked = pv[order] * m / np.arange(1, m + 1)
q = np.minimum.accumulate(ranked[::-1])[::-1]
q_bh = np.empty(m); q_bh[order] = np.minimum(q, 1)
q_sm = multipletests(pv, alpha=ALPHA, method="fdr_bh")[1]
print(f"own Benjamini-Hochberg vs statsmodels: max |difference| in q = {np.abs(q_bh - q_sm).max():.2e}")
p_cut = pv[q_bh <= ALPHA].max()
print(f"tests: {m:,}; p < 0.05 unadjusted: {np.sum(pv < ALPHA):,} ({np.mean(pv < ALPHA):.1%}); "
      f"expected false positives if no trend anywhere: {ALPHA * m:,.0f}")
print(f"BH q < 0.05: {np.sum(q_bh < ALPHA):,} ({np.mean(q_bh < ALPHA):.1%}); equivalent p-value threshold {p_cut:.4f}")
print(f"Bonferroni p < {ALPHA / m:.1e}: {np.sum(pv < ALPHA / m):,} ({np.mean(pv < ALPHA / m):.1%})")
qmap = xr.full_like(p_mk, np.nan); qmap.values[both] = q_bh
sl = slope.values[both]
green = (q_bh < ALPHA) & (sl > 0); brown = (q_bh < ALPHA) & (sl < 0)
print(f"significant greening: {green.sum():,} pixels ({green.mean():.1%}), median slope "
      f"{np.median(sl[green]):+.3f} NDVI/decade")
print(f"significant browning: {brown.sum():,} pixels ({brown.mean():.1%}), median slope "
      f"{np.median(sl[brown]):+.3f} NDVI/decade")
print(f"area-mean Theil-Sen slope: {np.mean(sl):+.4f} NDVI/decade; mean growing-season NDVI "
      f"{float(annual.sel(year=YEARS[0]).mean()):.3f} ({YEARS[0]}) -> {float(annual.sel(year=YEARS[-1]).mean()):.3f} ({YEARS[-1]})")

# significance is not magnitude
small_sig = np.sum((q_bh < ALPHA) & (np.abs(sl) < 0.02))
big_ns = np.sum((q_bh >= ALPHA) & (np.abs(sl) > 0.05))
print(f"significant but small (|slope| < 0.02/decade): {small_sig:,}; large (> 0.05/decade) but not significant: {big_ns:,}")

# sensor sensitivity: Terra's orbit has drifted since 2020
section("sensitivity to the end of the series (Terra orbit drift after 2020)")
st19 = mk_theil_sen(np.moveaxis(annual.sel(year=slice(2001, 2019)).values, 0, -1)[both], t_years[:19])
sig19 = st19[:, 3] < ALPHA
print(f"2001-2019: p < 0.05 for {sig19.mean():.1%} of pixels (full series {np.mean(pv < ALPHA):.1%}); "
      f"slope correlation with {YEARS[0]}-{YEARS[-1]} "
      f"{np.corrcoef(st19[:, 4] * 10, sl)[0, 1]:.2f}; same sign {np.mean(np.sign(st19[:, 4]) == np.sign(sl)):.1%}")

# ---------------------------------------------------------------- Figures
section("figures")


def to_geo(da, resampling=rasterio.enums.Resampling.nearest):
    da = da.rio.write_crs(crs).rio.set_spatial_dims(x_dim="x", y_dim="y")
    return da.rio.reproject("EPSG:4326", resampling=resampling).rio.clip_box(*BBOX)


def decorate(ax):
    for name, (lo, la) in CITIES.items():
        ax.plot(lo, la, "o", ms=4, mfc="white", mec="black", mew=0.8)
        ax.text(lo + 0.04, la + 0.03, name, fontsize=7.5,
                bbox=dict(boxstyle="round,pad=0.15", fc="white", ec="none", alpha=0.8))
    ax.set_xlim(BBOX[0], BBOX[2]); ax.set_ylim(BBOX[1], BBOX[3])
    ax.set_xticks(np.arange(BBOX[0], BBOX[2] + 0.01, 0.5))
    ax.set_aspect(1 / np.cos(np.radians(37.25)))
    ax.set_xlabel("Longitude (°E)"); ax.set_ylabel("Latitude (°N)")


def show(ax, da, **kw):
    g = to_geo(da)
    return ax.imshow(g.values, extent=(float(g.x.min()), float(g.x.max()), float(g.y.min()), float(g.y.max())),
                     interpolation="nearest", **kw)


# Figure 1: mean NDVI and valid observations
fig, axes = plt.subplots(1, 2, figsize=(10, 6.2), constrained_layout=True)
im = show(axes[0], mean_ndvi.where(mean_ndvi > -1), cmap="YlGn", vmin=0, vmax=0.8)
fig.colorbar(im, ax=axes[0], shrink=0.75, label=f"Mean growing-season NDVI, {YEARS[0]}–{YEARS[-1]}")
axes[0].set_title("Mean May–September NDVI (MOD13A1, 500 m)")
good_share = (rel.sel(time=rel.time.dt.month.isin(SEASON)) == 0).mean("time").where(inside)
im = show(axes[1], 100 * good_share, cmap="cividis", vmin=50, vmax=100)
fig.colorbar(im, ax=axes[1], shrink=0.75, extend="min",
             label="May–September composites rated 'good' (reliability 0), %")
axes[1].set_title("QA: share of good-quality composites")
for ax in axes:
    decorate(ax)
fig.savefig(OUT / "01-mean-ndvi-and-qa.png"); plt.close(fig)

# Figure 2: Theil-Sen slope map
fig, ax = plt.subplots(figsize=(6.4, 6.6), constrained_layout=True)
lim = 0.15
im = show(ax, slope, cmap="BrBG", norm=TwoSlopeNorm(0, -lim, lim))
cb = fig.colorbar(im, ax=ax, shrink=0.8, extend="both")
cb.set_label("Theil-Sen slope (NDVI units per decade)")
ax.set_title(f"Growing-season NDVI trend, {YEARS[0]}–{YEARS[-1]}")
decorate(ax)
fig.savefig(OUT / "02-theil-sen-slope.png"); plt.close(fig)

# Figure 3: significance after FDR control vs family-wise (Bonferroni) control
cat_cmap = ListedColormap(["#a6611a", "#dfc27d", "#e8e8e8", "#80cdc1", "#018571"])
fig, axes = plt.subplots(1, 2, figsize=(10.5, 6.2), constrained_layout=True)
for ax, pmap, cut, title in [(axes[0], qmap, ALPHA, "Benjamini-Hochberg FDR: q < 0.05"),
                             (axes[1], p_mk, ALPHA / m, "Bonferroni: p < 0.05 / m")]:
    sig_px = pmap < cut
    cat = xr.where(sig_px, xr.where(slope > 0, 4, 0), xr.where(slope > 0, 3, xr.where(slope < 0, 1, 2))).where(land)
    show(ax, cat, cmap=cat_cmap, norm=BoundaryNorm(np.arange(6) - 0.5, 5))
    ax.set_title(f"{title} ({float(sig_px.sum() / land.sum()):.1%} of pixels)")
    decorate(ax)
fig.legend(handles=[Patch(color=c, label=l) for c, l in zip(
    cat_cmap.colors[::-1], ["Significant greening", "Greening, not significant", "No change (slope 0)",
                            "Browning, not significant", "Significant browning"])],
    loc="outside lower center", ncols=5, fontsize=8)
fig.savefig(OUT / "03-significant-trends-fdr.png"); plt.close(fig)

# Figure 4: time series at example pixels
sl_full = slope.values; q_full = qmap.values
cands = {
    "Strong significant greening": np.nanargmax(np.where(q_full < 0.001, sl_full, np.nan)),
    "Significant browning": np.nanargmin(np.where(q_full < ALPHA, sl_full, np.nan)),
    "Large slope, not significant": np.nanargmax(np.where(q_full > 0.2, np.abs(sl_full), np.nan)),
    "Small slope, significant": np.nanargmin(np.where(q_full < 0.01, np.abs(sl_full), np.nan)),
}
fig, axes = plt.subplots(2, 2, figsize=(11, 6.4), sharex=True, constrained_layout=True)
lon_lat = {}
for ax, (title, flat) in zip(axes.flat, cands.items()):
    iy, ix = np.unravel_index(flat, sl_full.shape)
    lo, la = rasterio.warp.transform(crs, "EPSG:4326", [float(cube.x[ix])], [float(cube.y[iy])])
    lon_lat[title] = (lo[0], la[0])
    raw = ndvi[:, iy, ix]
    ax.plot(raw.time, raw, ".", color="#bbbbbb", ms=2, label="16-day composites (QA-passed)")
    a = annual[:, iy, ix]
    ax.plot(pd.to_datetime([f"{y}-07-15" for y in a.year.values]), a, "o-", color="#0072B2", ms=3.5,
            lw=1, label="May–September mean")
    st = stats.values[iy, ix]
    tt = t_years
    icpt = np.median(a.values - st[4] * tt)
    ax.plot(pd.to_datetime([f"{y}-07-15" for y in tt.astype(int)]), icpt + st[4] * tt, "--", color="#D55E00",
            lw=1.5, label="Theil-Sen line")
    ax.set_title(f"{title}\n{la[0]:.2f}°N {lo[0]:.2f}°E: slope {st[4] * 10:+.3f}/decade, "
                 f"p = {st[3]:.2g}, q = {q_full[iy, ix]:.2g}", fontsize=8.5)
    ax.set_ylabel("NDVI"); ax.grid(color="0.92")
h, l = axes[0, 0].get_legend_handles_labels()
fig.legend(h, l, loc="outside lower center", ncols=3, fontsize=8)
fig.savefig(OUT / "04-example-time-series.png"); plt.close(fig)
for k, v in lon_lat.items():
    print(f"example pixel '{k}': {v[1]:.3f} N, {v[0]:.3f} E")

# Figure 5: serial correlation and the modified tests
fig, axes = plt.subplots(1, 2, figsize=(11, 4.2), constrained_layout=True)
bins = np.linspace(-0.8, 1, 46)
axes[0].hist(r1_16day, bins=bins, color="#E69F00", alpha=0.8, label="16-day deseasonalised anomalies")
axes[0].hist(r1_annual, bins=bins, color="#0072B2", alpha=0.8, label="Annual growing-season means")
axes[0].axvspan(-crit, crit, color="0.85", zorder=0, label=f"|r$_1$| < {crit:.2f} (n = {len(YEARS)})")
axes[0].set_xlabel("Lag-1 autocorrelation of detrended series"); axes[0].set_ylabel("Pixels")
axes[0].legend(fontsize=7.5); axes[0].set_title("Aggregation removes most serial correlation")
names = list(sig_share)
axes[1].barh(range(len(names)), [100 * sig_share[n] for n in names], color="#56B4E9")
for i, n in enumerate(names):
    axes[1].text(100 * sig_share[n] + 0.5, i, f"{100 * sig_share[n]:.1f}%", va="center", fontsize=8)
axes[1].set_yticks(range(len(names)), names); axes[1].invert_yaxis()
axes[1].set_xlabel("Pixels with p < 0.05 (%, 2,000 random pixels)")
axes[1].set_xlim(0, max(sig_share.values()) * 115)
axes[1].set_title("Mann-Kendall variants on the annual series")
fig.savefig(OUT / "05-serial-correlation.png"); plt.close(fig)

# Figure 6: significance vs magnitude
fig, ax = plt.subplots(figsize=(7.5, 4.6), constrained_layout=True)
lp = -np.log10(np.maximum(q_bh, 1e-12))
ymax = np.ceil(lp.max()) + 0.5
hb = ax.hexbin(sl, lp, gridsize=70, bins="log", cmap="viridis", mincnt=1, extent=(-0.2, 0.25, 0, ymax))
ax.axhline(-np.log10(ALPHA), color="#D55E00", lw=1)
ax.text(-0.195, -np.log10(ALPHA) + 0.15, "q = 0.05", color="#D55E00", fontsize=8)
for xv in (-0.02, 0.02):
    ax.axvline(xv, color="0.4", ls=":", lw=0.8)
ax.set_xlim(-0.2, 0.25); ax.set_ylim(0, ymax)
ax.set_xlabel("Theil-Sen slope (NDVI units per decade)"); ax.set_ylabel("−log10(BH-adjusted q)")
fig.colorbar(hb, ax=ax, label="Pixels per hexagon")
ax.set_title("Significance is not magnitude: q depends on noise as well as slope")
fig.savefig(OUT / "06-significance-vs-magnitude.png"); plt.close(fig)
print("wrote 6 figures to output/")