"""Shoreline change on the Holderness coast (East Yorkshire) from Sentinel-2, 2017 to 2026.

Searches the Microsoft Planetary Computer STAC API for every Sentinel-2 L2A scene of tile 30UXE,
computes MNDWI, separates land from water with Otsu's threshold, extracts sub-pixel shorelines with
skimage.measure.find_contours, casts shore-normal transects every 100 m and measures where each
shoreline crosses each transect. Tide levels at each overpass come from a harmonic tide model fitted
to open Environment Agency tide-gauge data at Bridlington and Spurn Point. The beach slope is
estimated from the data, the shorelines are corrected to mean sea level, and the end-point rate (EPR)
and linear regression rate (LRR) are computed for each transect and compared with the East Riding of
Yorkshire Council's cliff-erosion monitoring.
Run with:  python coastal_erosion_shorelines.py
"""
import platform
import re
import time
from concurrent.futures import ThreadPoolExecutor
from importlib.metadata import version
from pathlib import Path

import geopandas as gpd
import matplotlib

matplotlib.use("Agg")
import matplotlib.dates as mdates
import matplotlib.patheffects as pe
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import planetary_computer
import pystac_client
import rasterio
import rasterio.features
import rasterio.transform
import requests
import skimage
import utide
from matplotlib.colors import TwoSlopeNorm
from pypdf import PdfReader
from pyproj import Transformer
from rasterio.enums import Resampling
from rasterio.windows import Window
from scipy import ndimage, stats
from shapely.geometry import LineString, MultiLineString, Point, box
from skimage.filters import threshold_otsu
from skimage.measure import find_contours

HERE = Path(__file__).parent
OUT, DATA, CACHE = HERE / "output", HERE / "data", HERE / "cache"
for d in (OUT, DATA, CACHE, CACHE / "shorelines"):
    d.mkdir(exist_ok=True)
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
pd.set_option("display.width", 140)
pd.set_option("display.max_columns", 20)

print("=== versions ===")
print(f"python {platform.python_version()} | pystac-client {pystac_client.__version__} | "
      f"planetary-computer {version('planetary-computer')} | rasterio {rasterio.__version__} | "
      f"scikit-image {skimage.__version__} | utide {version('utide')} | numpy {np.__version__} | pandas {pd.__version__}")


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


def get(url: str, **kw) -> requests.Response:
    for attempt in range(4):
        try:
            r = requests.get(url, timeout=120, **kw)
            if r.status_code < 500:
                return r
        except requests.RequestException as exc:
            print(f"retrying {url[:70]}: {exc.__class__.__name__}")
        time.sleep(10 * (attempt + 1))
    raise RuntimeError(f"failed: {url}")


CRS = "EPSG:32630"          # UTM 30N, the native grid of Sentinel-2 tile 30UXE
TILE = "30UXE"
PERIOD = "2017-01-01/2026-10-18"
to_utm = Transformer.from_crs("EPSG:27700", CRS, always_xy=True)
to_ll = Transformer.from_crs(CRS, "EPSG:4326", always_xy=True)

# --- Step 1: the council's erosion monitoring profiles -------------------------------------
# East Riding of Yorkshire Council publishes cliff-top erosion rates for its monitoring profiles twice a year.
ERYC_URL = "https://www.eastriding.gov.uk/download/east-riding-coastal-erosion-statistics/"
pdf_path = DATA / "eryc_cliff_erosion_table.pdf"
if not pdf_path.exists():
    r = get(ERYC_URL, headers={"User-Agent": "Mozilla/5.0 (TechnicalGIS tutorial script)"})   # blocks default agents
    r.raise_for_status()
    pdf_path.write_bytes(r.content)
text = "\n".join(p.extract_text() for p in PdfReader(pdf_path).pages)
NUM = r"(\d{6,7}\.\d)"
rows = []
for line in text.splitlines():
    m = re.match(rf"^(\d+)\s+(.+?)\s+{NUM}\s+{NUM}\s+{NUM}\s+{NUM}\s+(.*)$", line.strip())
    if not m:
        continue
    vals = re.findall(r"-?\d+\.\d+", m.group(7).replace("*", ""))
    rows.append({"profile": int(m.group(1)), "location": m.group(2),
                 "e1": float(m.group(3)), "n1": float(m.group(4)), "e2": float(m.group(5)), "n2": float(m.group(6)),
                 "defended": "Defended" in m.group(7),
                 # columns: loss Mar-Sep 2025, loss Sep 2025-Mar 2026, rate 1852-1989, rate 1989-2026, height, max loss
                 "rate_1989_2026": float(vals[3]) if len(vals) >= 6 else np.nan,
                 "loss_last_year": float(vals[0]) + float(vals[1]) if len(vals) >= 6 else np.nan})
eryc = pd.DataFrame(rows)
eryc["length_m"] = np.hypot(eryc.e2 - eryc.e1, eryc.n2 - eryc.n1)
bad = eryc[eryc.length_m > 3000]
eryc = eryc[eryc.length_m <= 3000].copy()
section("East Riding erosion profiles (PDF table)")
print(f"{len(rows)} profiles parsed; dropped {len(bad)} with a typing error in the coordinates: "
      f"{bad.profile.tolist()}")
STUDY = (22, 89)            # North Caravan Park, Skipsea to the north end of the Withernsea defences
eryc = eryc[eryc.profile.between(*STUDY)].reset_index(drop=True)
for a, b in (("e1", "n1"), ("e2", "n2")):
    eryc[f"{a}_utm"], eryc[f"{b}_utm"] = to_utm.transform(eryc[a].values, eryc[b].values)
undef = eryc[~eryc.defended & eryc.rate_1989_2026.notna()]
print(f"study stretch: profiles {STUDY[0]}-{STUDY[1]}, {len(eryc)} profiles, {len(undef)} with a 1989-2026 rate")
print(f"council cliff-top erosion rate 1989-2026: median {undef.rate_1989_2026.median():.2f} m/yr, "
      f"range {undef.rate_1989_2026.min():.2f}-{undef.rate_1989_2026.max():.2f} m/yr")
print(eryc[["profile", "location", "rate_1989_2026", "loss_last_year"]].head(8).to_string(index=False))

# --- Step 2: tide levels from open Environment Agency tide gauges ---------------------------
# The real-time API serves the last ~28 days; the daily archive dumps go back about a year.
GAUGES = {"Bridlington": ("L3001-level-tidal_level-i-15_min-m", 54.0798),
          "Spurn Point": ("L3370-level-tidal_level-i-15_min-mAOD", 53.5779)}
EA = "https://environment.data.gov.uk/flood-monitoring"
tide_csv = DATA / "ea_tide_gauges.csv"
if not tide_csv.exists():
    frames = []
    for name, (measure, _) in GAUGES.items():
        r = get(f"{EA}/id/measures/{measure}/readings", params={"since": "2026-09-01T00:00:00Z", "_sorted": "",
                                                                 "_limit": 10000}).json()
        frames.append(pd.DataFrame([{"gauge": name, "time": i["dateTime"], "level": i["value"]} for i in r["items"]]))
    # One day a month from the archive (each file is ~60 MB with every EA sensor, so we stream and filter).
    for day in pd.date_range("2025-11-01", "2026-08-31", freq="MS").strftime("%Y-%m-%d").tolist() + ["2026-06-15"]:
        r = requests.get(f"{EA}/archive/readings-{day}.csv", stream=True, timeout=300)
        if r.status_code != 200:
            print(f"archive {day}: HTTP {r.status_code}, skipped")
            continue
        keep = []
        for raw in r.iter_lines():
            line = raw.decode()
            for name, (measure, _) in GAUGES.items():
                if measure in line:
                    t, _, v = line.split(",")[:3]
                    keep.append({"gauge": name, "time": t, "level": v})
        frames.append(pd.DataFrame(keep))
        print(f"archive {day}: {len(keep)} readings")
    pd.concat(frames).to_csv(tide_csv, index=False)
tides = pd.read_csv(tide_csv)
tides["time"] = pd.to_datetime(tides["time"], utc=True, format="ISO8601")
tides["level"] = pd.to_numeric(tides["level"], errors="coerce")
tides = tides.dropna().drop_duplicates(["gauge", "time"]).sort_values(["gauge", "time"])
# Bridlington harbour dries out near low water: the gauge then sits flat at its floor (-1.61 m).
# Those readings are not the tide, so drop everything within 3 cm of each gauge's minimum.
floor = tides.groupby("gauge")["level"].transform("min")
dry = (tides.gauge == "Bridlington") & (tides.level <= floor + 0.03)
tides = tides[~dry]
tides["day"] = tides["time"].dt.strftime("%Y-%m-%d")
section("tide gauge readings (Environment Agency, 15-minute)")
print(tides.groupby("gauge").agg(readings=("level", "size"), first=("time", "min"), last=("time", "max"),
                                 days=("day", "nunique")).to_string())
print(f"dropped {dry.sum()} Bridlington readings sitting on the dry-harbour floor ({floor[dry].iloc[0]:.2f} m)")

# Fit a harmonic tide model per gauge. Hold out three archive days to test the predictions.
# A year of scattered days resolves the main constituents; asking for more would overfit.
CONSTITUENTS = ["M2", "S2", "N2", "K2", "K1", "O1", "P1", "Q1", "M4", "MS4", "MN4", "M6", "2N2", "MU2", "NU2", "L2"]
HOLDOUT = ["2026-01-01", "2026-04-01", "2026-07-01"]
models, val = {}, []
for name, (_, lat) in GAUGES.items():
    g = tides[tides.gauge == name]
    fit = g[~g.day.isin(HOLDOUT)]
    coef = utide.solve(fit.time.dt.tz_localize(None).values, fit.level.values - fit.level.mean(), lat=lat,
                       constit=CONSTITUENTS, method="ols", conf_int="none", trend=False, verbose=False)
    coef["offset"] = fit.level.mean()
    models[name] = coef
    test = g[g.day.isin(HOLDOUT)]
    pred = utide.reconstruct(test.time.dt.tz_localize(None).values, coef, verbose=False).h + coef["offset"]
    val.append(pd.DataFrame({"gauge": name, "time": test.time.values, "day": test.day.values,
                             "observed": test.level.values, "predicted": pred}))
val = pd.concat(val)
section("harmonic tide model (utide)")
for name, coef in models.items():
    top = pd.Series(coef["A"], index=coef["name"]).sort_values(ascending=False).head(5)
    print(f"{name:<12} {len(coef['name'])} constituents; largest amplitudes (m): "
          + ", ".join(f"{k} {v:.2f}" for k, v in top.items()))
val["error"] = val.predicted - val.observed
print("hold-out days, predicted minus observed (m):")
print(val.groupby(["gauge", "day"]).error.agg(rmse=lambda e: np.sqrt((e ** 2).mean()), bias="mean").round(2)
      .unstack(0).to_string())


# --- Step 3: search the Sentinel-2 archive ---------------------------------------------------
STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
lines_utm = [LineString([(r.e1_utm, r.n1_utm), (r.e2_utm, r.n2_utm)]) for r in eryc.itertuples()]
corridor = gpd.GeoSeries(lines_utm, crs=CRS).buffer(600).union_all()
bbox_ll = gpd.GeoSeries([corridor], crs=CRS).to_crs(4326).total_bounds
catalog = pystac_client.Client.open(STAC_URL)
items = list(catalog.search(collections=["sentinel-2-l2a"], bbox=bbox_ll, datetime=PERIOD,
                            query={"s2:mgrs_tile": {"eq": TILE}, "eo:cloud_cover": {"lt": 80}}).items())
best = {}
for it in items:   # a few dates exist twice (original and reprocessed): keep the newest processing baseline
    key = (it.datetime.date(), it.properties["sat:relative_orbit"])
    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("Sentinel-2 search")
print(f"{len(items)} scenes of tile {TILE} with tile cloud cover < 80%, {items[0].datetime:%Y-%m-%d} to "
      f"{items[-1].datetime:%Y-%m-%d}")
print(pd.Series([i.properties["sat:relative_orbit"] for i in items]).value_counts().rename("scenes by orbit")
      .to_string())

# --- Step 4: read only the 512 x 512 pixel blocks that the coast passes through ---------------
GRID_X0, GRID_Y0, RES, BLOCK = 600000.0, 6000000.0, 10.0, 512     # tile 30UXE, 10 m bands
blocks = [(r, c) for r in range(22) for c in range(22)
          if corridor.intersects(box(GRID_X0 + c * BLOCK * RES, GRID_Y0 - (r + 1) * BLOCK * RES,
                                     GRID_X0 + (c + 1) * BLOCK * RES, GRID_Y0 - r * BLOCK * RES))]
r0, r1 = min(b[0] for b in blocks), max(b[0] for b in blocks) + 1
c0, c1 = min(b[1] for b in blocks), max(b[1] for b in blocks) + 1
H, W = (r1 - r0) * BLOCK, (c1 - c0) * BLOCK
X0, Y0 = GRID_X0 + c0 * BLOCK * RES, GRID_Y0 - r0 * BLOCK * RES          # top-left of the mosaic


def rc_to_xy(rows, cols):
    """Mosaic pixel coordinates (fractional, pixel centres at integers) to UTM metres."""
    return X0 + (np.asarray(cols) + 0.5) * RES, Y0 - (np.asarray(rows) + 0.5) * RES


in_corridor = np.zeros((H, W), bool)
for r, c in blocks:
    in_corridor[(r - r0) * BLOCK:(r - r0 + 1) * BLOCK, (c - c0) * BLOCK:(c - c0 + 1) * BLOCK] = True
section("analysis mosaic")
print(f"{len(blocks)} blocks of 512 x 512 pixels at 10 m; mosaic {W} x {H} pixels "
      f"({W * RES / 1000:.1f} x {H * RES / 1000:.1f} km), {CRS}")

GDAL_ENV = dict(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", GDAL_HTTP_MAX_RETRY="5", GDAL_HTTP_RETRY_DELAY="3", VSI_CACHE="TRUE")


def read_band(item, band: str, scale: int, resampling=Resampling.nearest) -> np.ndarray:
    """Read one band into the 10 m mosaic. 20 m bands (scale=2) are resampled to 10 m."""
    out = np.zeros((H, W), np.uint16 if band != "SCL" else np.uint8)
    with rasterio.Env(**GDAL_ENV), rasterio.open(item.assets[band].href) as src:
        for r, c in blocks:
            win = Window(c * BLOCK // scale, r * BLOCK // scale, BLOCK // scale, BLOCK // scale)
            out[(r - r0) * BLOCK:(r - r0 + 1) * BLOCK, (c - c0) * BLOCK:(c - c0 + 1) * BLOCK] = src.read(
                1, window=win, out_shape=(BLOCK, BLOCK), resampling=resampling, boundless=True, fill_value=0)
    return out


def reflectance(item, dn: np.ndarray) -> np.ndarray:
    """Surface reflectance; baseline 04.00+ (from January 2022) adds 1000 to every digital number."""
    offset = 1000 if float(item.properties["s2:processing_baseline"]) >= 4.0 else 0
    return np.where(dn > 0, (dn.astype("float32") - offset) / 10000, np.nan)


BAD_SCL = [0, 1, 3, 8, 9, 10]           # no data, defective, cloud shadow, cloud medium/high, cirrus


def scene_mndwi(item, zone=None, min_clear=0.0):
    """MNDWI = (green - SWIR1) / (green + SWIR1) and a validity mask (clear, observed pixels).
    Reads the small SCL band first and gives up early if too little of `zone` is clear."""
    item = planetary_computer.sign(item)
    scl = read_band(item, "SCL", 2)
    bad = ndimage.binary_dilation(np.isin(scl, BAD_SCL) & in_corridor, iterations=3) | (scl == 0)
    if zone is not None and (~bad[zone]).mean() < min_clear:
        return None, None
    green = reflectance(item, read_band(item, "B03", 1))
    swir = reflectance(item, read_band(item, "B11", 2, Resampling.bilinear))
    with np.errstate(divide="ignore", invalid="ignore"):
        mndwi = (green - swir) / (green + swir)
    valid = ~bad & np.isfinite(mndwi) & in_corridor
    return mndwi, valid


# --- Step 5: a reference shoreline from a median composite ----------------------------------
ref_path = DATA / "reference_mndwi.npy"
clear = sorted(items, key=lambda i: i.properties["eo:cloud_cover"])[:20]
if not ref_path.exists():
    with ThreadPoolExecutor(4) as ex:
        stack = [np.where(v, m, np.nan).astype("float16") for m, v in ex.map(scene_mndwi, clear)]
    np.save(ref_path, np.nanmedian(np.stack(stack).astype("float32"), axis=0))
ref_mndwi = np.load(ref_path)
ref_thr = threshold_otsu(ref_mndwi[np.isfinite(ref_mndwi)])
contours = find_contours(np.nan_to_num(ref_mndwi, nan=-1), ref_thr, mask=np.isfinite(ref_mndwi))
longest = max(contours, key=len)
ref_line = LineString(np.column_stack(rc_to_xy(longest[:, 0], longest[:, 1])))
# Keep the stretch between the first and last council profile, then smooth it so transects are shore-normal.
ends = [ref_line.project(LineString(lines_utm[i]).intersection(ref_line)) if LineString(lines_utm[i])
        .intersects(ref_line) else ref_line.project(Point(lines_utm[i].coords[1])) for i in (0, -1)]
d = np.arange(min(ends), max(ends), 10.0)
pts = np.array([ref_line.interpolate(x).coords[0] for x in d])
smooth = np.column_stack([ndimage.uniform_filter1d(pts[:, k], 51, mode="nearest") for k in (0, 1)])
ref_smooth = LineString(smooth)
if ref_smooth.coords[0][1] < ref_smooth.coords[-1][1]:       # run the chainage from north (Skipsea) to south
    ref_smooth = ref_smooth.reverse()
section("reference shoreline")
print(f"median MNDWI of the {len(clear)} least cloudy scenes; Otsu threshold {ref_thr:.3f}")
print(f"reference shoreline {ref_smooth.length / 1000:.1f} km long, Skipsea to Withernsea")

# --- Step 6: shore-normal transects every 100 m ----------------------------------------------
SPACING, LAND, SEA = 100, 250, 350
transects = []
for s in np.arange(50, ref_smooth.length - 50, SPACING):
    p0, p1 = np.array(ref_smooth.interpolate(s - 25).coords[0]), np.array(ref_smooth.interpolate(s + 25).coords[0])
    t = (p1 - p0) / np.linalg.norm(p1 - p0)
    n = np.array([t[1], -t[0]])                       # right-hand normal; the sea is to the east
    if n[0] < 0:
        n = -n
    c = np.array(ref_smooth.interpolate(s).coords[0])
    transects.append({"chainage_km": s / 1000, "geometry": LineString([c - LAND * n, c + SEA * n]),
                      "x0": c[0] - LAND * n[0], "y0": c[1] - LAND * n[1], "nx": n[0], "ny": n[1]})
tr = gpd.GeoDataFrame(transects, crs=CRS)
tr["lon"], tr["lat"] = to_ll.transform(tr.geometry.centroid.x.values, tr.geometry.centroid.y.values)
print(f"{len(tr)} transects, {SPACING} m apart, {LAND} m landward and {SEA} m seaward of the reference line")

# --- Step 7: extract a sub-pixel shoreline from every scene ----------------------------------
ZONE_M, MIN_CLEAR = 400, 0.2
zone = rasterio.features.rasterize([ref_smooth.buffer(ZONE_M)], out_shape=(H, W),
                                   transform=rasterio.transform.from_origin(X0, Y0, RES, RES)).astype(bool)
# Pixels every 5 m along each transect, to test whether the transect is cloud-free in a scene.
samples = []
for g in tr.geometry:
    xs, ys = np.array([g.interpolate(s).coords[0] for s in np.arange(0, g.length, 5)]).T
    samples.append(((Y0 - ys) / RES).astype(int) * W + ((xs - X0) / RES).astype(int))


def extract_shoreline(item):
    """Otsu threshold on MNDWI in the shore zone, then the sub-pixel 'marching squares' contour."""
    path = CACHE / "shorelines" / f"{item.id}.npz"
    if path.exists():
        return dict(np.load(path))
    mndwi, valid = scene_mndwi(item, zone, MIN_CLEAR)
    if mndwi is None:
        res = {"status": np.array("cloudy")}
    else:
        ok = valid & zone
        thr = threshold_otsu(mndwi[ok])
        cs = find_contours(np.where(ok, mndwi, -1).astype("float32"), thr, mask=ok)
        cs = [np.column_stack(rc_to_xy(c[:, 0], c[:, 1])) for c in cs if len(c) >= 10]
        coords = np.concatenate([np.vstack([c, [np.nan, np.nan]]) for c in cs]) if cs else np.empty((0, 2))
        clear_tr = np.array([ok.ravel()[s].mean() for s in samples])
        res = {"status": np.array("ok"), "threshold": np.array(thr), "zone_clear": np.array(ok[zone].mean()),
               "contours": coords, "transect_clear": clear_tr}
    np.savez_compressed(path, **res)
    return res


t_start = time.time()
with ThreadPoolExecutor(6) as ex:
    results = list(ex.map(extract_shoreline, items))
section("shoreline extraction")
status = pd.Series([str(r["status"]) for r in results]).value_counts()
print(f"{len(items)} scenes processed in {time.time() - t_start:.0f} s: {status.get('ok', 0)} with a shoreline, "
      f"{status.get('cloudy', 0)} with less than {MIN_CLEAR:.0%} of the shore zone clear")
thr_all = np.array([float(r["threshold"]) for r in results if str(r["status"]) == "ok"])
print(f"Otsu MNDWI thresholds: median {np.median(thr_all):.3f}, 5th-95th percentile "
      f"{np.percentile(thr_all, 5):.3f} to {np.percentile(thr_all, 95):.3f}")

# --- Step 8: intersect every shoreline with every transect -----------------------------------
LAND = 250.0                # transect origin is 250 m landward of the reference line
MIN_TR_CLEAR, MAX_SPREAD = 0.9, 50.0
tr_geoms = tr.geometry.values
records = []
for item, r in zip(items, results):
    if str(r["status"]) != "ok" or len(r["contours"]) == 0:
        continue
    c = r["contours"]
    parts = np.split(c, np.where(np.isnan(c[:, 0]))[0])
    lines = MultiLineString([p[~np.isnan(p[:, 0])] for p in parts if (~np.isnan(p[:, 0])).sum() >= 2])
    for j, g in enumerate(tr_geoms):
        if r["transect_clear"][j] < MIN_TR_CLEAR:
            continue
        hit = g.intersection(lines)
        if hit.is_empty:
            continue
        d = np.array([g.project(p) for p in getattr(hit, "geoms", [hit])]) - LAND
        if d.max() - d.min() > MAX_SPREAD:      # several crossings far apart: whitewater, ords, cloud edge
            continue
        records.append((item.id, item.datetime, j, float(np.median(d))))
obs = pd.DataFrame(records, columns=["scene", "time", "transect", "x"])
obs["time"] = pd.to_datetime(obs.time, utc=True)
obs["year"] = (obs.time - pd.Timestamp("2017-01-01", tz="UTC")).dt.total_seconds() / (365.25 * 86400) + 2017

# Predicted tide at the moment of each overpass, interpolated between the gauges by latitude.
scene_t = obs.drop_duplicates("scene")[["scene", "time"]].reset_index(drop=True)
t_naive = scene_t.time.dt.tz_localize(None).values
scene_t["h_brid"] = utide.reconstruct(t_naive, models["Bridlington"], verbose=False).h
scene_t["h_spurn"] = utide.reconstruct(t_naive, models["Spurn Point"], verbose=False).h
obs = obs.merge(scene_t.drop(columns="time"), on="scene")
w = (tr.lat.values[obs.transect] - GAUGES["Spurn Point"][1]) / (GAUGES["Bridlington"][1] - GAUGES["Spurn Point"][1])
obs["tide"] = w * obs.h_brid + (1 - w) * obs.h_spurn
obs = obs.drop(columns=["h_brid", "h_spurn"])
section("shoreline positions on transects")
print(f"{len(obs):,} shoreline positions from {obs.scene.nunique()} scenes "
      f"({obs.time.min():%Y-%m-%d} to {obs.time.max():%Y-%m-%d})")
print(f"positions per transect: median {obs.groupby('transect').size().median():.0f}, "
      f"min {obs.groupby('transect').size().min()}")
print(f"overpass time {scene_t.time.dt.strftime('%H:%M').min()} to {scene_t.time.dt.strftime('%H:%M').max()} UTC; "
      f"predicted tide at overpass {obs.tide.min():+.2f} to {obs.tide.max():+.2f} m (median {obs.tide.median():+.2f} m)")

# --- Step 9: tide correction and change rates per transect -----------------------------------
def robust_fit(X, y):
    """Least squares, then refit without points more than 3 robust standard deviations away."""
    keep = np.ones(len(y), bool)
    for _ in range(2):
        beta, *_ = np.linalg.lstsq(X[keep], y[keep], rcond=None)
        res = y - X @ beta
        mad = 1.4826 * np.median(np.abs(res[keep] - np.median(res[keep])))
        keep = np.abs(res) < 3 * mad
    return beta, keep


def ols_rate(t, x):
    """Linear regression rate (LRR) in m/yr with its 95 % confidence half-width."""
    lr = stats.linregress(t, x)
    return lr.slope, stats.t.ppf(0.975, len(t) - 2) * lr.stderr


rows = []
obs["keep"] = False
for j, g in obs.groupby("transect"):
    X = np.column_stack([np.ones(len(g)), g.year - 2017, g.tide])
    beta, keep = robust_fit(X, g.x.values)
    obs.loc[g.index, "keep"] = keep
    rows.append({"transect": j, "k_raw": beta[2], "n": keep.sum()})
fit = pd.DataFrame(rows).set_index("transect")
# The tide coefficient k (metres of shoreline shift per metre of tide) is -1 / beach slope. It is noisy
# transect by transect, so smooth it alongshore with an 11-transect (1.1 km) running median.
fit["k"] = fit.k_raw.rolling(11, center=True, min_periods=1).median().clip(upper=-1)
fit["beach_slope"] = -1 / fit.k
obs = obs.join(fit.k, on="transect")
obs["x_corr"] = obs.x - obs.k * obs.tide          # shoreline position at mean sea level

out = []
t0, t1 = obs.time.min(), obs.time.max()
for j, g in obs[obs.keep].groupby("transect"):
    lrr, ci = ols_rate(g.year, g.x_corr)
    lrr_raw, ci_raw = ols_rate(g.year, g.x)
    first = g[g.time < t0 + pd.Timedelta(days=365)]
    last = g[g.time > t1 - pd.Timedelta(days=365)]
    epr = (last.x_corr.median() - first.x_corr.median()) / ((last.year.median() - first.year.median()))
    detr = lambda v: v - np.polyval(np.polyfit(g.year, v, 1), g.year)
    out.append({"transect": j, "lrr": lrr, "lrr_ci95": ci, "lrr_no_tide": lrr_raw, "lrr_no_tide_ci95": ci_raw,
                "epr": epr, "rmse_no_tide": detr(g.x).std(), "rmse_tide": detr(g.x_corr).std()})
rates = tr[["chainage_km", "lon", "lat"]].join(fit[["n", "k", "beach_slope"]]).join(pd.DataFrame(out).set_index("transect"))
section("tide correction")
print(f"tide coefficient k (m of shoreline shift per m of tide): median {rates.k.median():.1f}, "
      f"equivalent beach slope 1:{-rates.k.median():.0f} (5th-95th pct 1:{-rates.k.quantile(.95):.0f} to 1:{-rates.k.quantile(.05):.0f})")
print(f"scatter about the trend: median {rates.rmse_no_tide.median():.1f} m without tide correction, "
      f"{rates.rmse_tide.median():.1f} m with it")
print(f"outliers removed: {(~obs.keep).sum():,} of {len(obs):,} positions ({(~obs.keep).mean():.1%})")
section("shoreline change rates (negative = landward retreat)")
print(f"LRR, tide-corrected:  median {rates.lrr.median():+.2f} m/yr, IQR {rates.lrr.quantile(.25):+.2f} to {rates.lrr.quantile(.75):+.2f}")
print(f"LRR, no tide correct: median {rates.lrr_no_tide.median():+.2f} m/yr")
print(f"EPR (first vs last 12 months): median {rates.epr.median():+.2f} m/yr")
print(f"median 95% CI half-width of LRR: {rates.lrr_ci95.median():.2f} m/yr tide-corrected, "
      f"{rates.lrr_no_tide_ci95.median():.2f} m/yr uncorrected")
print(f"transects with significant retreat (CI excludes 0): {((rates.lrr + rates.lrr_ci95) < 0).sum()} of {len(rates)}")

# Why does the uncorrected rate differ? Because the tide sampled by a sun-synchronous satellite drifts.
tide_trend = stats.linregress(obs.drop_duplicates("scene").year, obs.drop_duplicates("scene").tide)
print(f"trend in predicted tide at overpass time: {tide_trend.slope * 100:+.1f} cm/yr; times the median k "
      f"that is a spurious {rates.k.median() * tide_trend.slope:+.2f} m/yr of apparent seaward shoreline motion")

# --- Step 10: compare with the East Riding cliff-top monitoring -------------------------------
DEFENCES = {"Hornsea": (41, 45), "Mappleton": (51, 52)}
ch = []
for line in lines_utm:
    hit = line.intersection(ref_smooth)
    ch.append(ref_smooth.project(hit if hit.geom_type == "Point" else hit.centroid) / 1000 if not hit.is_empty
              else np.nan)
eryc["chainage_km"] = ch
eryc_all_profiles = eryc.profile.values
eryc = eryc.dropna(subset=["chainage_km"])
# Average the satellite rates over the three transects (300 m) centred on each council profile.
for col in ("lrr", "lrr_no_tide", "epr"):
    eryc[f"s2_{col}"] = [rates.loc[(rates.chainage_km - c).abs() <= 0.15, col].mean() for c in eryc.chainage_km]
cmp_ = eryc.dropna(subset=["rate_1989_2026", "s2_lrr"])
section("comparison with East Riding of Yorkshire Council cliff-top rates")
print(f"{len(cmp_)} council profiles matched to transects")
print(f"council 1989-2026 cliff-top retreat: median {cmp_.rate_1989_2026.median():.2f} m/yr")
print(f"council Mar 2025-Mar 2026 cliff loss: median {cmp_.loss_last_year.median():.2f} m")
print(f"Sentinel-2 2017-2026 waterline retreat (tide-corrected LRR): median {-cmp_.s2_lrr.median():.2f} m/yr")
print(f"Sentinel-2 without tide correction: median {-cmp_.s2_lrr_no_tide.median():.2f} m/yr")
rho = stats.spearmanr(cmp_.rate_1989_2026, -cmp_.s2_lrr)
print(f"Spearman rank correlation, council rate vs Sentinel-2 retreat: {rho.statistic:.2f} (p = {rho.pvalue:.3f})")
diff = -cmp_.s2_lrr - cmp_.rate_1989_2026
print(f"difference (Sentinel-2 minus council): median {diff.median():+.2f} m/yr, MAD {stats.median_abs_deviation(diff):.2f}")
labels = ["Skipsea-Atwick (22-32)", "Atwick-Hornsea (33-46)", "Hornsea-Mappleton (47-59)",
          "Aldbrough (60-74)", "Aldbrough-Withernsea (75-89)"]
cmp_ = cmp_.assign(section=pd.cut(cmp_.profile, [21, 32, 46, 59, 74, 89], labels=labels))
tab = cmp_.groupby("section", observed=True).agg(profiles=("profile", "size"),
                                                council_1989_2026=("rate_1989_2026", "median"),
                                                s2_retreat=("s2_lrr", lambda v: -v.median()),
                                                s2_no_tide=("s2_lrr_no_tide", lambda v: -v.median()))
print(tab.round(2).to_string())

# --- Step 11: save the results ----------------------------------------------------------------
rates.to_csv(OUT / "transect_rates.csv", float_format="%.3f")
gpd.GeoDataFrame(rates, geometry=tr.geometry, crs=CRS).to_crs(4326).to_file(OUT / "transect_rates.geojson")
cmp_.drop(columns=["e1_utm", "n1_utm", "e2_utm", "n2_utm"]).to_csv(OUT / "eryc_comparison.csv", index=False,
                                                                   float_format="%.3f")
obs.to_csv(OUT / "shoreline_positions.csv", index=False, float_format="%.3f")

# --- Step 12: figures -------------------------------------------------------------------------
PLACES = {"Skipsea": [27], "Hornsea": [41, 45], "Mappleton": [51, 52], "Aldbrough": [63], "Withernsea": [89]}
# Alongshore position of every council profile (where it is nearest to the reference line).
prof_ch = pd.Series([ref_smooth.project(l.interpolate(0.5, normalized=True)) / 1000 for l in lines_utm],
                    index=eryc_all_profiles)
place_ch = {k: prof_ch.reindex(v).dropna().mean() for k, v in PLACES.items()}
defence_ch = {k: (prof_ch.get(a), prof_ch.get(b)) for k, (a, b) in DEFENCES.items()}
RED_BLUE = "RdBu"
halo = [pe.withStroke(linewidth=2.5, foreground="white")]

# A clear, recent scene for the true-colour backdrop and the worked example.
ok_scenes = [(i, r) for i, r in zip(items, results) if str(r["status"]) == "ok"]
demo_item, demo_res = max((p for p in ok_scenes if p[0].datetime.year >= 2025),
                          key=lambda p: (round(float(p[1]["zone_clear"]), 2), p[0].datetime))
demo_path = DATA / f"demo_{demo_item.id}.npz"
if not demo_path.exists():
    m, v = scene_mndwi(demo_item)
    signed = planetary_computer.sign(demo_item)
    rgb = np.dstack([reflectance(signed, read_band(signed, b, 1)) for b in ("B04", "B03", "B02")])
    np.savez_compressed(demo_path, mndwi=m.astype("float16"), valid=v,
                        rgb=(np.clip(rgb / 0.25, 0, 1) * 255).astype("uint8"))
demo = dict(np.load(demo_path))
demo_t = pd.Timestamp(demo_item.datetime)
demo_tide = scene_t.set_index("scene").loc[demo_item.id]
print(f"\nexample scene {demo_item.id} ({demo_t:%Y-%m-%d %H:%M} UTC), "
      f"{float(demo_res['zone_clear']):.0%} of shore zone clear, Otsu threshold {float(demo_res['threshold']):.3f}")
ext = (X0, X0 + W * RES, Y0 - H * RES, Y0)
rgb_bg = np.where(in_corridor[..., None], demo["rgb"], 255)


def place_labels(ax, dx=-900):
    for name, c in place_ch.items():
        p = ref_smooth.interpolate(c * 1000)
        ax.plot(p.x, p.y, "o", ms=4, mfc="white", mec="k", zorder=5)
        ax.text(p.x + dx, p.y, name, ha="right", va="center", fontsize=8, path_effects=halo, zorder=6)


def km_axes(ax):
    ax.set_xlabel("Easting (km, UTM 30N)")
    ax.set_ylabel("Northing (km, UTM 30N)")
    ax.xaxis.set_major_formatter(lambda v, _: f"{v / 1000:.0f}")
    ax.yaxis.set_major_formatter(lambda v, _: f"{v / 1000:.0f}")


# Figure 1: study area.
ne_zip = DATA / "ne_50m_admin_0_countries.zip"
if not ne_zip.exists():
    ne_zip.write_bytes(get("https://naciscdn.org/naturalearth/50m/cultural/ne_50m_admin_0_countries.zip",
                           headers={"User-Agent": "Mozilla/5.0"}).content)
world = gpd.read_file(ne_zip)
fig = plt.figure(figsize=(9, 7.6))
ax0 = fig.add_axes([0.0, 0.52, 0.30, 0.45])
world[world.ADMIN.isin(["United Kingdom", "Ireland", "France", "Netherlands", "Belgium"])].plot(
    ax=ax0, color="#e6e6e6", edgecolor="#777", lw=0.4)
ax0.add_patch(plt.Rectangle((bbox_ll[0], bbox_ll[1]), bbox_ll[2] - bbox_ll[0], bbox_ll[3] - bbox_ll[1],
                            fill=False, ec="#d7301f", lw=1.5))
ax0.set_xlim(-8.5, 3); ax0.set_ylim(49.8, 59); ax0.set_xticks([]); ax0.set_yticks([])
ax0.set_xlabel(""); ax0.set_ylabel("")
ax0.set_title("Holderness, East Yorkshire", fontsize=9)
ax1 = fig.add_axes([0.40, 0.04, 0.58, 0.93])
ax1.imshow(rgb_bg, extent=ext)
for g in tr.geometry[::5]:
    ax1.plot(*g.xy, color="#ffd92f", lw=0.8)
ax1.plot(*ref_smooth.xy, color="k", lw=0.8)
for l in lines_utm:
    ax1.plot(*l.xy, color="#e7298a", lw=0.8)
place_labels(ax1)
ax1.plot([], [], color="#ffd92f", label="every 5th Sentinel-2 transect (100 m spacing)")
ax1.plot([], [], color="#e7298a", label="East Riding cliff monitoring profiles")
ax1.plot([], [], color="k", label="reference shoreline (median 2017-2026)")
ax1.legend(loc="lower left", fontsize=7.5, framealpha=0.9)
xb, yb = ref_smooth.bounds[0::2], ref_smooth.bounds[1::2]
ax1.set_xlim(xb[0] - 4000, xb[1] + 2500); ax1.set_ylim(yb[0] - 1500, yb[1] + 1500)
km_axes(ax1)
ax1.set_title(f"Sentinel-2 true colour, {demo_t:%d %B %Y}", fontsize=9)
fig.savefig(OUT / "01-holderness-study-area.png"); plt.close(fig)

# Figure 2: MNDWI, Otsu and the sub-pixel contour for the example scene, zoomed on Aldbrough.
cx, cy = ref_smooth.interpolate(place_ch["Aldbrough"] * 1000).coords[0]
half = 1200
zoom = dict(xlim=(cx - half, cx + half), ylim=(cy - half, cy + half))
dm = np.where(demo["valid"], demo["mndwi"].astype("float32"), np.nan)
c = demo_res["contours"]
fig, axs = plt.subplots(1, 3, figsize=(12.5, 4.3), gridspec_kw={"width_ratios": [1, 1.2, 1.1], "wspace": 0.3})
axs[0].imshow(rgb_bg, extent=ext)
im = axs[1].imshow(dm, extent=ext, cmap="RdBu", vmin=-0.6, vmax=0.6)
fig.colorbar(im, ax=axs[1], shrink=0.8, pad=0.03, label="MNDWI")
for ax in axs[:2]:
    ax.plot(c[:, 0], c[:, 1], color="#ffd92f" if ax is axs[0] else "k", lw=1.2)
    for g in tr.geometry:
        if g.distance(Point(cx, cy)) < half * 1.5:
            ax.plot(*g.xy, color="w" if ax is axs[0] else "#555", lw=0.5)
    ax.set(**zoom); ax.set_xticks([]); ax.set_yticks([])
axs[0].set_title("True colour, shoreline (yellow), transects", fontsize=9)
axs[1].set_title(f"MNDWI and the {float(demo_res['threshold']):.2f} contour", fontsize=9)
zv = dm[zone & np.isfinite(dm)]
axs[2].hist(zv, bins=120, range=(-0.8, 0.8), color="#888")
axs[2].axvline(float(demo_res["threshold"]), color="#d7301f", lw=1.5,
               label=f"Otsu threshold {float(demo_res['threshold']):.3f}")
axs[2].set(xlabel="MNDWI (pixels within 400 m of the coast)", ylabel="pixels")
axs[2].legend(fontsize=8)
axs[2].yaxis.tick_right(); axs[2].yaxis.set_label_position("right")
axs[2].set_title("Peaks: fields, beach sand, sea", fontsize=9)
fig.suptitle(f"Shoreline extraction near Aldbrough, {demo_t:%d %B %Y} {demo_t:%H:%M} UTC", fontsize=10)
fig.savefig(OUT / "02-mndwi-otsu-shoreline.png"); plt.close(fig)

# Figure 3: tide model check and the tides that Sentinel-2 actually samples.
hourly = pd.date_range("2017-01-01", "2026-10-18", freq="h")
mid_w = (tr.lat.median() - GAUGES["Spurn Point"][1]) / (GAUGES["Bridlington"][1] - GAUGES["Spurn Point"][1])
h_all = mid_w * utide.reconstruct(hourly.values, models["Bridlington"], verbose=False).h + \
    (1 - mid_w) * utide.reconstruct(hourly.values, models["Spurn Point"], verbose=False).h
mid_tr = (tr.lat - tr.lat.median()).abs().idxmin()
h_s2 = obs[obs.transect == mid_tr].drop_duplicates("scene").tide
fig = plt.figure(figsize=(11, 6.6))
gs = fig.add_gridspec(2, 3, height_ratios=[1, 1.1], hspace=0.45, wspace=0.25)
for k, day in enumerate(HOLDOUT):
    ax = fig.add_subplot(gs[0, k])
    rms = []
    for name, col in (("Bridlington", "#1b9e77"), ("Spurn Point", "#7570b3")):
        v = val[(val.gauge == name) & (val.day == day)]
        grid = pd.date_range(day, periods=97, freq="15min")
        curve = utide.reconstruct(grid.values, models[name], verbose=False).h + models[name]["offset"]
        ax.plot(v.time, v.observed, ".", ms=3, color=col, label=f"{name} observed")
        ax.plot(grid, curve, "-", lw=1, color=col, label=f"{name} predicted")
        rms.append(f"{name.split()[0]} {np.sqrt((v.error ** 2).mean()):.2f} m")
    ax.xaxis.set_major_locator(mdates.HourLocator(byhour=[0, 6, 12, 18]))
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
    ax.set_xlabel("UTC")
    ax.set_title(f"Hold-out day {day}\nRMSE: " + ", ".join(rms), fontsize=9)
    if k == 0:
        ax.set_ylabel("water level (m, gauge datum)")
        fig.legend(*ax.get_legend_handles_labels(), loc="upper center", ncol=4, fontsize=8, frameon=False,
                   bbox_to_anchor=(0.5, 1.0))
ax = fig.add_subplot(gs[1, :])
bins_ = np.arange(-3.6, 3.0, 0.2)
ax.hist(h_all, bins=bins_, density=True, color="#bbb", label="every hour, 2017-2026 (model)")
ax.hist(h_s2, bins=bins_, density=True, histtype="step", lw=2, color="#d7301f",
        label=f"at the {len(h_s2)} Sentinel-2 overpasses used (~11:10 UTC)")
ax.axvline(np.median(h_s2), color="#d7301f", ls="--", lw=1)
ax.axvline(np.median(h_all), color="#555", ls="--", lw=1)
ax.set(xlabel="predicted tide relative to mean level (m), mid-Holderness", ylabel="density")
ax.legend(fontsize=8)
ax.set_title("The satellite samples the tide unevenly: overpasses fall at a fixed time of day", fontsize=9)
fig.savefig(OUT / "03-tide-model-and-overpass-tides.png"); plt.close(fig)
print(f"predicted tide: median {np.median(h_all):+.2f} m over all hours, {np.median(h_s2):+.2f} m at overpasses; "
      f"5th-95th pct {np.percentile(h_all, 5):+.2f}/{np.percentile(h_all, 95):+.2f} vs "
      f"{np.percentile(h_s2, 5):+.2f}/{np.percentile(h_s2, 95):+.2f} m")

# Figure 4: one transect, before and after tide correction.
ex_j = int((rates.chainage_km - place_ch["Aldbrough"]).abs().idxmin())
g = obs[(obs.transect == ex_j) & obs.keep]
rj = rates.loc[ex_j]
yrs = np.array([g.year.min(), g.year.max()])
fig, axs = plt.subplots(1, 3, figsize=(12, 3.8), gridspec_kw={"width_ratios": [1.4, 1, 1.4]})
sc = axs[0].scatter(g.year, g.x, c=g.tide, cmap="viridis", s=9)
fig.colorbar(sc, ax=axs[0], label="tide (m)")
lr0 = stats.linregress(g.year, g.x)
axs[0].plot(yrs, lr0.intercept + lr0.slope * yrs, "k", lw=1.2)
axs[0].set_title(f"Raw waterline: LRR {rj.lrr_no_tide:+.2f} m/yr, scatter {rj.rmse_no_tide:.0f} m", fontsize=9)
axs[1].scatter(g.tide, g.x - lr0.slope * (g.year - 2017), s=9, color="#555")
tt = np.array([g.tide.min(), g.tide.max()])
axs[1].plot(tt, np.median(g.x - lr0.slope * (g.year - 2017) - rj.k * g.tide) + rj.k * tt, color="#d7301f",
            label=f"k = {rj.k:.1f} m per m (slope 1:{-rj.k:.0f})")
axs[1].set(xlabel="predicted tide (m)", title="Higher tide, waterline further landward")
axs[1].title.set_fontsize(9)
axs[1].legend(fontsize=8)
axs[2].scatter(g.year, g.x_corr, s=9, color="#2c7fb8")
lr1 = stats.linregress(g.year, g.x_corr)
axs[2].plot(yrs, lr1.intercept + lr1.slope * yrs, "k", lw=1.2)
axs[2].set_title(f"Tide-corrected: LRR {rj.lrr:+.2f} ± {rj.lrr_ci95:.2f} m/yr, scatter {rj.rmse_tide:.0f} m",
                 fontsize=9)
for ax in (axs[0], axs[2]):
    ax.set(xlabel="year", ylabel="shoreline position (m, + seaward)")
fig.suptitle(f"Transect {ex_j} at Aldbrough (chainage {rj.chainage_km:.1f} km), {len(g)} Sentinel-2 shorelines",
             fontsize=10)
fig.tight_layout()
fig.savefig(OUT / "04-transect-tide-correction.png"); plt.close(fig)
print(f"example transect {ex_j} ({rj.chainage_km:.1f} km): LRR {rj.lrr:+.2f} ± {rj.lrr_ci95:.2f} m/yr "
      f"(uncorrected {rj.lrr_no_tide:+.2f}), k {rj.k:.1f}, scatter {rj.rmse_no_tide:.0f} -> {rj.rmse_tide:.0f} m")

# Figure 5: along-coast rates against the council's cliff-top rates.
fig = plt.figure(figsize=(11, 7.4))
gs = fig.add_gridspec(2, 2, height_ratios=[1.15, 1], hspace=0.35, wspace=0.25)
ax = fig.add_subplot(gs[0, :])
for name, (a, b) in defence_ch.items():
    if a is not None and b is not None:
        ax.axvspan(min(a, b) - 0.1, max(a, b) + 0.1, color="#ddd", zorder=0)
        ax.text((a + b) / 2, 4.6, f"{name}\ndefences", ha="center", va="top", fontsize=7, color="#555")
ax.fill_between(rates.chainage_km, -(rates.lrr - rates.lrr_ci95), -(rates.lrr + rates.lrr_ci95),
                color="#2c7fb8", alpha=0.25, lw=0)
ax.plot(rates.chainage_km, -rates.lrr, color="#2c7fb8", lw=1.2, label="Sentinel-2 waterline 2017-2026, tide-corrected (95% CI)")
ax.plot(rates.chainage_km, -rates.lrr_no_tide, color="#999", lw=0.8, label="Sentinel-2, no tide correction")
ax.plot(cmp_.chainage_km, cmp_.rate_1989_2026, "o", ms=4, color="#d7301f", label="East Riding cliff top, 1989-2026")
ax.axhline(0, color="k", lw=0.6)
ax.set(xlabel="distance south from Skipsea (km)", ylabel="retreat (m/yr, + landward)", ylim=(-3, 5),
       xlim=(0, rates.chainage_km.max() + 0.05))
for name, c in place_ch.items():
    ax.text(c, -2.8, name, ha="center", fontsize=8, path_effects=halo)
ax.legend(fontsize=7.5, loc="lower left", bbox_to_anchor=(0, 1.01), ncol=3, frameon=False)
ax2 = fig.add_subplot(gs[1, 0])
ax2.errorbar(cmp_.rate_1989_2026, -cmp_.s2_lrr, fmt="o", ms=3.5, color="#2c7fb8")
ax2.plot([0, 4], [0, 4], "k--", lw=0.8, label="1:1")
ax2.set(xlabel="East Riding cliff-top rate 1989-2026 (m/yr)", ylabel="Sentinel-2 waterline retreat (m/yr)",
        xlim=(0, 3.6), ylim=(-1.5, 4))
ax2.set_title(f"62 profiles, Spearman ρ = {rho.statistic:.2f}", fontsize=9)
ax2.legend(fontsize=8)
ax3 = fig.add_subplot(gs[1, 1])
yy_ = np.arange(len(tab))
ax3.barh(yy_ + 0.27, tab.council_1989_2026, 0.27, color="#d7301f", label="East Riding 1989-2026")
ax3.barh(yy_, tab.s2_retreat, 0.27, color="#2c7fb8", label="Sentinel-2 tide-corrected")
ax3.barh(yy_ - 0.27, tab.s2_no_tide, 0.27, color="#999", label="Sentinel-2 uncorrected")
ax3.set_yticks(yy_, [s.split(" (")[0] for s in tab.index], fontsize=8)
ax3.invert_yaxis(); ax3.axvline(0, color="k", lw=0.6)
ax3.set(xlabel="median retreat (m/yr)", xlim=(-0.8, 3.8))
ax3.legend(fontsize=7, loc="upper right")
fig.savefig(OUT / "05-rates-vs-east-riding.png"); plt.close(fig)

# Figure 6: map of the tide-corrected rates.
fig, ax = plt.subplots(figsize=(6.4, 8.6))
ax.imshow(rgb_bg.mean(axis=2), extent=ext, cmap="gray", vmin=0, vmax=255)
norm = TwoSlopeNorm(vmin=-3, vcenter=0, vmax=3)
cm_ = plt.get_cmap(RED_BLUE)
for j, g_ in enumerate(tr.geometry):
    ax.plot(*g_.xy, color=cm_(norm(rates.lrr.iloc[j])), lw=1.6, solid_capstyle="butt")
place_labels(ax)
sm = plt.cm.ScalarMappable(norm=norm, cmap=cm_)
fig.colorbar(sm, ax=ax, shrink=0.6, label="shoreline change 2017-2026 (m/yr; red = retreat)")
ax.set_xlim(xb[0] - 4000, xb[1] + 2500); ax.set_ylim(yb[0] - 1500, yb[1] + 1500)
km_axes(ax)
ax.set_title("Tide-corrected linear regression rate per transect", fontsize=9)
fig.savefig(OUT / "06-shoreline-change-map.png"); plt.close(fig)
print(f"\nfigures and CSV files written to {OUT.name}{chr(92)}")