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