Source code: Snow cover mapping with NDSI in Python: Sentinel-2 and MODIS over the Cairngorms
Runnable Python scripts, requirements and console output for the tutorial “Snow cover mapping with NDSI in Python: Sentinel-2 and MODIS over the…
Last updated
These are the exact scripts used to produce every figure and number in Snow cover mapping with NDSI in Python: Sentinel-2 and MODIS over the Cairngorms. Download them, install the requirements into a fresh virtual environment and run them from this folder.
- snow_cover_ndsi.py
- requirements.txt
- console.txt: the console output recorded when we ran it
python -m venv .venv.venv\Scripts\activate # macOS/Linux: source .venv/bin/activatepip install -r requirements.txtpython snow_cover_ndsi.pysnow_cover_ndsi.py
"""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 Cairngormsfrom October 2025 to June 2026, reads only the pixels of the study area, computes the NormalisedDifference Snow Index (NDSI), maps snow with a threshold, checks it against the Scene ClassificationLayer (SCL) to show where snow and cloud are confused, tracks snow-covered area through the wintertogether with the daily MODIS MOD10A1/MYD10A1 snow product, and uses the Copernicus GLO-30 DEM tomeasure snow cover by elevation, the snowline and the effect of aspect.Run with: python snow_cover_ndsi.py"""import platformfrom concurrent.futures import ThreadPoolExecutorfrom importlib.metadata import versionfrom pathlib import Path
import matplotlib
matplotlib.use("Agg")import matplotlib.dates as mdatesimport matplotlib.pyplot as pltimport numpy as npimport pandas as pdimport planetary_computerimport pystac_clientimport rasteriofrom matplotlib.colors import BoundaryNorm, ListedColormap, LightSourcefrom pyproj import Transformerfrom rasterio.enums import Resamplingfrom rasterio.transform import from_originfrom rasterio.warp import reproject, transform_boundsfrom rasterio.windows import from_bounds
HERE = Path(__file__).parentOUT, 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/latCRS, RES = "EPSG:32630", 20.0 # UTM 30N, the native grid of tile 30VVJ; 20 m = B11 pixelsSEASON = "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 + 20W, 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 slopeslope = np.degrees(np.arctan(np.hypot(gx, gy)))aspect = (np.degrees(np.arctan2(-gx, gy)) + 360) % 360 # direction the slope faces, clockwise from northhill = 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, snowCLOUD_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 clearconf[MC["snow_any"] & np.isin(M["scl"], CLOUD_SCL)] = 2 # NDSI snow, SCL cloudconf[~MC["snow_any"] & np.isin(M["scl"], CLOUD_SCL)] = 3 # cloud onlyaxs[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 gapsmed = 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 histogramaxs[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")