"""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")