"""Snow cover mapping with NDSI in Python: the Cairngorms, winter 2025-26, from Sentinel-2 and MODIS.

Searches the Microsoft Planetary Computer STAC API for every Sentinel-2 L2A scene of the Cairngorms
from October 2025 to June 2026, reads only the pixels of the study area, computes the Normalised
Difference Snow Index (NDSI), maps snow with a threshold, checks it against the Scene Classification
Layer (SCL) to show where snow and cloud are confused, tracks snow-covered area through the winter
together with the daily MODIS MOD10A1/MYD10A1 snow product, and uses the Copernicus GLO-30 DEM to
measure snow cover by elevation, the snowline and the effect of aspect.
Run with:  python snow_cover_ndsi.py
"""
import platform
from concurrent.futures import ThreadPoolExecutor
from importlib.metadata import version
from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import planetary_computer
import pystac_client
import rasterio
from matplotlib.colors import BoundaryNorm, ListedColormap, LightSource
from pyproj import Transformer
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.warp import reproject, transform_bounds
from rasterio.windows import from_bounds

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

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


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


# --- Step 1: the study area and a 20 m grid -----------------------------------------------------
# The high Cairngorms: Cairn Gorm, Ben Macdui, Braeriach, Cairn Toul, Beinn a' Bhuird and Ben Avon,
# plus the glens around them, from Glenmore and Rothiemurchus in the north-west to Glen Derry.
BBOX = (-3.90, 56.98, -3.35, 57.20)                  # lon/lat
CRS, RES = "EPSG:32630", 20.0                        # UTM 30N, the native grid of tile 30VVJ; 20 m = B11 pixels
SEASON = "2025-10-01/2026-06-30"
x0, y0, x1, y1 = transform_bounds("EPSG:4326", CRS, *BBOX)
# Snap to the Sentinel-2 20 m grid (tile origins are at ...9960 E / ...0020 N, i.e. multiples of 20 m minus 40).
x0, y1 = np.floor((x0 + 40) / RES) * RES - 40, np.ceil((y1 - 20) / RES) * RES + 20
W, H = int(np.ceil((x1 - x0) / RES)), int(np.ceil((y1 - y0) / RES))
TRANSFORM = from_origin(x0, y1, RES, RES)
BOUNDS = (x0, y1 - H * RES, x0 + W * RES, y1)
section("study area")
print(f"bbox {BBOX} -> {CRS} {W} x {H} pixels of 20 m ({W * H * RES ** 2 / 1e6:.0f} km2)")

catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1",
                                    modifier=planetary_computer.sign_inplace)

# --- Step 2: the Copernicus GLO-30 DEM on the same grid -----------------------------------------
dem_path = CACHE / "dem_20m.npy"
if not dem_path.exists():
    dem = np.full((H, W), np.nan, dtype="float32")
    for item in catalog.search(collections=["cop-dem-glo-30"], bbox=BBOX).items():
        with rasterio.open(item.assets["data"].href) as src:
            tmp = np.full((H, W), np.nan, dtype="float32")
            reproject(rasterio.band(src, 1), tmp, dst_transform=TRANSFORM, dst_crs=CRS,
                      dst_nodata=np.nan, resampling=Resampling.bilinear)
        dem = np.where(np.isnan(dem), tmp, dem)
    np.save(dem_path, dem)
dem = np.load(dem_path)
gy, gx = np.gradient(dem, RES)                       # rows run north to south, so -gy is the northward slope
slope = np.degrees(np.arctan(np.hypot(gx, gy)))
aspect = (np.degrees(np.arctan2(-gx, gy)) + 360) % 360   # direction the slope faces, clockwise from north
hill = LightSource(azdeg=315, altdeg=35).hillshade(dem, dx=RES, dy=RES, vert_exag=1.5)
print(f"DEM: {np.nanmin(dem):.0f}-{np.nanmax(dem):.0f} m, median {np.nanmedian(dem):.0f} m; "
      f"{100 * np.mean(dem > 900):.0f}% of the area is above 900 m")

# --- Step 3: search Sentinel-2 and read only the study area --------------------------------------
items = catalog.search(collections=["sentinel-2-l2a"], bbox=BBOX, datetime=SEASON,
                       query={"s2:mgrs_tile": {"eq": "30VVJ"}}).item_collection()
section("Sentinel-2 search")
print(f"{len(items)} Sentinel-2 L2A scenes of tile 30VVJ, {SEASON}")
BANDS = ["B02", "B03", "B04", "B08", "B11"]
CLEAR_SCL = [2, 4, 5, 6, 7, 11]                      # dark area, vegetation, bare, water, unclassified, snow
CLOUD_SCL = [3, 8, 9, 10]                            # cloud shadow, cloud medium, cloud high, thin cirrus


def read(href: str, resampling=Resampling.average) -> np.ndarray:
    with rasterio.open(href) as src:
        return src.read(1, window=from_bounds(*BOUNDS, src.transform), out_shape=(H, W),
                        resampling=resampling, boundless=True, fill_value=0)


def load_scene(item) -> dict | None:
    f = CACHE / "s2" / f"{item.id}.npz"
    if not f.exists():
        scl = read(item.assets["SCL"].href, Resampling.nearest)
        arrs = {"scl": scl}
        if np.isin(scl, CLEAR_SCL + [11]).mean() >= 0.15:          # skip the band reads for overcast scenes
            # Processing baseline 04.00 (January 2022) and later add 1000 to every digital number.
            offset = 1000 if float(item.properties["s2:processing_baseline"]) >= 4.0 else 0
            for b in BANDS:
                dn = read(item.assets[b].href).astype("float32")
                arrs[b] = np.where(dn > 0, (dn - offset) / 10000, np.nan).astype("float32")
        np.savez_compressed(f, **arrs)
    z = dict(np.load(f))
    z.update(id=item.id, date=pd.Timestamp(item.datetime).tz_localize(None),
             sun_elev=90 - float(item.properties["s2:mean_solar_zenith"]))
    return z


with ThreadPoolExecutor(8) as pool:
    scenes = list(pool.map(load_scene, items))
scenes.sort(key=lambda s: s["date"])

# --- Step 4: NDSI and a snow mask -------------------------------------------------------------------
NDSI_T, NIR_T, GREEN_T = 0.4, 0.11, 0.10             # thresholds from the MODIS SNOMAP algorithm (Hall et al.)


def classify(s: dict, t=NDSI_T) -> dict:
    g, swir, nir = s["B03"], s["B11"], s["B08"]
    with np.errstate(divide="ignore", invalid="ignore"):
        ndsi = (g - swir) / (g + swir)
    snow = (ndsi > t) & (nir > NIR_T) & (g > GREEN_T)
    clear = np.isin(s["scl"], CLEAR_SCL) & np.isfinite(ndsi)
    return {"ndsi": ndsi, "snow": snow & clear, "clear": clear, "snow_any": snow}


rows = []
for s in scenes:
    r = {"date": s["date"], "id": s["id"], "sun_elev": s["sun_elev"],
         "clear_frac": np.isin(s["scl"], CLEAR_SCL).mean(), "scl_snow_frac": (s["scl"] == 11).mean()}
    if "B03" in s:
        c = classify(s)
        nc = c["clear"].sum()
        r.update(clear_frac=nc / c["clear"].size, sca=c["snow"].sum() / max(nc, 1),
                 snow_km2=c["snow"].sum() * RES ** 2 / 1e6,
                 # how the SCL labels the pixels our NDSI test calls snow (cloud-free or not)
                 ndsi_snow_scl_cloud=np.isin(s["scl"], CLOUD_SCL)[c["snow_any"]].mean() if c["snow_any"].any() else np.nan,
                 scl_snow_agree=(c["snow"] == (s["scl"] == 11))[c["clear"]].mean() if nc else np.nan)
    rows.append(r)
s2 = pd.DataFrame(rows)
s2["day"] = s2.date.dt.normalize()
use = s2[s2.clear_frac >= 0.5].copy()
section("Sentinel-2 scenes and snow-covered area")
print(f"scenes with any data: {len(s2)}; with >= 15% clear land (bands read): {s2.sca.notna().sum()}; "
      f"with >= 50% clear (used): {len(use)}")
print(f"sun elevation at overpass: {s2.sun_elev.min():.1f} deg (December) to {s2.sun_elev.max():.1f} deg (June)")
print(use[["date", "sun_elev", "clear_frac", "sca", "snow_km2", "scl_snow_agree", "ndsi_snow_scl_cloud"]]
      .assign(date=use.date.dt.strftime("%Y-%m-%d %H:%M")).round(3).to_string(index=False))

# Showcase scene: the clearest scene with plenty of snow.
best = use[use.sca > 0.3].sort_values("clear_frac", ascending=False).iloc[0]
S = next(s for s in scenes if s["id"] == best.id)
C = classify(S)
section(f"showcase scene {best.date:%Y-%m-%d}")
print(f"{best.id}: clear {100 * best.clear_frac:.1f}%, snow-covered area {100 * best.sca:.1f}% "
      f"({best.snow_km2:.1f} km2), sun elevation {best.sun_elev:.1f} deg")
nd = C["ndsi"][C["clear"]]
print(f"NDSI of clear pixels: 5th/50th/95th percentile {np.nanpercentile(nd, 5):.2f} / {np.nanpercentile(nd, 50):.2f} / "
      f"{np.nanpercentile(nd, 95):.2f}")
sens = {t: (classify(S, t)["snow"].sum() / C["clear"].sum()) for t in np.round(np.arange(0.1, 0.75, 0.05), 2)}
print("snow-covered area vs NDSI threshold: " + ", ".join(f"{t:.2f}: {100 * v:.1f}%" for t, v in sens.items()))
ndsi_only = (C["ndsi"] > NDSI_T) & C["clear"]
print(f"pixels with NDSI > 0.4 removed by the NIR > 0.11 / green > 0.10 tests: {(ndsi_only & ~C['snow']).sum():,} "
      f"({100 * (ndsi_only & ~C['snow']).sum() / max(ndsi_only.sum(), 1):.1f}%), of which in water (SCL 6): "
      f"{(ndsi_only & ~C['snow'] & (S['scl'] == 6)).sum():,}")


def rgb(s: dict, gamma=0.7) -> np.ndarray:
    a = np.dstack([s["B04"], s["B03"], s["B02"]])
    return np.clip(np.nan_to_num(a / 0.9), 0, 1) ** gamma


def false_colour(s: dict) -> np.ndarray:                # SWIR1, NIR, red: snow cyan/blue, cloud white
    a = np.dstack([s["B11"], s["B08"], s["B04"]])
    return np.clip(np.nan_to_num(a / 0.8), 0, 1) ** 0.8


ext = [BOUNDS[0], BOUNDS[2], BOUNDS[1], BOUNDS[3]]
to_utm = Transformer.from_crs("EPSG:4326", CRS, always_xy=True)
PEAKS = {"Cairn Gorm": (-3.644, 57.117), "Ben Macdui": (-3.669, 57.070), "Braeriach": (-3.728, 57.078),
         "Beinn a' Bhuird": (-3.499, 57.087), "Ben Avon": (-3.434, 57.099), "Aviemore": (-3.829, 57.193)}


def label_peaks(ax, colour="black"):
    for name, (lo, la) in PEAKS.items():
        px, py = to_utm.transform(lo, la)
        if BOUNDS[0] < px < BOUNDS[2] and BOUNDS[1] < py < BOUNDS[3]:
            ax.plot(px, py, "^", ms=4, color=colour)
            ax.annotate(name, (px, py), xytext=(3, 3), textcoords="offset points", fontsize=7, color=colour)


def km_axes(ax):
    ax.set_xticks([]); ax.set_yticks([])


# --- Figure 1: true colour, NDSI and snow mask ---------------------------------------------------
fig, axs = plt.subplots(1, 3, figsize=(12, 3.9), layout="constrained")
axs[0].imshow(rgb(S), extent=ext); axs[0].set_title(f"True colour, {best.date:%d %B %Y}")
im = axs[1].imshow(C["ndsi"], extent=ext, cmap="RdBu", vmin=-0.8, vmax=1)
fig.colorbar(im, ax=axs[1], shrink=0.8, label="NDSI = (B03 − B11) / (B03 + B11)")
axs[1].set_title("NDSI")
axs[2].imshow(hill, extent=ext, cmap="gray", vmin=0, vmax=1)
cls = np.where(~C["clear"], 2, np.where(C["snow"], 1, 0)).astype(float)
axs[2].imshow(cls, extent=ext, cmap=ListedColormap(["#a6611a", "#2c7bb6", "#bbbbbb"]), alpha=0.6,
              norm=BoundaryNorm([-0.5, 0.5, 1.5, 2.5], 3), interpolation="nearest")
axs[2].set_title(f"Snow (blue): {100 * best.sca:.0f}% of clear land; grey = cloud/no data")
for ax in axs:
    km_axes(ax)
label_peaks(axs[2])
fig.savefig(OUT / "01-cairngorms-ndsi-snow-map.png"); plt.close(fig)

# --- Figure 2: choosing the NDSI threshold -------------------------------------------------------
fig, axs = plt.subplots(1, 2, figsize=(11, 3.8), layout="constrained")
bins = np.linspace(-0.6, 1, 81)
scl_c = S["scl"][C["clear"]]
for code, lab, col in [(4, "SCL vegetation", "#1a9641"), (5, "SCL bare soil / rock", "#a6611a"),
                       (11, "SCL snow", "#2c7bb6"), (6, "SCL water", "#7b3294"), (2, "SCL dark area (shadow)", "#555555")]:
    v = nd[scl_c == code]
    if v.size:
        axs[0].hist(v, bins=bins, histtype="step", lw=1.4, color=col, label=f"{lab} ({v.size:,})")
axs[0].axvline(NDSI_T, color="#d7191c", ls="--"); axs[0].text(NDSI_T + 0.02, axs[0].get_ylim()[1] * 0.9, "0.4", color="#d7191c")
axs[0].set_xlabel("NDSI"); axs[0].set_ylabel("pixels"); axs[0].set_yscale("log")
axs[0].set_title("NDSI of clear pixels by Sen2Cor scene class"); axs[0].legend(fontsize=7, frameon=False)
axs[1].plot(list(sens), [100 * v for v in sens.values()], "o-", color="#2c7bb6", mec="black", mew=0.3)
axs[1].axvline(NDSI_T, color="#d7191c", ls="--")
axs[1].set_xlabel("NDSI threshold"); axs[1].set_ylabel("snow-covered area (% of clear land)")
axs[1].set_title("Sensitivity of snow-covered area to the threshold"); axs[1].grid(alpha=0.3)
fig.savefig(OUT / "02-ndsi-threshold.png"); plt.close(fig)

# --- Step 5: snow and cloud confusion -------------------------------------------------------------
# A scene with both snow and cloud: the largest share of NDSI snow that the SCL calls cloud.
mixed = s2[(s2.clear_frac.between(0.3, 0.8)) & (s2.sca > 0.2)].sort_values("ndsi_snow_scl_cloud", ascending=False)
M = next(s for s in scenes if s["id"] == mixed.iloc[0].id)
MC = classify(M)
section(f"snow/cloud confusion, {mixed.iloc[0].date:%Y-%m-%d}")
scl_names = {0: "no data", 1: "saturated", 2: "dark area", 3: "cloud shadow", 4: "vegetation", 5: "bare",
             6: "water", 7: "unclassified", 8: "cloud medium", 9: "cloud high", 10: "thin cirrus", 11: "snow"}
tab = pd.Series(M["scl"][MC["snow_any"]]).map(scl_names).value_counts(normalize=True).mul(100).round(1)
print("SCL class of the pixels the NDSI test calls snow (%):")
print(tab.to_string())
cloud = np.isin(M["scl"], [8, 9])
print(f"NDSI of SCL cloud (medium/high) pixels: median {np.nanmedian(MC['ndsi'][cloud]):.2f}; "
      f"{100 * np.nanmean(MC['ndsi'][cloud] > NDSI_T):.1f}% exceed 0.4")
print(f"B11 (SWIR) reflectance: SCL snow median {np.nanmedian(M['B11'][M['scl'] == 11]):.3f}, "
      f"SCL cloud median {np.nanmedian(M['B11'][cloud]):.3f}")
print("\nacross all scenes read: share of NDSI-snow pixels that the SCL labels cloud or shadow: "
      f"median {100 * s2.ndsi_snow_scl_cloud.median():.1f}%, max {100 * s2.ndsi_snow_scl_cloud.max():.1f}%")
print(f"agreement of our snow mask with SCL class 11 on clear pixels (scenes used): "
      f"median {100 * use.scl_snow_agree.median():.1f}%, min {100 * use.scl_snow_agree.min():.1f}%")

scl_cols = {0: "#000000", 1: "#ff0000", 2: "#404040", 3: "#7f5f3f", 4: "#1a9641", 5: "#d8b365", 6: "#2c7bb6",
            7: "#999999", 8: "#cccccc", 9: "#ffffff", 10: "#f4a6d7", 11: "#5ee0f2"}
fig, axs = plt.subplots(1, 4, figsize=(14, 3.7), layout="constrained")
axs[0].imshow(rgb(M), extent=ext); axs[0].set_title(f"True colour, {mixed.iloc[0].date:%d %b %Y}")
axs[1].imshow(false_colour(M), extent=ext); axs[1].set_title("False colour B11/B08/B04:\nsnow blue-cyan, cloud white")
cmap_scl = ListedColormap([scl_cols[i] for i in range(12)])
axs[2].imshow(M["scl"], extent=ext, cmap=cmap_scl, vmin=-0.5, vmax=11.5, interpolation="nearest")
axs[2].set_title("Sen2Cor scene classification (SCL)")
handles = [plt.Rectangle((0, 0), 1, 1, fc=scl_cols[i], ec="#555555", lw=0.4) for i in (2, 4, 5, 8, 9, 10, 11)]
axs[2].legend(handles, [scl_names[i] for i in (2, 4, 5, 8, 9, 10, 11)], fontsize=6, loc="lower left", framealpha=0.85)
conf = np.zeros(M["scl"].shape, dtype=float)
conf[MC["snow_any"] & ~np.isin(M["scl"], CLOUD_SCL)] = 1                   # NDSI snow, SCL clear
conf[MC["snow_any"] & np.isin(M["scl"], CLOUD_SCL)] = 2                    # NDSI snow, SCL cloud
conf[~MC["snow_any"] & np.isin(M["scl"], CLOUD_SCL)] = 3                   # cloud only
axs[3].imshow(hill, extent=ext, cmap="gray", vmin=0, vmax=1)
axs[3].imshow(np.where(conf == 0, np.nan, conf), extent=ext, interpolation="nearest", alpha=0.75,
              cmap=ListedColormap(["#2c7bb6", "#d7191c", "#bbbbbb"]), norm=BoundaryNorm([0.5, 1.5, 2.5, 3.5], 3))
axs[3].set_title("NDSI snow vs SCL cloud")
axs[3].legend([plt.Rectangle((0, 0), 1, 1, fc=c) for c in ("#2c7bb6", "#d7191c", "#bbbbbb")],
              ["NDSI snow, SCL clear", "NDSI snow, SCL cloud", "SCL cloud only"], fontsize=6, loc="lower left", framealpha=0.85)
for ax in axs:
    km_axes(ax)
fig.savefig(OUT / "03-snow-cloud-confusion.png"); plt.close(fig)

# --- Step 6: the daily MODIS snow product -----------------------------------------------------------
section("MODIS MOD10A1 / MYD10A1 daily snow cover")
SIN = 'PROJCS["sinusoidal",GEOGCS["sphere",DATUM["sphere",SPHEROID["sphere",6371007.181,0]],PRIMEM["Greenwich",0],' \
      'UNIT["degree",0.0174532925199433]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],' \
      'PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1]]'
modis_items = catalog.search(collections=["modis-10A1-061"], bbox=BBOX, datetime=SEASON).item_collection()
mb = transform_bounds(CRS, SIN, *BOUNDS)


def load_modis(item) -> tuple:
    f = CACHE / "modis" / f"{item.id}.npy"
    if not f.exists():
        with rasterio.open(item.assets["NDSI_Snow_Cover"].href) as src:
            win = from_bounds(*mb, src.transform).round_offsets().round_lengths()
            np.save(f, src.read(1, window=win))
            if not (CACHE / "modis_transform.npy").exists():
                np.save(CACHE / "modis_transform.npy", np.array(src.window_transform(win))[:6])
    return item.id, pd.Timestamp(item.properties["start_datetime"]).tz_localize(None).normalize(), np.load(f)


with ThreadPoolExecutor(16) as pool:
    mods = list(pool.map(load_modis, modis_items))
mt = rasterio.Affine(*np.load(CACHE / "modis_transform.npy"))
shape_m = mods[0][2].shape
# DEM averaged onto the 463 m MODIS grid, and a mask of the pixels inside the study area.
dem_m = np.full(shape_m, np.nan, dtype="float32")
reproject(dem, dem_m, src_transform=TRANSFORM, src_crs=CRS, dst_transform=mt, dst_crs=SIN,
          src_nodata=np.nan, dst_nodata=np.nan, resampling=Resampling.average)
terra = {d: a for i, d, a in mods if i.startswith("MOD")}
aqua = {d: a for i, d, a in mods if i.startswith("MYD")}
print(f"{len(terra)} Terra and {len(aqua)} Aqua days on a {shape_m[1]} x {shape_m[0]} grid of 463 m pixels")
mrows, mdaily = [], {}
for d in sorted(set(terra) | set(aqua)):
    t = terra.get(d, np.full(shape_m, 250, "uint8"))
    a = aqua.get(d, np.full(shape_m, 250, "uint8"))
    comb = np.where(t <= 100, t, np.where(a <= 100, a, t))           # Terra first, Aqua where Terra is cloudy
    ok = (comb <= 100) & np.isfinite(dem_m)
    snow = ok & (comb >= 40)                                           # NDSI >= 0.4 (values are NDSI x 100)
    mdaily[d] = (comb, ok, snow)
    mrows.append({"day": d, "clear_terra": np.mean((t <= 100)[np.isfinite(dem_m)]),
                  "clear": ok.sum() / np.isfinite(dem_m).sum(), "sca": snow.sum() / max(ok.sum(), 1)})
modis = pd.DataFrame(mrows)
mclear = modis[modis.clear >= 0.7]
print(f"days with >= 70% clear (Terra or Aqua): {len(mclear)} of {len(modis)}; Terra alone would give "
      f"{(modis.clear_terra >= 0.7).sum()}")
print("monthly mean snow-covered area on clear days (%):")
print(mclear.set_index("day").sca.resample("MS").agg(["mean", "count"]).assign(mean=lambda d: (100 * d["mean"]).round(1))
      .rename(index=lambda i: i.strftime("%Y-%m")).T.to_string())

both = use.merge(modis, on="day", suffixes=("_s2", "_modis"))
both = both[both.clear >= 0.7]
print(f"\nsame-day Sentinel-2 and MODIS (both >= 50%/70% clear): {len(both)} days; "
      f"SCA difference MODIS - S2: mean {100 * (both.sca_modis - both.sca_s2).mean():+.1f} points, "
      f"r = {both[['sca_s2', 'sca_modis']].corr().iloc[0, 1]:.2f}")

# --- Figure 4: snow-covered area through the winter ------------------------------------------------
fig, ax = plt.subplots(figsize=(11, 4.2))
ax.bar(modis.day, 100 * modis.clear, width=1, color="#dddddd", label="MODIS clear-sky share of the area")
ax.plot(mclear.day, 100 * mclear.sca, ".", color="#fdae61", ms=5, label="MODIS daily SCA (≥ 70% clear)")
med = mclear.set_index("day").sca.rolling("15D").median()
gaps = med.index.to_series().diff() > pd.Timedelta("15D")                 # break the line across long cloudy gaps
med = pd.concat([med, pd.Series(np.nan, index=med.index[gaps] - pd.Timedelta("1D"))]).sort_index()
ax.plot(med.index, 100 * med, color="#d7191c", lw=1.3, label="MODIS 15-day median")
ax.plot(use.date, 100 * use.sca, "o", color="#2c7bb6", mec="black", mew=0.4, ms=6, label="Sentinel-2 SCA (≥ 50% clear)")
ax.set_ylabel("snow-covered area (% of clear land)"); ax.set_ylim(0, 100)
ax.set_title("Snow-covered area of the Cairngorms study area, winter 2025–26")
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b\n%Y")); ax.legend(fontsize=7, ncol=2, loc="upper right", frameon=True)
fig.savefig(OUT / "04-snow-covered-area-season.png"); plt.close(fig)

# --- Step 7: snow cover by elevation and the snowline ------------------------------------------------
section("snow cover by elevation and the snowline")
EDGES = np.arange(200, 1351, 50)
MID = (EDGES[:-1] + EDGES[1:]) / 2


def snow_by_band(snow, clear, z, edges=EDGES, min_px=200):
    zb = np.digitize(z, edges) - 1
    out = np.full(len(edges) - 1, np.nan)
    for i in range(len(edges) - 1):
        m = clear & (zb == i)
        if m.sum() >= min_px:
            out[i] = snow[m].mean()
    return out


def snowline(frac, mids=MID, level=0.5):
    """Lowest elevation above which the snow fraction stays at or above `level` (linear interpolation)."""
    ok = np.isfinite(frac)
    f, z = frac[ok], mids[ok]
    if len(f) < 3 or f[-1] < level:
        return np.nan                                  # no snowline: the top band is not snow-covered
    above = f >= level
    i = len(f) - 1
    while i > 0 and above[i - 1]:
        i -= 1
    if i == 0:
        return z[0]                                    # snow down to the lowest band with data
    return z[i - 1] + (level - f[i - 1]) / (f[i] - f[i - 1]) * (z[i] - z[i - 1])


north = (aspect >= 315) | (aspect < 45)
south = (aspect >= 135) & (aspect < 225)
srows, profiles = [], {}
for _, r in use.iterrows():
    s = next(x for x in scenes if x["id"] == r.id)
    c = classify(s)
    frac = snow_by_band(c["snow"], c["clear"], dem)
    profiles[r.date] = frac
    srows.append({"date": r.date, "sca": r.sca, "snowline": snowline(frac),
                  "snowline_N": snowline(snow_by_band(c["snow"], c["clear"] & north & (slope > 5), dem, min_px=100)),
                  "snowline_S": snowline(snow_by_band(c["snow"], c["clear"] & south & (slope > 5), dem, min_px=100))})
sl = pd.DataFrame(srows)
mprof = []
EDGES_M = np.arange(200, 1351, 100)
for d in mclear.day:
    comb, ok, snow = mdaily[d]
    frac = snow_by_band(snow, ok, dem_m, EDGES_M, min_px=8)
    mprof.append({"day": d, "snowline": snowline(frac, (EDGES_M[:-1] + EDGES_M[1:]) / 2)})
msl = pd.DataFrame(mprof)
print(sl.assign(date=sl.date.dt.strftime("%Y-%m-%d")).round({"sca": 3, "snowline": 0, "snowline_N": 0, "snowline_S": 0})
      .to_string(index=False))
d_ns = (sl.snowline_S - sl.snowline_N).dropna()
print(f"\nsouth-facing minus north-facing snowline: median {d_ns.median():.0f} m over {len(d_ns)} scenes "
      f"(range {d_ns.min():.0f} to {d_ns.max():.0f} m)")
print(f"MODIS snowline (100 m bands) on {msl.snowline.notna().sum()} clear days; monthly median (m):")
print(msl.set_index("day").snowline.resample("MS").median().round(0).rename(index=lambda i: i.strftime("%Y-%m"))
      .to_frame().T.to_string())

# --- Figure 5: snow fraction against elevation, and aspect ------------------------------------------
pick = [d for d in profiles if np.isfinite(profiles[d]).sum() >= 10]
pick = [pick[i] for i in np.unique(np.linspace(0, len(pick) - 1, min(7, len(pick))).round().astype(int))]
fig, axs = plt.subplots(1, 2, figsize=(11, 4.3), layout="constrained", width_ratios=[1.2, 1])
cm = plt.get_cmap("viridis")
for i, d in enumerate(pick):
    axs[0].plot(100 * profiles[d], MID, "-o", ms=3, color=cm(i / max(len(pick) - 1, 1)), label=f"{d:%d %b %Y}")
axs[0].axvline(50, color="black", lw=0.8, ls="--")
axs[0].set_xlabel("snow-covered share of clear pixels (%)"); axs[0].set_ylabel("elevation (m, Copernicus DEM)")
axs[0].set_title("Snow cover by 50 m elevation band (Sentinel-2)"); axs[0].legend(fontsize=7, frameon=False)
axs[0].grid(alpha=0.3)
dem_valid = dem[np.isfinite(dem)]
ax2 = axs[0].twiny()
ax2.hist(dem_valid, bins=EDGES, orientation="horizontal", color="#cccccc", alpha=0.5, zorder=0)
ax2.set_xlim(ax2.get_xlim()[1] * 4, 0); ax2.set_xticks([])
axs[0].set_zorder(ax2.get_zorder() + 1); axs[0].patch.set_visible(False)    # profiles above the histogram
axs[1].plot(sl.date, sl.snowline_N, "o", color="#2c7bb6", mec="black", mew=0.3, label="north-facing (315–45°)")
axs[1].plot(sl.date, sl.snowline_S, "s", color="#d7191c", mec="black", mew=0.3, label="south-facing (135–225°)")
for _, r in sl.dropna(subset=["snowline_N", "snowline_S"]).iterrows():
    axs[1].plot([r.date, r.date], [r.snowline_N, r.snowline_S], color="#999999", lw=0.8, zorder=0)
axs[1].set_ylabel("snowline elevation (m)"); axs[1].set_title("Snowline on north- and south-facing slopes")
axs[1].xaxis.set_major_formatter(mdates.DateFormatter("%b\n%Y")); axs[1].legend(fontsize=7, frameon=False); axs[1].grid(alpha=0.3)
fig.savefig(OUT / "05-snow-by-elevation-aspect.png"); plt.close(fig)

# --- Figure 6: the snowline through the winter ------------------------------------------------------
fig, ax = plt.subplots(figsize=(11, 4))
ax.plot(msl.day, msl.snowline, ".", color="#fdae61", ms=5, label="MODIS daily snowline (100 m bands)")
ax.plot(msl.set_index("day").snowline.rolling("15D", min_periods=2).median(), color="#d7191c", lw=1.3,
        label="MODIS 15-day median")
ax.plot(sl.date, sl.snowline, "o", color="#2c7bb6", mec="black", mew=0.4, ms=6, label="Sentinel-2 snowline (50 m bands)")
ax.axhline(np.nanmax(dem), color="#555555", lw=0.8, ls=":")
ax.text(msl.day.min(), np.nanmax(dem) + 15, f"highest point in the DEM, {np.nanmax(dem):.0f} m", fontsize=7, color="#555555")
ax.set_ylabel("snowline elevation (m)"); ax.set_ylim(150, 1400)
ax.set_title("Cairngorms snowline (lowest elevation where at least half of the clear land is snow-covered)")
ax.xaxis.set_major_formatter(mdates.DateFormatter("%b\n%Y")); ax.legend(fontsize=7, frameon=True, loc="lower center", ncol=3)
fig.savefig(OUT / "06-snowline-season.png"); plt.close(fig)

# --- Step 8: save the results ------------------------------------------------------------------------
s2.drop(columns=["day"]).to_csv(OUT / "sentinel2_scenes.csv", index=False, float_format="%.4f")
modis.to_csv(OUT / "modis_daily_sca.csv", index=False, float_format="%.4f")
sl.to_csv(OUT / "snowline_sentinel2.csv", index=False, float_format="%.1f")
with rasterio.open(OUT / f"snow_{best.date:%Y%m%d}.tif", "w", driver="GTiff", height=H, width=W, count=1, dtype="uint8",
                   crs=CRS, transform=TRANSFORM, nodata=255, compress="deflate") as dst:
    dst.write(np.where(C["clear"], C["snow"].astype("uint8"), 255), 1)
print(f"\nwrote output/sentinel2_scenes.csv, modis_daily_sca.csv, snowline_sentinel2.csv and snow_{best.date:%Y%m%d}.tif")