Source code: Burn severity mapping with Sentinel-2 and dNBR in Python (Saddleworth Moor 2018)

"""Wildfire burn severity mapping with dNBR in Python: the 2018 Saddleworth Moor fire.

Streams Sentinel-2 L2A scenes from the Microsoft Planetary Computer, masks cloud, shadow and water
with the scene classification layer (SCL), computes the Normalised Burn Ratio (NBR) before and after
the fire, the difference (dNBR) with an offset from unburned ground, the USGS (Key and Benson, 2006)
severity classes and their area in hectares, a burn perimeter, the fire's progression, the relative
dNBR (RdNBR) and an eight-year recovery time series.
Run with:  python burn_severity_dnbr.py
"""
import platform
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.pyplot as plt
import numpy as np
import pandas as pd
import planetary_computer
import pystac_client
import rasterio
import rasterio.features
from matplotlib.colors import BoundaryNorm, ListedColormap
from matplotlib.patches import Patch
from pyproj import Transformer
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.windows import from_bounds
from scipy import ndimage

HERE = Path(__file__).parent
OUT, CACHE = HERE / "output", HERE / "cache"
for d in (OUT, CACHE):
    d.mkdir(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__}")


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


# --- Step 1: study area, analysis grid and the pre- and post-fire scenes ------------------
STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
AOI = (-2.06, 53.475, -1.90, 53.545)   # lon/lat: Saddleworth Moor, Buckton Moor and Dovestone, Peak District
CRS = "EPSG:32630"                     # UTM zone 30N, the native grid of Sentinel-2 tile 30UWE
TILE = "30UWE"
RES = 20                               # B8A and B12 are both 20 m bands
PRE, POST = "2018-06-24", "2018-07-14"
PROGRESSION = ["2018-06-26", "2018-06-29", "2018-07-04", "2018-07-14"]

to_utm = Transformer.from_crs("EPSG:4326", CRS, always_xy=True)
xs, ys = to_utm.transform([AOI[0], AOI[2]], [AOI[1], AOI[3]])
x0, y0 = np.floor(min(xs) / RES) * RES, np.floor(min(ys) / RES) * RES   # snap to the native 20 m grid
x1, y1 = np.ceil(max(xs) / RES) * RES, np.ceil(max(ys) / RES) * RES
WIDTH, HEIGHT = int((x1 - x0) / RES), int((y1 - y0) / RES)
TRANSFORM = from_origin(x0, y1, RES, RES)
BOUNDS = (x0, y0, x1, y1)
HA_PER_PX = RES * RES / 10_000          # one 20 m pixel = 0.04 ha

catalog = pystac_client.Client.open(STAC_URL)
search = catalog.search(collections=["sentinel-2-l2a"], bbox=AOI, datetime="2018-06-01/2018-08-31",
                        query={"s2:mgrs_tile": {"eq": TILE}, "eo:cloud_cover": {"lt": 50}})
summer = {it.datetime.strftime("%Y-%m-%d"): it for it in search.items()}
section("STAC search")
print(f"grid: {WIDTH} x {HEIGHT} pixels of {RES} m ({(x1 - x0) / 1000:.1f} x {(y1 - y0) / 1000:.1f} km), {CRS}")
print(f"{len(summer)} scenes of tile {TILE} with < 50% cloud, June to August 2018")
for d in [PRE] + PROGRESSION:
    p = summer[d].properties
    print(f"{d}  {p['platform']}  scene cloud {p['eo:cloud_cover']:5.1f}%  baseline {p['s2:processing_baseline']}"
          f"  sun zenith {p['s2:mean_solar_zenith']:.1f}")

# --- Step 2: read bands, convert to reflectance and mask cloud, shadow and water -------------
BAD_SCL = [0, 1, 3, 8, 9, 10, 11]   # no data, saturated, cloud shadow, cloud medium/high, cirrus, snow
WATER_SCL = 6
GDAL_ENV = dict(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", GDAL_HTTP_MAX_RETRY=5, GDAL_HTTP_RETRY_DELAY=2)


def read_band(item, band: str, bounds=BOUNDS) -> np.ndarray:
    """Read one band for a window of the native 20 m grid (10 m bands are averaged to 20 m)."""
    href = planetary_computer.sign(item.assets[band].href)   # sign at read time: tokens expire
    shape = (int(round((bounds[3] - bounds[1]) / RES)), int(round((bounds[2] - bounds[0]) / RES)))
    for attempt in range(4):
        try:
            with rasterio.Env(**GDAL_ENV), rasterio.open(href) as src:
                return src.read(1, window=from_bounds(*bounds, src.transform), out_shape=shape, boundless=True,
                                fill_value=0, resampling=Resampling.nearest if band == "SCL" else Resampling.average)
        except rasterio.errors.RasterioIOError:
            time.sleep(5 * (attempt + 1))
    raise RuntimeError(f"could not read {band} of {item.id}")


def to_reflectance(dn: np.ndarray, item) -> np.ndarray:
    """L2A digital numbers to surface reflectance. Baseline 04.00+ (from 25 Jan 2022) adds a -1000 offset."""
    offset = -1000 if float(item.properties["s2:processing_baseline"]) >= 4.0 else 0
    return np.where(dn > 0, (dn.astype("float32") + offset) / 10_000, np.nan)


def nbr(nir: np.ndarray, swir2: np.ndarray) -> np.ndarray:
    """Normalised Burn Ratio = (NIR - SWIR2) / (NIR + SWIR2), with B8A and B12 for Sentinel-2."""
    with np.errstate(divide="ignore", invalid="ignore"):
        return (nir - swir2) / (nir + swir2)


def load_scene(date: str, bands=("B8A", "B12", "B04", "B03", "B02")) -> dict:
    """Reflectance, NBR and a clear-land mask for one date, cached as .npz."""
    f = CACHE / f"s2_{date}_{len(bands)}.npz"
    if not f.exists():
        item = summer[date]
        arr = {b: to_reflectance(read_band(item, b), item) for b in bands}
        arr["SCL"] = read_band(item, "SCL")
        np.savez_compressed(f, **arr)
    s = dict(np.load(f))
    cloud = ndimage.binary_dilation(np.isin(s["SCL"], BAD_SCL[1:]), iterations=3)   # 60 m buffer round cloud
    s["clear"] = ~cloud & (s["SCL"] != 0) & (s["SCL"] != WATER_SCL) & np.isfinite(s["B12"])
    s["NBR"] = nbr(s["B8A"], s["B12"])
    return s


pre, post = load_scene(PRE), load_scene(POST)
valid = pre["clear"] & post["clear"]
section("cloud and water masking")
for name, s in (("pre ", pre), ("post", post)):
    scl = s["SCL"]
    print(f"{name} SCL: vegetation {np.mean(scl == 4):.1%}, not vegetated {np.mean(scl == 5):.1%}, "
          f"water {np.mean(scl == 6):.1%}, cloud/shadow/cirrus {np.isin(scl, BAD_SCL[1:]).mean():.1%}")
print(f"pixels clear on both dates: {valid.sum():,} of {valid.size:,} ({valid.mean():.1%})")

# --- Step 3: NBR and dNBR ------------------------------------------------------------------
dnbr_raw = np.where(valid, pre["NBR"] - post["NBR"], np.nan)
section("NBR and dNBR")
print(f"median NBR: pre-fire {np.nanmedian(pre['NBR'][valid]):.3f}, post-fire {np.nanmedian(post['NBR'][valid]):.3f}")
print(f"dNBR 0.5th-99.5th percentile: {np.nanpercentile(dnbr_raw, 0.5):.3f} to {np.nanpercentile(dnbr_raw, 99.5):.3f}")

# A first, generous burn mask finds the scar so that we can measure the unburned background around it.
seed = ndimage.binary_opening(np.nan_to_num(dnbr_raw) > 0.27, iterations=1)
lab, n = ndimage.label(seed)
sizes = ndimage.sum(seed, lab, index=np.arange(1, n + 1))
scar = np.isin(lab, 1 + np.flatnonzero(sizes * HA_PER_PX >= 50))          # patches of at least 50 ha
near = ndimage.binary_dilation(scar, iterations=50)                       # within 1 km of the scar
ring = near & ~ndimage.binary_dilation(scar, iterations=15) & valid       # 300 m to 1 km outside it
OFFSET = float(np.nanmedian(dnbr_raw[ring]))
dnbr = dnbr_raw - OFFSET
print(f"unburned reference ring: {ring.sum() * HA_PER_PX:,.0f} ha, median dNBR {OFFSET:+.3f} -> subtracted as the offset")

# --- Step 4: USGS severity classes and burned area in hectares --------------------------------
CLASSES = [  # Key and Benson (2006) dNBR ranges, as published by USGS and UN-SPIDER
    ("Enhanced regrowth, high", -np.inf, -0.251, "#7a8737"),
    ("Enhanced regrowth, low", -0.251, -0.101, "#acbe4d"),
    ("Unburned", -0.101, 0.099, "#0ae042"),
    ("Low severity", 0.099, 0.269, "#fff70b"),
    ("Moderate-low severity", 0.269, 0.439, "#ffaf38"),
    ("Moderate-high severity", 0.439, 0.659, "#ff641b"),
    ("High severity", 0.659, np.inf, "#a41fd6"),
]
edges = [c[2] for c in CLASSES[:-1]]
cls = np.where(np.isfinite(dnbr), np.digitize(dnbr, edges, right=True), 255).astype("uint8")   # 0..6, 255 = masked

# Burn perimeter: low severity or worse, cleaned, in patches connected to the main scar.
burned = ndimage.binary_opening((cls >= 3) & (cls < 255), iterations=1)
lab, n = ndimage.label(burned)
keep = np.unique(lab[scar & burned])
perimeter = np.isin(lab, keep[keep > 0])
perimeter = ndimage.binary_closing(perimeter, iterations=2) & np.isfinite(dnbr)

rows = []
for k, (name, lo, hi, _) in enumerate(CLASSES):
    rng = f"< {hi:.3f}" if lo == -np.inf else (f">= {lo:.3f}" if hi == np.inf else f"{lo:.3f} to {hi:.3f}")
    rows.append({"class": name, "dNBR range": rng,
                 "whole AOI (ha)": round((cls == k).sum() * HA_PER_PX, 1),
                 "inside perimeter (ha)": round(((cls == k) & perimeter).sum() * HA_PER_PX, 1)})
area = pd.DataFrame(rows)
inside = area["inside perimeter (ha)"]
area["% of burn"] = (100 * inside / inside[3:].sum()).round(1).where(area.index >= 3, np.nan)
area.to_csv(OUT / "severity_area_ha.csv", index=False)
section("burn severity area (USGS classes)")
print(area.to_string(index=False))
burn_ha = perimeter.sum() * HA_PER_PX
print(f"burn perimeter: {burn_ha:,.0f} ha ({burn_ha / 100:.1f} km2); "
      f"moderate-low or worse: {((cls >= 4) & (cls < 255) & perimeter).sum() * HA_PER_PX:,.0f} ha")
print(f"whole AOI at low severity or worse: {(np.nan_to_num(dnbr_raw, nan=-9) > 0.099).sum() * HA_PER_PX:,.0f} ha "
      f"without the offset, {((cls >= 3) & (cls < 255)).sum() * HA_PER_PX:,.0f} ha with it")

# Save GeoTIFFs and the perimeter as polygons.
profile = dict(driver="GTiff", width=WIDTH, height=HEIGHT, count=1, crs=CRS, transform=TRANSFORM, compress="deflate")
with rasterio.open(OUT / "dnbr_2018.tif", "w", dtype="float32", nodata=np.nan, **profile) as dst:
    dst.write(dnbr.astype("float32"), 1)
with rasterio.open(OUT / "severity_classes_2018.tif", "w", dtype="uint8", nodata=255, **profile) as dst:
    dst.write(cls, 1)
    dst.write_colormap(1, {k: tuple(int(c[3][i:i + 2], 16) for i in (1, 3, 5)) + (255,) for k, c in enumerate(CLASSES)})
shapes = rasterio.features.shapes(perimeter.astype("uint8"), mask=perimeter, transform=TRANSFORM)
poly = gpd.GeoDataFrame.from_features([{"geometry": g, "properties": {"fire": "Saddleworth Moor 2018"}} for g, _ in shapes],
                                      crs=CRS).dissolve("fire").reset_index()
poly["area_ha"] = poly.area / 10_000
poly.to_file(OUT / "burn_perimeter_2018.gpkg", driver="GPKG")
print(f"perimeter polygon: {poly['area_ha'].iloc[0]:,.0f} ha, {poly.geometry.iloc[0].length / 1000:.1f} km of edge")

# --- Step 5: RdNBR, the relative version ------------------------------------------------------
pre_nbr = np.where(valid, pre["NBR"], np.nan)
with np.errstate(divide="ignore", invalid="ignore"):
    rdnbr = dnbr / np.sqrt(np.clip(np.abs(pre_nbr), 0.001, None))   # Miller and Thode (2007), unscaled units
section("RdNBR")
m = perimeter & np.isfinite(rdnbr)
print(f"pre-fire NBR inside the perimeter: median {np.median(pre_nbr[m]):.3f}, "
      f"5th-95th pct {np.percentile(pre_nbr[m], 5):.3f} to {np.percentile(pre_nbr[m], 95):.3f}")
print(f"inside perimeter: median dNBR {np.median(dnbr[m]):.3f}, median RdNBR {np.median(rdnbr[m]):.3f}")
low_pre = m & (pre_nbr < np.percentile(pre_nbr[m], 25))
high_pre = m & (pre_nbr > np.percentile(pre_nbr[m], 75))
print(f"lowest quarter of pre-fire NBR:  median dNBR {np.median(dnbr[low_pre]):.3f}, RdNBR {np.median(rdnbr[low_pre]):.3f}")
print(f"highest quarter of pre-fire NBR: median dNBR {np.median(dnbr[high_pre]):.3f}, RdNBR {np.median(rdnbr[high_pre]):.3f}")
RDNBR_EDGES = [0.069, 0.316, 0.641]   # Miller and Thode (2007): 69, 316, 641 on their x1000 scale
rd_cls = np.digitize(rdnbr[m], RDNBR_EDGES)
dn_cls = np.digitize(dnbr[m], [0.1, 0.27, 0.66])
print("hectares inside the perimeter by four-class scheme:")
for k, name in enumerate(["unchanged", "low", "moderate", "high"]):
    print(f"  {name:<9}  RdNBR {np.sum(rd_cls == k) * HA_PER_PX:6,.0f} ha | dNBR {np.sum(dn_cls == k) * HA_PER_PX:6,.0f} ha")

# --- Step 6: the fire's progression from successive images ----------------------------------------
section("fire progression")
first_seen = np.zeros((HEIGHT, WIDTH), dtype="uint8")
prog = []
zone = ndimage.binary_dilation(perimeter, iterations=5)
for k, d in enumerate(PROGRESSION, 1):
    s = load_scene(d, bands=("B8A", "B12"))
    ok = pre["clear"] & s["clear"]
    burnt_now = ok & zone & ((pre["NBR"] - s["NBR"] - OFFSET) > 0.27)
    first_seen[burnt_now & (first_seen == 0)] = k
    prog.append({"date": d, "clear_in_perimeter_%": round(100 * ok[perimeter].mean(), 1),
                 "burnt_ha": int(round(burnt_now.sum() * HA_PER_PX))})
prog = pd.DataFrame(prog)
print(prog.to_string(index=False))

# --- Step 7: recovery, 2017 to 2025 ----------------------------------------------------------------
section("recovery time series")
iy, ix = np.nonzero(perimeter)
pad = 60
r0, r1 = max(iy.min() - pad, 0), min(iy.max() + pad, HEIGHT)
c0, c1 = max(ix.min() - pad, 0), min(ix.max() + pad, WIDTH)
SUB = (x0 + c0 * RES, y1 - r1 * RES, x0 + c1 * RES, y1 - r0 * RES)
sub_burn = ((cls >= 4) & (cls < 255) & perimeter)[r0:r1, c0:c1]       # moderate-low severity or worse
sub_ctrl = ring[r0:r1, c0:c1]
sub_ll = gpd.GeoSeries.from_xy([SUB[0], SUB[2]], [SUB[1], SUB[3]], crs=CRS).to_crs(4326).total_bounds
long_items = [it for it in catalog.search(collections=["sentinel-2-l2a"], bbox=tuple(sub_ll),
                                          datetime="2017-01-01/2025-12-31",
                                          query={"s2:mgrs_tile": {"eq": TILE}, "eo:cloud_cover": {"lt": 70}}).items()
              if 5 <= it.datetime.month <= 9]


def recovery_row(item) -> dict:
    """Median NBR of the burned and control pixels in one scene, if both are at least 80% clear."""
    scl = read_band(item, "SCL", SUB)
    ok = ~np.isin(scl, BAD_SCL + [WATER_SCL])
    row = {"date": item.datetime.strftime("%Y-%m-%d"), "baseline": item.properties["s2:processing_baseline"],
           "clear_burn": float(ok[sub_burn].mean()), "clear_ctrl": float(ok[sub_ctrl].mean())}
    if row["clear_burn"] >= 0.8 and row["clear_ctrl"] >= 0.8:
        v = nbr(to_reflectance(read_band(item, "B8A", SUB), item), to_reflectance(read_band(item, "B12", SUB), item))
        row.update(nbr_burn=float(np.nanmedian(v[sub_burn & ok])), nbr_ctrl=float(np.nanmedian(v[sub_ctrl & ok])))
    return row


rec_file = CACHE / "recovery.csv"
if not rec_file.exists():
    with ThreadPoolExecutor(max_workers=8) as pool:
        pd.DataFrame(list(pool.map(recovery_row, long_items))).to_csv(rec_file, index=False)
rec = pd.read_csv(rec_file).dropna(subset=["nbr_burn"]).sort_values("date").drop_duplicates("date")
rec["date"] = pd.to_datetime(rec["date"])
rec["gap"] = rec["nbr_ctrl"] - rec["nbr_burn"]          # how far the scar is below its unburned surroundings
rec.to_csv(OUT / "recovery_nbr.csv", index=False)
print(f"{len(long_items)} May-September scenes 2017-2025, {len(rec)} with >= 80% clear burn and control pixels")
rec["period"] = np.where(rec["date"] < pd.Timestamp(PRE) + pd.Timedelta(days=1), "pre-fire",
                         rec["date"].dt.year.astype(str))
yearly = rec.groupby("period", sort=False).agg(scenes=("gap", "size"), nbr_burn=("nbr_burn", "median"),
                                               nbr_ctrl=("nbr_ctrl", "median"), gap=("gap", "median")).round(3)
print(yearly.to_string())
post_years = yearly.drop("pre-fire")
recovered = post_years[post_years["gap"] <= yearly.loc["pre-fire", "gap"]]
print(f"pre-fire gap {yearly.loc['pre-fire', 'gap']:.3f}; first post-fire summer at or below it: "
      f"{recovered.index[0] if len(recovered) else 'none yet'}")

# --- Figures -----------------------------------------------------------------------------------
EXTENT = (x0, x1, y0, y1)
SEV_CMAP = ListedColormap([c[3] for c in CLASSES])
SEV_NORM = BoundaryNorm(np.arange(-0.5, 7.5), SEV_CMAP.N)


def km_axes(ax, title=None, xlim=(x0, x1), ylim=(y0, y1)):
    ax.set_xticks([]), ax.set_yticks([])
    ax.set_xlim(*xlim), ax.set_ylim(*ylim)
    if title:
        ax.set_title(title, fontsize=9)
    xr, yb = xlim[1], ylim[0]
    ax.plot([xr - 1300, xr - 300], [yb + 300] * 2, color="black", lw=3, solid_capstyle="butt")
    ax.text(xr - 800, yb + 420, "1 km", ha="center", va="bottom", fontsize=8,
            bbox=dict(fc="white", ec="none", alpha=0.75, pad=1))


def rgb(s, bands, gain):
    img = np.dstack([s[b] for b in bands]) * gain
    return np.clip(np.nan_to_num(img), 0, 1) ** (1 / 1.4)


fig, axes = plt.subplots(1, 3, figsize=(13, 3.9), constrained_layout=True)
axes[0].imshow(rgb(pre, ["B04", "B03", "B02"], 3.5), extent=EXTENT)
km_axes(axes[0], f"True colour, {PRE} (before)")
axes[1].imshow(rgb(post, ["B04", "B03", "B02"], 3.5), extent=EXTENT)
km_axes(axes[1], f"True colour, {POST} (after)")
axes[2].imshow(rgb(post, ["B12", "B8A", "B04"], 2.2), extent=EXTENT)
km_axes(axes[2], f"False colour B12/B8A/B4, {POST}")
fig.savefig(OUT / "01-true-false-colour.png")
plt.close(fig)

fig, axes = plt.subplots(1, 3, figsize=(13, 4.0), constrained_layout=True)
for ax, arr, t in ((axes[0], np.where(pre["clear"], pre["NBR"], np.nan), f"NBR before ({PRE})"),
                   (axes[1], np.where(post["clear"], post["NBR"], np.nan), f"NBR after ({POST})")):
    im = ax.imshow(arr, cmap="BrBG", vmin=-0.4, vmax=0.6, extent=EXTENT)
    km_axes(ax, t)
fig.colorbar(im, ax=axes[:2], shrink=0.8, label="NBR")
im = axes[2].imshow(dnbr, cmap="RdYlGn_r", vmin=-0.3, vmax=0.9, extent=EXTENT)
km_axes(axes[2], f"dNBR (pre − post, offset {OFFSET:+.3f} removed)")
fig.colorbar(im, ax=axes[2], shrink=0.8, label="dNBR")
fig.savefig(OUT / "02-nbr-dnbr.png")
plt.close(fig)

fig = plt.figure(figsize=(12.5, 5.0), constrained_layout=True)
gs = fig.add_gridspec(1, 2, width_ratios=[1.7, 1])
ax = fig.add_subplot(gs[0])
ax.imshow(np.mean(rgb(post, ["B04", "B03", "B02"], 3.5), axis=2), cmap="gray", vmin=0, vmax=1, extent=EXTENT)
ax.imshow(np.ma.masked_where(~perimeter, cls), cmap=SEV_CMAP, norm=SEV_NORM, extent=EXTENT, interpolation="nearest")
ax.contour(perimeter, levels=[0.5], extent=(x0, x1, y1, y0), colors="black", linewidths=0.7)
km_axes(ax, f"Burn severity (USGS dNBR classes), Saddleworth Moor, {POST}")
ax2 = fig.add_subplot(gs[1])
sev = area.iloc[2:]
bars = ax2.barh(sev["class"].str.replace(" severity", ""), sev["inside perimeter (ha)"],
                color=[c[3] for c in CLASSES[2:]], edgecolor="black", lw=0.5)
ax2.bar_label(bars, labels=[f"{v:,.0f} ha" for v in sev["inside perimeter (ha)"]], padding=3)
ax2.invert_yaxis()
ax2.set_xlabel("Area inside the burn perimeter (ha)")
ax2.set_xlim(0, sev["inside perimeter (ha)"].max() * 1.35)
ax2.spines[["top", "right"]].set_visible(False)
ax2.set_title(f"Burn perimeter: {burn_ha:,.0f} ha", fontsize=9)
fig.savefig(OUT / "03-severity-map.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(8.5, 5.4), constrained_layout=True)
prog_cols = ["#fde725", "#f98e09", "#bc3754", "#57106e"]
ax.imshow(np.mean(rgb(pre, ["B04", "B03", "B02"], 3.5), axis=2), cmap="gray", vmin=0, vmax=1, extent=EXTENT)
ax.imshow(np.ma.masked_equal(first_seen, 0), cmap=ListedColormap(prog_cols),
          norm=BoundaryNorm(np.arange(0.5, 5.5), 4), extent=EXTENT, interpolation="nearest")
km_axes(ax, "First Sentinel-2 image in which each pixel shows as burnt (dNBR > 0.27)")
ax.legend(handles=[Patch(color=c, label=f"{pd.Timestamp(d):%d %b}: {h:,} ha burnt by then")
                   for c, d, h in zip(prog_cols, PROGRESSION, prog["burnt_ha"])],
          loc="upper left", fontsize=8, framealpha=0.9)
fig.savefig(OUT / "04-fire-progression.png")
plt.close(fig)

fig, axes = plt.subplots(1, 3, figsize=(13, 4.0), constrained_layout=True, gridspec_kw={"width_ratios": [1, 1, 0.95]})
zx, zy = (x0 + c0 * RES, x0 + c1 * RES), (y1 - r1 * RES, y1 - r0 * RES)
for ax, arr, t, vmax in ((axes[0], np.where(perimeter, dnbr, np.nan), "dNBR inside the perimeter", 0.9),
                         (axes[1], np.where(perimeter, rdnbr, np.nan), "RdNBR inside the perimeter", 1.5)):
    im = ax.imshow(arr, cmap="inferno_r", vmin=0, vmax=vmax, extent=EXTENT)
    km_axes(ax, t, zx, zy)
    fig.colorbar(im, ax=ax, shrink=0.8)
hb = axes[2].hexbin(pre_nbr[m], dnbr[m], gridsize=45, cmap="Blues", mincnt=1, bins="log")
axes[2].set_xlabel("pre-fire NBR")
axes[2].set_ylabel("dNBR")
axes[2].set_title("Burned pixels: dNBR against pre-fire NBR", fontsize=9)
fig.colorbar(hb, ax=axes[2], shrink=0.8, label="pixels (log scale)")
fig.savefig(OUT / "05-rdnbr-comparison.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(10, 4.2), constrained_layout=True)
ax.plot(rec["date"], rec["nbr_ctrl"], "o", ms=3.5, color="#1d6fa3", alpha=0.7, label="unburned moorland around the scar")
ax.plot(rec["date"], rec["nbr_burn"], "o", ms=3.5, color="#d9480f", alpha=0.7, label="burned, moderate-low severity or worse")
yr_mid = pd.to_datetime([f"{y}-08-01" for y in post_years.index])
ax.plot(yr_mid, post_years["nbr_ctrl"], "-", color="#1d6fa3", lw=1.5, label="median of each summer after the fire")
ax.plot(yr_mid, post_years["nbr_burn"], "-", color="#d9480f", lw=1.5)
for col, c in (("nbr_ctrl", "#1d6fa3"), ("nbr_burn", "#d9480f")):
    ax.hlines(yearly.loc["pre-fire", col], rec["date"].min(), pd.Timestamp(PRE), colors=c, lw=1.5, ls=":")
ax.axvline(pd.Timestamp(PRE), color="black", lw=1, ls="--")
ax.annotate("fire, 24 June 2018", (pd.Timestamp(PRE), 0.02), xytext=(6, 0), textcoords="offset points",
            xycoords=("data", "axes fraction"), fontsize=8)
ax.set_ylabel("median NBR, May–September scenes")
ax.xaxis.set_major_locator(mdates.YearLocator())
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
ax.spines[["top", "right"]].set_visible(False)
ax.legend(loc="lower right", fontsize=8)
fig.savefig(OUT / "06-recovery-timeseries.png")
plt.close(fig)

section("outputs")
for f in sorted(OUT.iterdir()):
    if f.name != "console.txt":
        print(f"{f.name:<30} {f.stat().st_size / 1024:8.0f} KB")