Source code: Satellite timelapse in Python: a Sentinel-2 GIF, WebP and MP4 of Hinkley Point C

"""Satellite timelapse in Python: eleven summers of Hinkley Point C from Sentinel-2.

Hinkley Point C, the nuclear power station being built on the Somerset coast, went from a
freshly levelled earthworks platform in 2016 to reactor buildings, a marine jetty and a
finished southern compound by 2026.
This script builds one cloud-free April-September composite per year and turns them into
an animated GIF, an animated WebP and an MP4, with a date label and a scale bar on every frame.

Steps
  1. find every Sentinel-2 L2A scene of the area for April-September, 2016-2026 (STAC)
  2. measure the clear fraction of each scene inside the area from its SCL band; keep the clear ones
  3. load the kept scenes at 10 m and mask cloud and shadow with SCL (plus a small buffer)
  4. take a per-pixel median for each summer
  5. apply one stretch to every year so that frames are comparable
  6. draw labels and a scale bar and export GIF, WebP and MP4
  7. measure the change: area inside the construction site that differs from 2016, year by year

Run:  python timelapse.py      (a few minutes; it reads only small windows of each scene from the Planetary Computer)
"""
import io
import os
import sys
import time
from pathlib import Path

for var in ("PROJ_LIB", "PROJ_DATA", "GDAL_DATA"):   # stray installs (e.g. PostGIS) break the wheels' own copies
    os.environ.pop(var, None)

import imageio.v3 as iio
import imageio_ffmpeg
import matplotlib
import numpy as np
import odc.stac
import pandas as pd
import planetary_computer as pc
import PIL
import pystac_client
import xarray as xr
from odc.geo.geobox import GeoBox
from PIL import Image, ImageDraw, ImageFont
from scipy import ndimage

matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import font_manager
from matplotlib.patches import Rectangle

HERE = Path(__file__).resolve().parent
OUT = HERE / "output"
OUT.mkdir(exist_ok=True)
T_START = time.perf_counter()
plt.rcParams.update({"savefig.dpi": 150, "font.size": 10, "axes.spines.top": False, "axes.spines.right": False})

STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
CRS = "EPSG:32630"                                    # UTM 30N, the native grid of tile 30UVB
AOI = (487600, 5670700, 493600, 5675200)              # 6.0 x 4.5 km: Hinkley Point A, B and C, Stolford, the foreshore
SITE = (489200, 5671650, 490500, 5673150)             # 1.3 x 1.5 km box around Hinkley Point C, west of Hinkley Point B
TILE = "30UVB"
YEARS = list(range(2016, 2027))
SEASON = ("04-01", "09-30")                           # April to September: long days, low sun angle problems avoided
MIN_CLEAR = 0.4                                       # keep a scene if at least 40% of the area is clear
SCL_BAD = [0, 1, 3, 8, 9, 10]                         # no data, saturated, cloud shadow, cloud medium/high, cirrus
BANDS = ["B02", "B03", "B04", "B08"]

print("=== versions ===")
print(f"python {sys.version.split()[0]} | odc-stac {odc.stac.__version__} | xarray {xr.__version__} | numpy {np.__version__} | "
      f"pillow {PIL.__version__} | imageio-ffmpeg {imageio_ffmpeg.__version__} (ffmpeg {imageio_ffmpeg.get_ffmpeg_version()})")

# ------------------------------------------------------------------ Step 1: every candidate scene
print("\n=== Step 1: Sentinel-2 L2A scenes over Hinkley Point, April-September 2016-2026 ===")
catalog = pystac_client.Client.open(STAC_URL, modifier=pc.sign_inplace)
from pyproj import Transformer

to_ll = Transformer.from_crs(CRS, "EPSG:4326", always_xy=True)
lon0, lat0 = to_ll.transform(AOI[0], AOI[1])
lon1, lat1 = to_ll.transform(AOI[2], AOI[3])
bbox_ll = [lon0, lat0, lon1, lat1]
items = []
for y in YEARS:
    found = catalog.search(collections=["sentinel-2-l2a"], bbox=bbox_ll, datetime=f"{y}-{SEASON[0]}/{y}-{SEASON[1]}",
                           query={"s2:mgrs_tile": {"eq": TILE}, "eo:cloud_cover": {"lt": 80}}).item_collection()
    # keep the latest processing of each acquisition
    latest = {}
    for it in found:
        k = it.properties["datetime"][:10] + it.id.split("_")[3]          # date + relative orbit
        if k not in latest or it.properties["s2:generation_time"] > latest[k].properties["s2:generation_time"]:
            latest[k] = it
    items += list(latest.values())
    print(f"{y}: {len(latest):>3} scenes with tile cloud < 80%")
items = sorted(items, key=lambda i: i.datetime)
print(f"{len(items)} candidate scenes in total")

# ------------------------------------------------------------------ Step 2: clear fraction from SCL
print("\n=== Step 2: clear fraction inside the area, from each scene's SCL band ===")
gb20 = GeoBox.from_bbox(AOI, crs=CRS, resolution=20)
t0 = time.perf_counter()
scl = odc.stac.load(items, bands=["SCL"], geobox=gb20, resampling="nearest", pool=16)["SCL"]
print(f"read SCL for {scl.sizes['time']} scenes on a {gb20.shape[1]} x {gb20.shape[0]} grid in {time.perf_counter() - t0:.0f} s")
scl_np = scl.values
has_data = (scl_np != 0).mean(axis=(1, 2))
clear = (~np.isin(scl_np, SCL_BAD)).mean(axis=(1, 2))
scenes = pd.DataFrame({"id": [i.id for i in items], "date": pd.to_datetime([i.datetime.date() for i in items]),
                       "tile_cloud": [i.properties["eo:cloud_cover"] for i in items],
                       "baseline": [i.properties["s2:processing_baseline"] for i in items],
                       "coverage": has_data, "clear": clear})
scenes["keep"] = (scenes["coverage"] > 0.99) & (scenes["clear"] >= MIN_CLEAR)
scenes["year"] = scenes["date"].dt.year
summary = scenes.groupby("year").agg(candidates=("id", "size"), kept=("keep", "sum"),
                                     median_clear=("clear", "median"))
print(summary.to_string(formatters={"median_clear": "{:.0%}".format}))
print(f"kept {scenes['keep'].sum()} of {len(scenes)} scenes; correlation of tile cloud with area clear fraction: "
      f"{np.corrcoef(scenes['tile_cloud'], scenes['clear'])[0, 1]:.2f}")
scenes.drop(columns="id").to_csv(OUT / "scenes.csv", index=False)

fig, ax = plt.subplots(figsize=(13, 4.4))
kept = scenes[scenes["keep"]]
rej = scenes[~scenes["keep"]]
ax.scatter(rej["date"], rej["clear"] * 100, s=16, color="#bbbbbb", label=f"rejected ({len(rej)})", zorder=2)
ax.scatter(kept["date"], kept["clear"] * 100, s=24, color="#3288bd", edgecolor="black", lw=0.3,
           label=f"kept: ≥{MIN_CLEAR:.0%} of the area clear and full coverage ({len(kept)})", zorder=3)
ax.axhline(MIN_CLEAR * 100, color="#d53e4f", lw=1, ls="--")
for y in YEARS:
    n = int(summary.loc[y, "kept"])
    ax.text(pd.Timestamp(f"{y}-07-01"), 104, f"{n}", ha="center", fontsize=9, color="#3288bd", fontweight="bold")
    ax.axvspan(pd.Timestamp(f"{y}-{SEASON[0]}"), pd.Timestamp(f"{y}-{SEASON[1]}"), color="#fee08b", alpha=0.25, lw=0)
ax.set_ylim(0, 110)
ax.set_ylabel("clear pixels in the 6 × 4.5 km area (%)")
ax.set_title("Scene selection: every April–September Sentinel-2 scene, with the number kept per summer above each year")
ax.legend(loc="lower left", fontsize=8.5, ncol=2)
fig.tight_layout()
fig.savefig(OUT / "01-scene-selection.png", bbox_inches="tight", facecolor="white")
plt.close(fig)

# ------------------------------------------------------------------ Step 3: load and mask
print("\n=== Step 3: load the kept scenes at 10 m and mask cloud and shadow ===")
keep_items = [it for it, k in zip(items, scenes["keep"]) if k]
gb10 = GeoBox.from_bbox(AOI, crs=CRS, resolution=10)
t0 = time.perf_counter()
ds = odc.stac.load(keep_items, bands=BANDS + ["SCL"], geobox=gb10,
                   resampling={"SCL": "nearest", "*": "bilinear"}, pool=16)
print(f"loaded {ds.sizes['time']} scenes x {len(BANDS) + 1} bands on a {gb10.shape[1]} x {gb10.shape[0]} grid "
      f"in {time.perf_counter() - t0:.0f} s")

# Processing baseline 04.00 (25 January 2022) added an offset of 1000 to every DN. Without
# removing it, 2022+ frames would look brighter than 2016-2021 ones.
offset = np.array([1000 if float(b) >= 4.0 else 0 for b in scenes.loc[scenes["keep"], "baseline"]], dtype="float32")
print("baselines:", ", ".join(f"{b} ({n})" for b, n in scenes.loc[scenes['keep'], 'baseline'].value_counts().sort_index().items()))
refl = {}
for b in BANDS:
    dn = ds[b].values.astype("float32")
    refl[b] = np.where(dn == 0, np.nan, (dn - offset[:, None, None]) / 10000)
bad = np.isin(ds["SCL"].values, SCL_BAD)
bad_buf = np.stack([ndimage.binary_dilation(m, iterations=3) for m in bad])      # 30 m buffer catches hazy cloud edges
print(f"SCL masks {bad.mean() * 100:.1f}% of kept pixels; with a 30 m buffer {bad_buf.mean() * 100:.1f}%")
for b in BANDS:
    refl[b][bad_buf] = np.nan
dates = pd.to_datetime(ds["time"].values)
years = dates.year.values


def rgb(r, g, b, lo, hi, gamma=1.0):
    """Linear stretch between lo and hi (per band), then gamma; returns floats in 0..1."""
    x = np.dstack([r, g, b])
    x = np.clip((x - lo) / (hi - lo), 0, 1) ** (1 / gamma)
    return np.nan_to_num(x, nan=0.0)


SHOW_YEAR = 2021
idx_year = np.where(years == SHOW_YEAR)[0]
cloudiest = idx_year[np.argmax(bad[idx_year].mean(axis=(1, 2)))]
d = dates[cloudiest].date()
raw = [(ds[b].values[cloudiest].astype("float32") - offset[cloudiest]) / 10000 for b in ["B04", "B03", "B02"]]
lo_show, hi_show = np.array([0.01, 0.02, 0.01]), np.array([0.16, 0.14, 0.12])
fig, axes = plt.subplots(1, 3, figsize=(13, 3.9))
axes[0].imshow(rgb(*raw, lo_show, hi_show, 1.4))
axes[0].set_title(f"{d}: true colour", fontsize=10)
ov = rgb(*raw, lo_show, hi_show, 1.4) * 0.55
ov[bad[cloudiest]] = [0.99, 0.68, 0.38]
ov[bad_buf[cloudiest] & ~bad[cloudiest]] = [0.84, 0.24, 0.31]
axes[1].imshow(ov)
axes[1].set_title(f"SCL cloud and shadow, {bad[cloudiest].mean():.0%} (orange); 30 m buffer (red)", fontsize=10)
masked = rgb(*[np.where(bad_buf[cloudiest], np.nan, r) for r in raw], lo_show, hi_show, 1.4)
masked[bad_buf[cloudiest]] = 1.0
axes[2].imshow(masked)
axes[2].set_title("Masked: only clear pixels go into the median", fontsize=10)
for ax in axes:
    ax.set_xticks([])
    ax.set_yticks([])
fig.tight_layout()
fig.savefig(OUT / "02-cloud-mask.png", bbox_inches="tight", facecolor="white")
plt.close(fig)

# ------------------------------------------------------------------ Step 4: median per summer
print("\n=== Step 4: one median composite per summer ===")
comp, nclear = {}, {}
for y in YEARS:
    sel = years == y
    with np.errstate(all="ignore"):
        comp[y] = {b: np.nanmedian(refl[b][sel], axis=0) for b in BANDS}
    nclear[y] = np.isfinite(refl["B04"][sel]).sum(axis=0)
    gap = nclear[y] == 0
    # Pixels that were never clear (bright concrete and the sea wall, which SCL often calls cloud, or
    # persistent cloud) fall back to the median of the unmasked scenes so the frame has no holes.
    if gap.any():
        for b in BANDS:
            dn = ds[b].values[sel].astype("float32")
            unmasked = np.where(dn == 0, np.nan, (dn - offset[sel][:, None, None]) / 10000)
            with np.errstate(all="ignore"):
                comp[y][b][gap] = np.nanmedian(unmasked[:, gap], axis=0)
    print(f"{y}: {sel.sum():>2} scenes; clear observations per pixel min {nclear[y].min()}, median {np.median(nclear[y]):.0f}; "
          f"never clear {gap.mean() * 100:.2f}% (filled from unmasked scenes)")

lo_tmp, hi_tmp = np.array([0.01, 0.02, 0.01]), np.array([0.16, 0.14, 0.12])
fig, axes = plt.subplots(2, 4, figsize=(13, 5.6))
for ax, i in zip(axes.flat[:6], idx_year[np.linspace(0, len(idx_year) - 1, 6).round().astype(int)]):
    img = rgb(refl["B04"][i], refl["B03"][i], refl["B02"][i], lo_tmp, hi_tmp, 1.4)
    img[np.isnan(refl["B04"][i])] = 1.0
    ax.imshow(img)
    ax.set_title(f"{dates[i].date()} ({np.isnan(refl['B04'][i]).mean():.0%} masked)", fontsize=9)
axes.flat[6].imshow(rgb(comp[SHOW_YEAR]["B04"], comp[SHOW_YEAR]["B03"], comp[SHOW_YEAR]["B02"], lo_tmp, hi_tmp, 1.4))
axes.flat[6].set_title(f"{SHOW_YEAR} median of {len(idx_year)} masked scenes", fontsize=9, fontweight="bold")
imc = axes.flat[7].imshow(nclear[SHOW_YEAR], cmap="Spectral", vmin=0, vmax=len(idx_year))
axes.flat[7].set_title("clear observations per pixel", fontsize=9)
fig.colorbar(imc, ax=axes.flat[7], shrink=0.8)
for ax in axes.flat:
    ax.set_xticks([])
    ax.set_yticks([])
fig.suptitle(f"Six of the {len(idx_year)} masked scenes from summer {SHOW_YEAR} (white = masked) and their median", y=1.0)
fig.tight_layout()
fig.savefig(OUT / "03-median-composite.png", bbox_inches="tight", facecolor="white")
plt.close(fig)

# ------------------------------------------------------------------ Step 5: one stretch for all years
print("\n=== Step 5: a consistent stretch ===")
pool = np.stack([np.stack([comp[y][b] for b in ["B04", "B03", "B02"]], axis=-1) for y in YEARS])
LO = np.nanpercentile(pool, 1, axis=(0, 1, 2))
HI = np.nanpercentile(pool, 99, axis=(0, 1, 2))
GAMMA = 1.4
print("global stretch (reflectance) R/G/B low " + " / ".join(f"{v:.3f}" for v in LO)
      + ", high " + " / ".join(f"{v:.3f}" for v in HI) + f", gamma {GAMMA}")
lims = {}
frames_global, frames_auto = {}, {}
for y in YEARS:
    r, g, b = comp[y]["B04"], comp[y]["B03"], comp[y]["B02"]
    frames_global[y] = rgb(r, g, b, LO, HI, GAMMA)
    stack = np.dstack([r, g, b])
    lo_y, hi_y = np.nanpercentile(stack, 1, axis=(0, 1)), np.nanpercentile(stack, 99, axis=(0, 1))
    frames_auto[y] = rgb(r, g, b, lo_y, hi_y, GAMMA)
    lims[y] = (lo_y, hi_y)
hi_red = np.array([lims[y][1][0] for y in YEARS])
lo_red = np.array([lims[y][0][0] for y in YEARS])
print(f"per-frame red limits move between years: high {hi_red.min():.3f}-{hi_red.max():.3f}, "
      f"low {lo_red.min():.3f}-{lo_red.max():.3f} (the global stretch fixes them at {HI[0]:.3f} and {LO[0]:.3f})")
# The same surface reflectance therefore maps to a different screen brightness each year:
probe = 0.10
shown = [np.clip((probe - lo_red[i]) / (hi_red[i] - lo_red[i]), 0, 1) ** (1 / GAMMA) for i in range(len(YEARS))]
print(f"a red reflectance of {probe} is drawn at {min(shown):.2f}-{max(shown):.2f} of full brightness with per-frame "
      f"stretches, and at {np.clip((probe - LO[0]) / (HI[0] - LO[0]), 0, 1) ** (1 / GAMMA):.2f} in every frame with one stretch")

fig = plt.figure(figsize=(13, 4.6))
gs = fig.add_gridspec(2, 5, width_ratios=[1, 1, 1, 1, 1.5], wspace=0.06, hspace=0.18)
show = [2016, 2019, 2022, 2026]
for j, y in enumerate(show):
    for i, (frames, lab) in enumerate([(frames_auto, "per-frame stretch"), (frames_global, "one stretch")]):
        ax = fig.add_subplot(gs[i, j])
        ax.imshow(frames[y])
        ax.set_xticks([])
        ax.set_yticks([])
        if i == 0:
            ax.set_title(str(y), fontsize=10)
        if j == 0:
            ax.set_ylabel(lab, fontsize=10)
ax = fig.add_subplot(gs[:, 4])
ax.plot(YEARS, hi_red, "o-", color="#d53e4f", lw=2, label="per-frame 99th percentile")
ax.plot(YEARS, lo_red, "o-", color="#f46d43", lw=2, label="per-frame 1st percentile")
ax.axhline(HI[0], color="#3288bd", ls="--", lw=2, label="one stretch, all years")
ax.axhline(LO[0], color="#3288bd", ls="--", lw=2)
ax.set_ylabel("red-band stretch limit (reflectance)")
ax.yaxis.tick_right()
ax.yaxis.set_label_position("right")
ax.set_xticks(YEARS[::2])
ax.set_title("Per-frame limits move every year", fontsize=10)
ax.legend(fontsize=8, loc="center right")
fig.savefig(OUT / "04-consistent-stretch.png", bbox_inches="tight", facecolor="white")
plt.close(fig)

# ------------------------------------------------------------------ Step 6: frames and export
print("\n=== Step 6: label the frames and export GIF, WebP and MP4 ===")
SCALE = 4 / 3                                       # 600 x 450 px -> 800 x 600 px
W, H = int(gb10.shape[1] * SCALE), int(gb10.shape[0] * SCALE)
BAR = 40                                            # caption strip; 600 + 40 px is a multiple of 16 for H.264
font_path = font_manager.findfont("DejaVu Sans")
font_bold = font_manager.findfont(font_manager.FontProperties(family="DejaVu Sans", weight="bold"))
f_year = ImageFont.truetype(font_bold, 30)
f_small = ImageFont.truetype(font_path, 14)
f_scale = ImageFont.truetype(font_bold, 15)
n_kept = scenes[scenes["keep"]].groupby("year").size()


def make_frame(y):
    img = Image.fromarray((frames_global[y] * 255).round().astype("uint8")).resize((W, H), Image.Resampling.LANCZOS)
    canvas = Image.new("RGB", (W, H + BAR), (20, 22, 28))
    canvas.paste(img, (0, 0))
    d = ImageDraw.Draw(canvas)
    # date label: a dark box in the top-left corner
    d.rounded_rectangle((14, 14, 124, 60), radius=6, fill=(20, 22, 28))
    d.text((24, 18), str(y), font=f_year, fill=(255, 255, 255))
    # scale bar: 1 km = 100 px at 10 m, times the display scale
    px_km = 1000 / 10 * SCALE
    x0, y0 = W - 24 - px_km, H - 34
    d.rounded_rectangle((x0 - 12, y0 - 26, x0 + px_km + 12, y0 + 18), radius=6, fill=(20, 22, 28))
    d.rectangle((x0, y0, x0 + px_km / 2, y0 + 8), fill=(255, 255, 255))
    d.rectangle((x0 + px_km / 2, y0, x0 + px_km, y0 + 8), fill=(20, 22, 28), outline=(255, 255, 255), width=2)
    d.text((x0, y0 - 22), "0", font=f_scale, fill=(255, 255, 255))
    d.text((x0 + px_km - 30, y0 - 22), "1 km", font=f_scale, fill=(255, 255, 255))
    # north arrow
    d.polygon([(W - 34, 20), (W - 44, 46), (W - 24, 46)], fill=(255, 255, 255))
    d.text((W - 40, 48), "N", font=f_scale, fill=(255, 255, 255))
    # caption strip
    d.text((14, H + 3), f"Hinkley Point C, Somerset · Sentinel-2 median, April–September {y} ({n_kept[y]} clear scenes)",
           font=f_small, fill=(255, 255, 255))
    d.text((14, H + 21), "Contains modified Copernicus Sentinel data 2016–2026 · technicalgis.com",
           font=f_small, fill=(190, 190, 190))
    return np.asarray(canvas)


frames = [make_frame(y) for y in YEARS]
Image.fromarray(frames[-1]).save(OUT / "timelapse_poster.png", optimize=True)
DUR = [1100] * (len(frames) - 1) + [2200]           # ms per frame; hold the last year longer

# The GIF is 600 px wide: 256 colours per frame and LZW compression make GIFs several times larger than WebP or MP4.
gif_frames = [np.asarray(Image.fromarray(f).resize((600, 480), Image.Resampling.LANCZOS)) for f in frames]
iio.imwrite(OUT / "hinkley_timelapse.gif", gif_frames, duration=DUR, loop=0)
gif_mb = (OUT / "hinkley_timelapse.gif").stat().st_size / 1e6
pil_frames = [Image.fromarray(f) for f in frames]
pil_frames[0].save(OUT / "hinkley_timelapse.webp", save_all=True, append_images=pil_frames[1:], duration=DUR,
                   loop=0, quality=80, method=6)
webp_mb = (OUT / "hinkley_timelapse.webp").stat().st_size / 1e6

# MP4: 12 frames per second with a short cross-fade between years, H.264 in yuv420p so every browser plays it.
FPS, HOLD, FADE = 12, 12, 5
video = []
for i, f in enumerate(frames):
    video += [f] * (HOLD if i < len(frames) - 1 else HOLD * 2)
    if i < len(frames) - 1:
        a, b = f.astype("float32"), frames[i + 1].astype("float32")
        video += [((1 - t) * a + t * b).round().astype("uint8") for t in np.linspace(0, 1, FADE + 2)[1:-1]]
writer = imageio_ffmpeg.write_frames(str(OUT / "hinkley_timelapse.mp4"), (W, H + BAR), fps=FPS, codec="libx264",
                                     pix_fmt_out="yuv420p", quality=None,
                                     output_params=["-crf", "26", "-preset", "slow", "-movflags", "+faststart"])
writer.send(None)
for f in video:
    writer.send(np.ascontiguousarray(f))
writer.close()
mp4_mb = (OUT / "hinkley_timelapse.mp4").stat().st_size / 1e6
print(f"{len(frames)} frames of {W} x {H + BAR} px")
print(f"GIF  {gif_mb:5.2f} MB ({len(frames)} frames at 600 x 480 px, 256 colours per frame)")
print(f"WebP {webp_mb:5.2f} MB (lossy, quality 80)")
print(f"MP4  {mp4_mb:5.2f} MB ({len(video)} frames at {FPS} fps, {len(video) / FPS:.1f} s, H.264 CRF 26)")
print(f"poster PNG {(OUT / 'timelapse_poster.png').stat().st_size / 1e6:.2f} MB")

fig, axes = plt.subplots(3, 4, figsize=(13, 8.4))
for ax, y in zip(axes.flat, YEARS):
    ax.imshow(frames_global[y])
    ax.set_title(f"{y} ({n_kept[y]} scenes)", fontsize=10)
    ax.set_xticks([])
    ax.set_yticks([])
axes.flat[-1].axis("off")
axes.flat[-1].text(0.02, 0.5, "Hinkley Point C, Somerset\nApril–September medians\n6 × 4.5 km, 10 m pixels\n"
                   "one stretch for all years", fontsize=11, va="center")
fig.tight_layout()
fig.savefig(OUT / "05-frames-contact-sheet.png", dpi=110, bbox_inches="tight", facecolor="white")
plt.close(fig)

# ------------------------------------------------------------------ Step 7: measure the change
print("\n=== Step 7: how much of the site changes from year to year ===")
c0, r0 = int((SITE[0] - AOI[0]) / 10), int((AOI[3] - SITE[3]) / 10)
c1, r1 = int((SITE[2] - AOI[0]) / 10), int((AOI[3] - SITE[1]) / 10)
site_ha = (r1 - r0) * (c1 - c0) * 0.01
RGB = ["B04", "B03", "B02"]
ndvi = {y: (comp[y]["B08"] - comp[y]["B04"]) / (comp[y]["B08"] + comp[y]["B04"]) for y in YEARS}
diff = {y: np.sqrt(sum((comp[y][b] - comp[YEARS[0]][b]) ** 2 for b in RGB)) for y in YEARS}
CHANGE = 0.05                                          # visible-band reflectance distance that counts as changed
rows = []
for y in YEARS:
    rows.append({"year": y,
                 "changed_since_2016_ha": float((diff[y][r0:r1, c0:c1] > CHANGE).sum() * 0.01),
                 "unvegetated_ha": float((ndvi[y][r0:r1, c0:c1] < 0.3).sum() * 0.01),
                 "site_mean_ndvi": float(np.nanmean(ndvi[y][r0:r1, c0:c1]))})
change = pd.DataFrame(rows)
print(change.to_string(index=False, formatters={"changed_since_2016_ha": "{:.0f}".format,
                                                 "unvegetated_ha": "{:.0f}".format, "site_mean_ndvi": "{:.2f}".format}))
print(f"site box {site_ha:.0f} ha; changed since 2016: {change.changed_since_2016_ha.iloc[-1]:.0f} ha by {YEARS[-1]} "
      f"({change.changed_since_2016_ha.iloc[-1] / site_ha:.0%}); unvegetated share {change.unvegetated_ha.min() / site_ha:.0%}-"
      f"{change.unvegetated_ha.max() / site_ha:.0%} in every year")
change.to_csv(OUT / "site_change.csv", index=False)

fig, axes = plt.subplots(1, 4, figsize=(14, 4.2), gridspec_kw={"width_ratios": [1, 1, 1.2, 1.2]})
pad = 30
crop = (slice(max(r0 - pad, 0), r1 + pad), slice(max(c0 - pad, 0), c1 + pad))
for ax, y in zip(axes[:2], [YEARS[0], YEARS[-1]]):
    ax.imshow(frames_global[y][crop])
    ax.set_title(f"Site, summer {y}", fontsize=10)
im = axes[2].imshow(diff[YEARS[-1]][crop], cmap="Spectral_r", vmin=0, vmax=0.2)
axes[2].set_title(f"Change {YEARS[0]} to {YEARS[-1]}", fontsize=10)
fig.colorbar(im, ax=axes[2], shrink=0.8, label="reflectance distance")
for ax in axes[:3]:
    ax.add_patch(Rectangle((c0 - crop[1].start, r0 - crop[0].start), c1 - c0, r1 - r0, fill=False, ec="black", lw=1.2))
    ax.set_xticks([])
    ax.set_yticks([])
ax = axes[3]
ax.bar(change["year"], change["changed_since_2016_ha"], color="#f46d43", edgecolor="black", lw=0.4,
       label=f"changed since {YEARS[0]} (> {CHANGE})")
ax.plot(change["year"], change["unvegetated_ha"], "o-", color="#3288bd", lw=2, label="NDVI < 0.3")
ax.axhline(site_ha, color="#5e4fa2", ls="--", lw=1)
ax.text(YEARS[0] - 0.4, site_ha * 1.03, f"site box {site_ha:.0f} ha", color="#5e4fa2", fontsize=8.5)
ax.set_ylabel("hectares in the site box")
ax.set_ylim(0, site_ha * 1.4)
ax.set_xticks(YEARS[::2])
ax.legend(fontsize=8, loc="upper left", bbox_to_anchor=(0, 0.93))
ax.set_title("Change versus bare ground", fontsize=10)
fig.tight_layout()
fig.savefig(OUT / "06-site-change.png", bbox_inches="tight", facecolor="white")
plt.close(fig)

print(f"\nfinished in {(time.perf_counter() - T_START) / 60:.1f} min")