From 8a2e032bbcb459bd12ebd054102746f244067c23 Mon Sep 17 00:00:00 2001 From: Scott W McCoy <106191879+scottwmccoy@users.noreply.github.com> Date: Fri, 14 Aug 2026 23:54:54 -0700 Subject: [PATCH] Add near-real-time burn severity from CIMSS BRISK New stormscape.burn module and burn subcommand: screen the BRISK archive for fires intersecting an AOI, cache and mosaic the scenes, classify severity, and drape it over the hillshade. BRISK maps every large US fire daily from a nine-satellite dNBR composite, which is what makes it usable while a fire is still burning -- BAER soil burn severity lands days to weeks after containment, MTBS a year later. The portal's WMTS serves rendered PNG capped at zoom 7, so this reads the per-fire GeoTIFFs from the open archive behind it. Finding which fires touch an AOI costs no pixels: GDAL range-reads the headers (~60 fires in 1.4 s) and the footprints are memoised. Maps draw in the BAER class palette, banded into severity classes with the class names on the colour bar, so a stormscape map can sit beside a BAER product. That palette is read from the products themselves -- all 77 of the 2025 BAER soil-burn-severity rasters embed the same four class colours, and they are the four BRISK publishes. Validated against the BAER teams' own dNBR on the 39 fires carrying both (examples/brisk_vs_baer.py): BRISK computes the same quantity -- median slope 1.013, bias +0.003, r = 0.938 over 1.24 M cells. The disagreements are compositing latency, not a different algorithm: every poor performer read low on the BAER date and recovered to r = 0.80-0.96 given a scene 5-21 days later. So scenes now carry an age, immature composites are flagged, and --min-age turns that into a filter. dNBR is a vegetation index, not soil burn severity, and the USGS debris-flow models are calibrated on the latter -- flagged in the README, the CLI help and the module docstring. Also: drape_i15 gains norm/cbar_ticks/cbar_ticklabels for classed maps; export learns that severity is categorical and scales burn fields on their own range; brisk_cache/ joins the reserved event-root directories. 66 offline tests (398 total). Co-Authored-By: Claude Opus 5 --- .gitignore | 1 + CLAUDE.md | 144 ++++++- README.md | 87 +++- docs/cli.md | 62 ++- examples/brisk_vs_baer.py | 224 +++++++++++ stormscape/__init__.py | 5 + stormscape/burn.py | 713 +++++++++++++++++++++++++++++++++ stormscape/cli.py | 258 ++++++++++++ stormscape/export.py | 12 +- stormscape/layout.py | 2 +- stormscape/plot.py | 29 +- tests/test_burn.py | 808 ++++++++++++++++++++++++++++++++++++++ 12 files changed, 2332 insertions(+), 13 deletions(-) create mode 100644 examples/brisk_vs_baer.py create mode 100644 stormscape/burn.py create mode 100644 tests/test_burn.py diff --git a/.gitignore b/.gitignore index 22aa218..fa74e3c 100644 --- a/.gitignore +++ b/.gitignore @@ -34,3 +34,4 @@ RainGaugeData/ VirtualGaugeFigures/ # Sample figures can still be committed with `git add -f *.png` +baer_study/ diff --git a/CLAUDE.md b/CLAUDE.md index 78c586d..27de13e 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -15,7 +15,10 @@ vector file, or a shapely geometry) plus a **storm-day date** into: labelled vector overlays (NHD streams, TIGER roads, GNIS places) and gauge points, plus a **radar-vs-gauge comparison** — per-gauge residuals + skill stats (`stormscape/plot.py`, `stormscape/compare.py`, `stormscape/refdata.py`); -5. raw **single-radar NEXRAD Level II** tilts (reflectivity/velocity) for the +5. **near-real-time burn severity** — CIMSS **BRISK** daily multi-satellite + dNBR composites (+ the BAER soil-burn-severity archive) over the AOI, cached, + classified and mapped (`stormscape/burn.py`); +6. raw **single-radar NEXRAD Level II** tilts (reflectivity/velocity) for the radar nearest the AOI — gridded over the AOI or sampled at the gauges, with a Z–R diagnostic (`stormscape/nexrad.py`); the underlying radar behind the MRMS mosaic, and the way to reach pre-2020 events. @@ -92,10 +95,11 @@ in `README.md`. - `stormscape/atlas14.py` — NOAA Atlas 14 gridded climatology → intensity fields + anomaly - `stormscape/smoothing.py` — NaN-aware field smoothing (gaussian/uniform/median/idw) + radar-gauge skill sweep - `stormscape/export.py` — georeferenced exports for GIS/CalTopo: EPSG:3857 GeoTIFFs (raw float + colorized RGBA) + GeoPDF figures +- `stormscape/burn.py` — near-real-time burn severity (BRISK dNBR / BAER SBS) - `stormscape/refdata.py` — AOI-scoped NHD streams / TIGER roads / GNIS places - `stormscape/plot.py` — drape i15 over hillshade + basemap/vector/gauge overlays - `stormscape/data/` — bundled tables (`nexrad_sites.csv` NCEI HOMR, `atlas14_regions.csv`) -- `stormscape/cli.py` — `stormscape {dem,i15,map,run,gauges,compare,nexrad,panels,vgauge,zoom,pick,climate,smooth,recurrence,export}` +- `stormscape/cli.py` — `stormscape {dem,i15,map,run,gauges,compare,nexrad,panels,vgauge,zoom,pick,climate,smooth,recurrence,export,burn}` - `environment.yml`, `pyproject.toml`, `examples/` ## Conventions & gotchas (learned during development) @@ -906,6 +910,142 @@ in `README.md`. **Gotcha for anyone testing this:** RQI/SHSR are read with `fetch` (singular), not `fetch_many` — mock both or the "offline" test hits S3 and CI's `-m "not network"` run breaks. +- **Near-real-time burn severity (`burn.py`, CLI `burn`, 2026-08-14) — CIMSS + BRISK.** Puts a scar under the rain *while the fire is still burning*, which the + authoritative products cannot: BAER soil burn severity lands days-to-weeks after + containment and only for assessed fires, MTBS a year+. BRISK is a Google Earth + Engine **dNBR data-fusion composite over nine satellites** (GOES-E/W ABI, SNPP / + NOAA-20 / NOAA-21 VIIRS, Landsat 8/9, Sentinel-2a/b), mapping every large + (>~5,000 acre) US fire **daily**. + **Access — the portal is a decoy.** `cimss.ssec.wisc.edu/brisk` is a RealEarth + viewer whose WMTS (`re-brisk.ssec.wisc.edu/wmts/BurnScars-dNBR.xml`) serves + **rendered PNG tiles, capped at zoom 7** — pretty, not data. The raw field is in + the open Apache-indexed archive behind it, **one GeoTIFF per fire per day**: + `bin.ssec.wisc.edu/pub/realearth/brisk//--dNBR__235959.tif` + (2025→present; 8k+ scenes / ~590 fires in 2026 alone, current through *today*). + So `burn.py` scrapes the directory index into a cached catalog and downloads + only the scenes that intersect the AOI. + **AOI screening is header-only.** GDAL range-reads a remote GeoTIFF header via + `/vsicurl/` in ~0.02 s across a 12-thread pool — a whole day's ~60 fires screen + in **1.4 s without downloading a pixel**; footprints are memoised in + `brisk_cache/bounds_dnbr.json`, so the first full-archive screen (~650 scenes, + ~25 s) makes every later one **0.04 s**. Index cache TTL is **6 h for the + current year and infinite for past years** (closed years never gain scenes), + and a dead server falls back to the stale index instead of failing. + **Gotchas, all measured:** + (a) scenes are **EPSG:3857 at "60 m", which is not 60 m of ground** — Web + Mercator metres shrink with latitude, so the cell is **46.5 m at 39°N**; areas + computed on the raw transform run ~66% high there, hence the `cos(lat)` + correction in `cli._burn_class_table`; + (b) every scene tested lands on an **exact 60 m multiple**, so mosaicking + neighbouring fires is a paste, not a warp; + (c) **NaN marks outside-the-burn but the files tag no nodata**, so + `masked=True` masks nothing — test `np.isfinite` (typically only 10-19% of a + scene is valid); + (d) dNBR is **unscaled** (~-0.3 to 1.0), not the ×1000 integer form the + MTBS/USGS thresholds are quoted in — `SEVERITY_SCHEMES["usgs"]` is + 0.10/0.27/0.44/0.66, `["brisk"]` the portal's own 0.10/0.40/0.70; + (e) a fire's footprint **grows day to day**, so `find_scenes` takes the latest + scene per fire **on or before `--date`** — a storm-day map must not include + severity mapped after the storm. + **Mosaic rule: NaN-aware `np.fmax`.** Scenes are mostly NaN outside their own + fire, so a plain "first/last wins" merge blanks whichever fire is written second + wherever footprints overlap; `fmax` combines, is order-independent, and takes + the more severe value where two fires genuinely overlap. Each scene is + `reproject`d into one destination grid, so the same code path handles the + BAER scenes (which are **per-fire UTM at 20 m**, not 3857). + **BAER SBS (`--product sbs`, `baer-data/`)** is the authoritative soil product + but **sparse** (8 fires in 2026 vs BRISK's ~590) and classified uint8. Its + embedded palette colours 1-4 with BRISK's own four severity colours and paints + **0 and 5+ the same black** — i.e. the product itself treats anything outside + 1-4 as a mask (water/inholding/unmapped, ~3% of the NV scene checked), so + `SBS_VALID = (1, 4)` reads the rest as missing rather than charting it as a + class. Categorical mosaics use **nearest** (averaging class 1 and 3 into 2 + would be a fabricated severity). + **Display — the BAER palette is the default, and it IS the BRISK palette + (2026-08-14).** Measured, not assumed: **all 77** of the 2025 BAER + soil-burn-severity rasters — written by many different BAER teams, in ERDAS + Imagine — embed an *identical* class palette, `1 (0,128,128) teal / 2 + (82,204,204) cyan / 3 (255,232,32) yellow / 4 (168,0,0) dark red` (only the + class-5 mask colour varies: black/gray/white). Those are exactly the four + colours in BRISK's `qgis_BRISK_dNBR_colorscale_v2.txt`, so the portal and the + BAER deliverables already share one scheme — `BAER_CLASS_COLORS` / + `BAER_ANCHORS` (`BRISK_ANCHORS` is an alias; `register_brisk_cmap` → + `register_baer_cmap`). Maps are **classed by default**, the way BAER + publishes: `burn.severity_colors(scheme)` returns + `(ListedColormap, BoundaryNorm, ticks, labels)` so the field is banded at the + severity breaks and the **colour bar is labelled with class names, not dNBR + numbers**. A 4-class scheme uses the official table exactly; the 5-class + `usgs` scheme has no official 5th colour so its colours are *sampled from the + same ramp* at class midpoints (documented, not smuggled). `--continuous` gives + the smooth ramp (same colours), an explicit `--cmap` opts out entirely. + This required `plot.drape_i15` to gain **`norm=` / `cbar_ticks=` / + `cbar_ticklabels=`** (a norm replaces the linear `vmin=0,vmax=` scale) — small + and general, usable by any future classed map. `_burn_display_defaults` + resolves the cut/scale **per product** — dNBR 0.10/1.0, SBS 1.5/4.0 — because + one shared default would clip the class map to its lowest class or paint + unburned ground. The drape reuses `drape_i15` (so `--reference`/`--clip`/north + arrow/ticks all just work); `--alpha` is worth raising from the project 0.32 + when the scar is the subject. + **BAER also publishes its own dNBR, and BRISK matches it (validated + 2026-08-14).** `baer-data//` holds `--prelim-dNBR__*.tif` + (**118 for 2025**, none yet for 2026) alongside the `-sbs_` rasters — note the + **`-prelim-` infix**, which `parse_name` strips. They are **int16 ×1000** + (the BARC convention; NOAA's own `BARC256 = dNBR×5 − 275` identity confirms the + scaling), **ESRI:102039** Albers, **20 m where the source was Sentinel-2 and + 30 m where it was Landsat**. **39 fires have both products.** Method: divide + BAER by 1000 and `Resampling.average` it **down onto the BRISK grid** — + aggregating the finer product to the coarser support rather than upsampling + the thing under test — then score only **burned cells (BAER dNBR ≥ 0.1)**, + because the huge unburned surround agreeing near zero inflates r by ~0.03–0.1. + **Result: BRISK computes the same dNBR.** Pooled over 8 well-matched fires, + **r = 0.938 on 1.24 M cells, slope 1.089**; across all 39, median **slope + 1.013, bias +0.003, ratio 1.011** — unity slope, no bias. **The disagreements + are compositing latency, not a different algorithm**: every poor performer + (middle-mesa-nm r 0.03, turkeyfeather-nm 0.12, blind-az 0.22, derby-co 0.27, + dillon-ca 0.58, laguna-nm 0.60, island-creek-id 0.61) recovers to **r 0.80–0.96 + when given a BRISK scene 5–21 days later**, and every one of them read *low* on + the BAER date — the composite had not yet ingested a clear post-fire overpass. + A pre-set (not per-fire-tuned) **+14 d** rule lifts the distribution: median r + 0.874→0.913, r≥0.90 **14→24** of 39, r<0.60 **5→1**, IQR 0.765–0.934 → 0.866–0.932. + Note it rescues the *tail* rather than improving the median fire (21 better, + 18 slightly worse) — a mature scar keeps darkening away from the BAER snapshot, + which is why the +14 d slope rises to ~1.05–1.09. **BAER 20 m (Sentinel-2) + agrees better than 30 m (Landsat)** — median r 0.923 vs 0.814, Mann-Whitney + **p = 0.0005** — most likely because BRISK's composite and the Sentinel-based + BAER share the same underlying acquisition. **Operational rule: trust BRISK's + *pattern* immediately, but give the composite ~2 weeks before trusting its + *magnitude*** — early scenes under-read. Study script + CSVs live in the + repo as **`examples/brisk_vs_baer.py`** (reruns the whole study from the two + catalogs; writes CSV + figure to `--out-dir`, default `baer_study/`, which is + gitignored). The BAER dNBR is reachable as **`--product baer_dnbr`** + (`PRODUCTS[...]["scale"] = 1000.0` divides it back to a plain index on read). + **Acted on in the tool:** `burn.MATURITY_DAYS = 14`; `find_scenes` adds an + **`age_days`** column (days since the fire's *first* appearance in the + catalog, computed **before** any `--date`/`--since` trim, or trimming early + scenes would make an old fire look new), `--list` prints it and stars + immature composites, `burn_severity` prints an advisory, and **`--min-age + DAYS`** hard-filters. The maturity screen runs **after** the AOI + intersection -- filtering the national catalog first named 593 irrelevant + fires in the drop message. Off by default: an immature scar still beats no + scar, provided the run says so. + **`dnbr`/`severity` are deliberately NOT in `plot._MASK_DRY`** — that mask cuts + below 0.5, which on a 0-1 dNBR field would erase everything short of high + severity, and it also drives smoothing, which a severity field does not want. + **THE CAVEAT THAT MATTERS FOR THIS GROUP'S SCIENCE: dNBR is a *vegetation* + index, not soil burn severity.** The USGS post-fire debris-flow models are + calibrated on **soil** burn severity — dNBR adjusted by field crews for + hydrophobicity, ground cover and duff consumption. BRISK is explicitly + **interim**: act on it early, then supersede with `--product sbs`. + Outputs are `mrms`-style result dicts (`save_fields` / `drape_i15` just work): + `_dnbr.tif`, `_severity.tif`, `_burn_classes.csv` (pixels + true + km² + fraction per class), `_burn_scenes.geojson` (provenance: which fire, + which date), `_burn.png`. Cache is `brisk_cache/` (added to + `layout.RESERVED`, matches `nexrad_cache/`; keep out of git). **No new deps.** + Validated live on Ward NV 20260814 (91.2 km² burned, 6.3% high severity), + a 15-fire central-Oregon mosaic, and Cottonwood-Peak NV SBS. + Tests: `tests/test_burn.py` (48, offline — synthetic listings + local GeoTIFFs). + - **Testing (`tests/`, pytest, added 2026-07-29) — run it before every push.** `pytest` (or `/opt/anaconda3/envs/GISMan/bin/python -m pytest`) — **232 tests, ~2 s, entirely offline** (no MRMS/NEXRAD/Synoptic/3DEP/Atlas 14 request, no diff --git a/README.md b/README.md index 4f621c1..edc3ef8 100644 --- a/README.md +++ b/README.md @@ -25,7 +25,10 @@ shapely geometry) and, for rainfall, a **storm-day date**, it will: 5. **Reach the raw single-radar NEXRAD Level II volumes** for the radar nearest the AOI — reflectivity / velocity at each elevation tilt, gridded over the AOI or sampled at the gauges (the underlying radar behind the MRMS mosaic; - `stormscape.nexrad`). + `stormscape.nexrad`); +6. **Map near-real-time burn severity** over the same AOI from **CIMSS BRISK** + daily multi-satellite dNBR — so a scar can be put under the rain while the + fire is still burning (`stormscape.burn`). It was lifted out of a post-fire debris-flow study and generalized: nothing here is tied to that project's fire perimeters or event inventory — you supply @@ -552,6 +555,77 @@ picks the format; `--streams-bbox` keeps whole flowlines over the bounding box; --- +### Near-real-time burn severity (BRISK dNBR) + +`burn` maps the **burn scar** over an AOI. It reads **CIMSS BRISK** — a daily, +nine-satellite dNBR composite (GOES ABI, VIIRS, Landsat 8/9, Sentinel-2) that +covers every large US fire — so severity is available *while the fire is still +burning*, which the authoritative products are not: BAER soil burn severity +arrives days to weeks after containment and only for fires that get an +assessment, MTBS a year or more later. + +```bash +# what has burned in this AOI? (headers only — downloads nothing) +python -m stormscape burn --aoi event_AOI.kmz --list +``` + +```bash +# fetch, cache, classify and map it +python -m stormscape burn --aoi event_AOI.kmz --date 20260814 \ + --out-dir ./out --key event --dem --reference --clip +``` + +That writes `event_dnbr.tif` (the composite index), `event_severity.tif` (the +class index), `event_burn_classes.csv` (pixels, **true ground area** and fraction +per class), `event_burn_scenes.geojson` (which fire each scene came from, and +when), and `event_burn.png` draped over the hillshade. + +`--date` asks for the scar **as of** that date — each fire is shown at its latest +scene on or before it, so a storm-day map does not include severity mapped after +the storm. `--since` drops older seasons, `--fire NAME` restricts to one incident, +and `--product sbs` switches to the BAER **soil** burn severity archive where an +assessment exists. `--scheme` picks the dNBR breaks: `usgs` (MTBS/USGS +0.10/0.27/0.44/0.66, the default) or `brisk` (0.10/0.40/0.70). + +Maps are drawn in the **BAER class palette** and **banded into severity classes**, +the way burn severity is published — the colour bar carries class names rather than +dNBR numbers, so a stormscape map can be laid beside a BAER product. (That palette +is verified from the products themselves: all 77 of the 2025 BAER soil-burn-severity +rasters embed the same four class colours, and they are the same four BRISK +publishes.) `--continuous` shades a smooth ramp in the same colours instead, and +`--alpha` makes the scar read more strongly over the terrain. + +**Does BRISK's dNBR match the BAER teams' own?** Yes. Of the fires with both +products in 2025 (39 of them), agreement over burned cells has median **slope +1.013, bias +0.003, ratio 1.011** — unity slope, no bias — and **r = 0.938 across +1.24 M cells** pooled over well-matched fires. The exceptions are **compositing +latency, not a different algorithm**: every poor performer read *low* on the BAER +date and recovered to r = 0.80–0.96 given a scene 5–21 days later. A fixed +14-day +rule takes fires at r ≥ 0.90 from 14/39 to 24/39 and fires below 0.60 from 5 to 1. +So **trust the pattern immediately, but give the composite about two weeks before +trusting the magnitude.** BAER products derived from Sentinel-2 (20 m) agree better +than Landsat-derived ones (30 m) — median r 0.923 vs 0.814, p = 0.0005. + +`burn` acts on that: every run reports each scene's **age** (days since the fire +entered the archive) and flags composites younger than 14 days, and `--min-age DAYS` +turns that into a hard filter, skipping too-young fires by name. It is off by +default — an immature scar still beats no scar, as long as the run says so. +`--product baer_dnbr` fetches the BAER teams' own dNBR for comparison, and +[`examples/brisk_vs_baer.py`](examples/brisk_vs_baer.py) reproduces the whole study. + +Scenes are cached in `brisk_cache/` and the archive index in the same place, so +re-running a map costs no network at all. The first AOI screen reads a few hundred +GeoTIFF *headers* over HTTP range requests (~20 s) and memoises the footprints; +after that it is instant. + +> **dNBR is a *vegetation* index, not soil burn severity.** The USGS post-fire +> debris-flow models are calibrated on **soil** burn severity, which is dNBR +> adjusted by field crews for hydrophobicity, ground cover and duff consumption. +> BRISK is explicitly **interim**: use it to act early, then supersede it with +> `--product sbs` when the BAER assessment lands. + +--- + ## Python API ```python @@ -658,6 +732,11 @@ kept, and 2-min `PrecipRate` is stacked over each contiguous wet run. - **NOAA Atlas 14** precipitation-frequency climatology via the gridded ASCII grids on the PFDS GIS server (`hdsc.nws.noaa.gov/pub/hdsc/data/...`); region extents bundled in `data/atlas14_regions.csv`. +- **Burn severity** from **CIMSS BRISK** (UW-Madison SSEC) — daily per-fire dNBR + composites in the open archive behind the + [BRISK portal](https://cimss.ssec.wisc.edu/brisk/) + (`bin.ssec.wisc.edu/pub/realearth/brisk/`), plus the companion **BAER soil burn + severity** archive (`.../baer-data/`). ### DEM resampling @@ -721,10 +800,11 @@ stormscape/ │ ├── atlas14.py NOAA Atlas 14 climatology grids -> intensity fields + anomaly │ ├── smoothing.py NaN-aware field smoothing + radar-gauge skill sweep │ ├── export.py EPSG:3857 GeoTIFFs, GeoPDF figures, NHD stream vectors +│ ├── burn.py near-real-time burn severity (CIMSS BRISK dNBR / BAER SBS) │ ├── refdata.py AOI-scoped NHD streams / TIGER roads / GNIS places │ ├── plot.py drape fields over hillshade + basemap/vector/gauge overlays │ ├── data/ bundled tables (nexrad_sites.csv, atlas14_regions.csv) -│ └── cli.py 15 subcommands (see docs/cli.md) +│ └── cli.py 16 subcommands (see docs/cli.md) ├── examples/ worked example + batch event templates ├── docs/cli.md full CLI reference ├── environment.yml conda-forge environment ("stormscape") @@ -753,7 +833,8 @@ science math against known answers rather than just importing things: | comparison | bias/RMSE/ratio arithmetic; the RQI and cadence screens; cadence must **not** filter the storm total | | merging | mean-field factor is Σgauge/Σradar; leave-one-out beats the raw field on bias | | export | EPSG:3857 output; dry and no-data cells are transparent; categorical fields use nearest-neighbour; the GeoPDF neatline bounds the map frame | -| CLI | all 15 subcommands build; documented flags still exist on every command | +| burn severity | dNBR breaks match the published x1000 thresholds; unburned stays distinct from unobserved; a scene's NaN surround never erases a neighbouring fire; class mosaics are never interpolated; a past-year archive index caches forever while today's expires | +| CLI | all 16 subcommands build; documented flags still exist on every command | Tests needing an optional dependency (Py-ART, wradlib, GDAL's PDF driver) are marked `optional_deps` and skip themselves. Anything hitting the network would be diff --git a/docs/cli.md b/docs/cli.md index 206905f..513a248 100644 --- a/docs/cli.md +++ b/docs/cli.md @@ -5,7 +5,8 @@ package is installed). This page mirrors `python -m stormscape --help` run that for the authoritative, always-current text. ``` -stormscape {dem, i15, map, run, gauges, compare, nexrad, panels, vgauge, zoom, pick} +stormscape {dem, i15, map, run, gauges, compare, nexrad, panels, vgauge, zoom, + pick, climate, smooth, recurrence, export, burn} ``` | Command | Purpose | @@ -24,6 +25,8 @@ stormscape {dem, i15, map, run, gauges, compare, nexrad, panels, vgauge, zoom, p | `climate` | NOAA Atlas 14 rainfall climatology vs observed I15/I30/I60 — comparison figure + anomaly maps (reuses rasters; no re-run) | | `smooth` | Smooth a radar field: methods×radii comparison figure + optional radar–gauge skill sweep (reuses rasters; no re-run) | | `recurrence` | Wet-gauge table: peak I15/I30/I60 + time-of-peak, anomaly + recurrence interval vs NOAA Atlas 14 (PFDS point query) | +| `export` | Georeferenced EPSG:3857 GeoTIFFs + GeoPDF figures + NHD stream vectors for GIS / CalTopo (reuses rasters; no re-run) | +| `burn` | Near-real-time burn severity over an AOI from CIMSS BRISK daily dNBR (or the BAER soil-burn-severity archive) | The AOI is **always** given as either `--bbox W S E N` (lon/lat degrees) or `--aoi ` (GeoJSON/SHP/GPKG/KMZ). Gauge steps need a free Synoptic @@ -496,6 +499,62 @@ python -m stormscape export --from-dir ./out --from-key HiddenValley \ --- +### `burn` — near-real-time burn severity (CIMSS BRISK dNBR) + +Maps the **burn scar** over an AOI from **CIMSS BRISK**: a daily, nine-satellite +dNBR composite (GOES ABI, VIIRS, Landsat 8/9, Sentinel-2) covering every large US +fire, so severity is available *while the fire is still burning*. Scenes are +per-fire, per-day GeoTIFFs in the open SSEC archive; this command screens the +archive for fires intersecting the AOI (reading only GeoTIFF *headers* over HTTP +range requests — it downloads nothing until it knows what it needs), caches the +scenes it does need, mosaics them, classifies severity and draws the map. + +> **dNBR is a *vegetation* index, not soil burn severity.** The USGS post-fire +> debris-flow models are calibrated on **soil** burn severity — dNBR adjusted by +> field crews for hydrophobicity, ground cover and duff consumption. BRISK is +> explicitly **interim**: act on it early, then supersede it with `--product sbs` +> when the BAER assessment lands. + +Groups **A + C**, plus: +| Arg | Type / default | Meaning | +|---|---|---| +| `--date YYYYMMDD` | str | the scar **as of** this date — each fire is shown at its latest scene on or before it (default: latest available) | +| `--since YYYYMMDD` | str | ignore scenes older than this (e.g. to drop last season's fires) | +| `--product {dnbr,sbs,baer_dnbr}` | `dnbr` | `dnbr` = BRISK daily composite; `sbs` = BAER **soil** burn severity, authoritative but only for assessed fires; `baer_dnbr` = the BAER teams' own dNBR (2025 only so far) | +| `--min-age DAYS` | float | require the composite to be at least DAYS old (since the fire entered the archive) and skip fires that are not, naming them. Off by default — an immature scar still beats none, and the run says so | +| `--scheme {usgs,brisk}` | `usgs` | dNBR severity breaks: USGS/MTBS `0.10/0.27/0.44/0.66`, or the portal's `0.10/0.40/0.70` | +| `--fire NAME…` | strs | restrict to these fires (names as printed by `--list`) | +| `--years Y…` | ints | archive years to search (default: this year and last) | +| `--list` | flag | list the fires intersecting the AOI and stop — **downloads nothing** | +| `--all-dates` | flag | with `--list`, show every scene per fire instead of only the latest | +| `--cache-dir` | path | scene + index cache (default `/brisk_cache`) | +| `--hillshade` | path | hillshade for the map backdrop (default `/_hillshade.tif`) | +| `--dem` | flag | fetch a DEM + hillshade for the AOI if none is found | +| `--resolution` | int, `10` | DEM resolution (m) for `--dem` | +| `--vmax` | float | colour-scale max (default `1.0` for dNBR, `4.0` for the `sbs` class field) | +| `--no-map` | flag | write the rasters + table only, no figure | +| `--workers` | int, `12` | parallel header reads / downloads | + +| `--continuous` | flag | shade dNBR as a smooth ramp instead of banding it into severity classes (same colours) | + +Maps default to the **BAER class palette**, banded into the severity classes BAER +publishes, with the colour bar labelled by class name rather than dNBR value — so a +stormscape map can be laid beside a BAER product. `--continuous` gives the smooth +ramp in the same colours, an explicit `--cmap` (e.g. `YlOrRd`) opts out of the class +colours entirely, and `--alpha` makes the scar read more strongly than the +project-wide 0.32 drape. + +```bash +# what has burned here? (headers only, no downloads) +python -m stormscape burn --aoi event_AOI.kmz --list + +# fetch, cache, classify and map it +python -m stormscape burn --aoi event_AOI.kmz --date 20260814 \ + --out-dir ./out --key event --dem --reference --clip +``` + +--- + ## Outputs at a glance | Command | Writes | @@ -508,6 +567,7 @@ python -m stormscape export --from-dir ./out --from-key HiddenValley \ | `nexrad` | `_refl.tif` + `_nexrad.png` (or `_i15max.tif` … with `--intensity`); volumes cached under `--cache-dir` | | `panels` | `_panels.png` | | `vgauge` | per-gauge CSVs in `RainGaugeData/`; `_vg_atlas.png` (`--atlas`); per-gauge figures in `VirtualGaugeFigures/` (`--detail`) | +| `burn` | `_dnbr.tif`, `_severity.tif`, `_burn_classes.csv` (pixels + true km² + fraction per class), `_burn_scenes.geojson` (which fire, which date), `_burn.png`; scenes cached under `brisk_cache/` | | `zoom` | `.png` (zoomed map) + `_panels.png` + (unless `--no-climate`) `_climate_compare.png`, `_clim_i{d}.tif`, `_anom_i{d}.png`/`.tif`; `--crop-rasters` adds cropped `_.tif`; `--refine-dem` adds `_dem/_hillshade.tif` | | `pick` | `_pick.html` (self-contained browser bbox picker) | | `climate` | `_clim_i{15,30,60}.tif`, `_climate_compare.png`, `_anom_i{15,30,60}.tif` + `.png`; grids cached under `atlas14_cache/` | diff --git a/examples/brisk_vs_baer.py b/examples/brisk_vs_baer.py new file mode 100644 index 0000000..fea515c --- /dev/null +++ b/examples/brisk_vs_baer.py @@ -0,0 +1,224 @@ +"""Validation study: does BRISK compute dNBR the same way the BAER teams do? + +BRISK gives a burn scar within a day; the BAER teams' own dNBR arrives days to +weeks later but is the product post-fire hazard work is built on. Both are +published in the same SSEC archive, so for every fire that has both we can put +them cell-for-cell on one grid and check. + +Run it:: + + python examples/brisk_vs_baer.py # all matched fires + python examples/brisk_vs_baer.py --years 2025 --out-dir ./baer_study + +It downloads nothing permanently -- both products are read straight over HTTP +with GDAL range requests -- and writes a per-fire CSV plus a two-panel figure. + +Method +------ +Four choices matter, and each is a place a careless comparison goes wrong: + +1. **Scaling.** BAER dNBR is int16 **x1000** (the BARC convention; NOAA's own + ``BARC256 = dNBR*5 - 275`` identity confirms it). BRISK is float32 unscaled. + Divide BAER by 1000 or the correlation is fine and every ratio is 1000x off. + +2. **Support.** BAER is 20 m (Sentinel-2 source) or 30 m (Landsat) in Albers; + BRISK is 60 m Web-Mercator (~46 m of ground). BAER is **averaged down onto + the BRISK grid**, aggregating the finer product to the coarser support -- + upsampling BRISK instead would invent detail in the thing under test. + +3. **Which cells.** Scoring every valid cell inflates r by ~0.03-0.10, because + the large unburned surround agrees near zero and does the work. Only + **burned cells (BAER dNBR >= 0.1)** answer the question. + +4. **Which date.** BRISK is a *running composite*, so it is scored twice: on the + BAER assessment date, and at a fixed **+14 days**. The offset is set in + advance rather than tuned per fire -- picking each fire's best-correlating + scene would fit noise and inflate the answer. + +What it found (39 fires, 2025) +------------------------------ +BRISK computes the same quantity: median **slope 1.013, bias +0.003, ratio +1.011** over burned cells, and **r = 0.938 across 1.24 M cells** pooled over +well-matched fires. The poor performers are **compositing latency, not a +different algorithm** -- every one read *low* on the BAER date and recovered to +r = 0.80-0.96 given a scene 5-21 days later. The +14 d rule takes fires at +r >= 0.90 from 14/39 to 24/39 and fires below 0.60 from 5 to 1. + +So: trust BRISK's *pattern* immediately, but give the composite about two weeks +before trusting its *magnitude* -- which is what ``burn --min-age`` enforces. +""" + +from __future__ import annotations + +import argparse +import datetime as dt +import os +import warnings +from concurrent.futures import ThreadPoolExecutor + +import numpy as np +import pandas as pd +import rasterio +from rasterio.enums import Resampling +from rasterio.warp import reproject + +from stormscape.burn import catalog + +warnings.filterwarnings("ignore") + +MATURED_OFFSET_DAYS = 14 # fixed in advance; see the module docstring +BURNED_MIN_DNBR = 0.10 # BAER's own unburned break +MIN_CELLS = 100 +GDAL_ENV = dict(GDAL_DISABLE_READDIR_ON_OPEN="EMPTY_DIR", + CPL_VSIL_CURL_ALLOWED_EXTENSIONS=".tif", + GDAL_HTTP_MAX_RETRY="2", GDAL_HTTP_RETRY_DELAY="1") + + +def _matched_fires(years): + """``{fire: (baer_rows, brisk_rows)}`` for fires carrying both products.""" + baer = catalog("baer_dnbr", years=years, verbose=False) + brisk = catalog("dnbr", years=years, verbose=False) + if not len(baer) or not len(brisk): + return {} + b = {f: g for f, g in baer.groupby("fire")} + k = {f: g for f, g in brisk.groupby("fire")} + return {f: (b[f], k[f]) for f in sorted(set(b) & set(k))} + + +def _read_pair(baer_url, brisk_url): + """BAER dNBR resampled onto the BRISK grid; returns ``(brisk, baer)`` arrays.""" + with rasterio.Env(**GDAL_ENV): + with rasterio.open(brisk_url) as kds: + k = kds.read(1).astype("float64") + dst_t, dst_crs = kds.transform, kds.crs + with rasterio.open(baer_url) as bds: + src = bds.read(1).astype("float64") + if bds.nodata is not None and np.isfinite(bds.nodata): + src = np.where(src == bds.nodata, np.nan, src) + src = src / 1000.0 # x1000 -> plain dNBR + b = np.full(k.shape, np.nan) + reproject(source=src, destination=b, + src_transform=bds.transform, src_crs=bds.crs, + dst_transform=dst_t, dst_crs=dst_crs, + src_nodata=np.nan, dst_nodata=np.nan, + resampling=Resampling.average) # fine -> coarse support + res_m = abs(bds.res[0]) + return k, b, res_m + + +def _stats(brisk, baer): + m = np.isfinite(brisk) & np.isfinite(baer) & (baer >= BURNED_MIN_DNBR) + if m.sum() < MIN_CELLS: + return None + x, y = brisk[m], baer[m] + if np.std(x) == 0 or np.std(y) == 0: + return None + slope, intercept = np.polyfit(y, x, 1) # BRISK = f(BAER) + d = x - y + return dict(n=int(m.sum()), r=float(np.corrcoef(x, y)[0, 1]), + slope=float(slope), intercept=float(intercept), + bias=float(d.mean()), rmse=float(np.sqrt((d ** 2).mean())), + ratio=float(np.median(x / y))) + + +def compare_fire(fire, baer_rows, brisk_rows): + """Score one fire on the BAER date and at +14 days.""" + br = baer_rows.sort_values("date").iloc[-1] # latest BAER assessment + cand = brisk_rows.sort_values("date") + out = dict(fire=fire, baer_date=str(br.date)) + for tag, off in (("same", 0), ("d14", MATURED_OFFSET_DAYS)): + want = br.date + dt.timedelta(days=off) + pick = cand.iloc[int(np.argmin([abs((d - want).days) for d in cand.date]))] + try: + k, b, res_m = _read_pair(br.url, pick.url) + except Exception as exc: # noqa: BLE001 + print(f" !! {fire} ({tag}): {exc}", flush=True) + continue + out[f"gap_{tag}"] = (pick.date - br.date).days + out.setdefault("baer_res_m", res_m) + st = _stats(k, b) + if st: + out.update({f"{key}_{tag}": v for key, v in st.items()}) + return out + + +def summarise(df): + for tag, name in (("same", "BRISK on the BAER date"), + ("d14", f"BRISK +{MATURED_OFFSET_DAYS} days (matured)")): + r = df.get(f"r_{tag}", pd.Series(dtype=float)).dropna() + if not len(r): + continue + print(f"\n{name}: n={len(r)} fires") + print(f" median r {r.median():.3f} " + f"IQR {r.quantile(.25):.3f}..{r.quantile(.75):.3f}") + print(f" median slope {df[f'slope_{tag}'].median():.3f} " + f"bias {df[f'bias_{tag}'].median():+.4f} " + f"ratio {df[f'ratio_{tag}'].median():.3f}") + print(f" r>=0.90: {(r >= .9).sum()}/{len(r)} " + f"r>=0.80: {(r >= .8).sum()}/{len(r)} r<0.60: {(r < .6).sum()}/{len(r)}") + if "baer_res_m" in df: + print("\nby BAER source resolution (20 m = Sentinel-2, 30 m = Landsat):") + for res, g in df.groupby("baer_res_m"): + rr = g["r_same"].dropna() + if len(rr): + print(f" {int(res)} m n={len(g):2d} median r {rr.median():.3f}") + + +def plot(df, out_png): + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + d = df.dropna(subset=["r_same"]).sort_values("r_same") + fig, ax = plt.subplots(figsize=(7.5, 0.22 * len(d) + 2.2)) + y = np.arange(len(d)) + ax.hlines(y, d.r_same, d.r_d14, color="0.75", lw=1.2, zorder=1) + ax.scatter(d.r_same, y, s=26, color="#d62728", + label="BRISK on the BAER date", zorder=2) + ax.scatter(d.r_d14, y, s=26, color="#2ca02c", + label=f"BRISK +{MATURED_OFFSET_DAYS} days", zorder=2) + ax.set_yticks(y) + ax.set_yticklabels(d.fire, fontsize=6.5) + ax.axvline(0.9, color="0.5", lw=0.7, ls=":") + ax.set_xlabel("correlation with BAER dNBR (burned cells)") + ax.set_title("BRISK vs BAER dNBR per fire:\nthe failures are compositing " + "latency, not a different algorithm", fontsize=11) + ax.legend(loc="upper left", fontsize=9, framealpha=0.95) + ax.set_xlim(-0.05, 1.02) + ax.set_ylim(-1, len(d)) + ax.grid(axis="x", alpha=0.3) + fig.tight_layout() + fig.savefig(out_png, dpi=170, bbox_inches="tight") + print(f"\nwrote {out_png}") + + +def main(argv=None): + ap = argparse.ArgumentParser(description=__doc__.split("\n")[0]) + ap.add_argument("--years", type=int, nargs="+", default=[2025], + help="archive years to search (default 2025 -- the only " + "year with BAER dNBR at the time of writing)") + ap.add_argument("--out-dir", default="baer_study") + ap.add_argument("--workers", type=int, default=8) + args = ap.parse_args(argv) + + os.makedirs(args.out_dir, exist_ok=True) + matched = _matched_fires(args.years) + print(f"{len(matched)} fires carry both a BRISK composite and a BAER dNBR") + if not matched: + return None + + with ThreadPoolExecutor(args.workers) as ex: + rows = list(ex.map(lambda kv: compare_fire(kv[0], *kv[1]), + matched.items())) + df = pd.DataFrame([r for r in rows if r]) + csv = os.path.join(args.out_dir, "brisk_vs_baer.csv") + df.to_csv(csv, index=False) + print(f"wrote {csv}") + summarise(df) + if "r_d14" in df: + plot(df, os.path.join(args.out_dir, "brisk_vs_baer.png")) + return df + + +if __name__ == "__main__": + main() diff --git a/stormscape/__init__.py b/stormscape/__init__.py index d80b218..786a303 100644 --- a/stormscape/__init__.py +++ b/stormscape/__init__.py @@ -33,6 +33,9 @@ from .aoi import load_aoi, pad_bounds from .atlas14 import (anomaly, climatology_field, fetch_grid, grid_url, pf_point, recurrence_interval, region_for_bounds) +from .burn import (SEVERITY_SCHEMES, burn_severity, catalog as burn_catalog, + classify as burn_classify, fetch_scene, find_scenes, + register_brisk_cmap) from .compare import (comparison_stats, compare_storm, gauge_recurrence_table, radar_vs_gauge, sample_raster_at_points) from .merge import (conditional_merge, local_bias, loo_cross_validate, @@ -94,4 +97,6 @@ "reproject_geotiff", "export_geotiffs", "export_streams", "figure_to_geopdf", "geopdf_supported", "DEFAULT_EXPORT_FIELDS", "streams", "roads", "places", "reference_layers", + "burn_severity", "find_scenes", "burn_catalog", "fetch_scene", + "burn_classify", "SEVERITY_SCHEMES", "register_brisk_cmap", ] diff --git a/stormscape/burn.py b/stormscape/burn.py new file mode 100644 index 0000000..c1d52ee --- /dev/null +++ b/stormscape/burn.py @@ -0,0 +1,713 @@ +"""Near-real-time burn severity (dNBR) from CIMSS **BRISK**. + +Where :mod:`stormscape.mrms` gives the rain that falls *on* a burn scar, this +module gives the scar. **BRISK** ("Burned area Rapid Interim Severity risK", +UW-Madison CIMSS/SSEC) maps every large (>~5,000 acre) US wildfire **daily**, +which is what makes it usable while a fire is still burning: the authoritative +products land far too late for that -- BAER soil burn severity arrives days to +weeks after containment and only for fires that get an assessment, MTBS a year +or more later. BRISK closes the gap that matters for post-fire debris-flow +hazard, where the rain can arrive before the fire is out. + +It is a Google Earth Engine data-fusion dNBR composite over nine satellites +(GOES-East/West ABI, S-NPP / NOAA-20 / NOAA-21 VIIRS, Landsat 8/9, Sentinel-2a/b), +so a scar gets a coarse-but-immediate GOES look that sharpens as Landsat and +Sentinel-2 overpasses accumulate. + +Access +------ +The portal (``cimss.ssec.wisc.edu/brisk``) is a RealEarth viewer, and RealEarth's +WMTS tiles are **rendered PNG** -- pretty, but not data. The raw field is in the +open archive behind it, one GeoTIFF per fire per day:: + + https://bin.ssec.wisc.edu/pub/realearth/brisk//--dNBR__235959.tif + +so this module reads the Apache directory index, parses those names into a +catalog, and downloads only the scenes that intersect the AOI. Finding which +fires those are is cheap because GDAL can range-read a GeoTIFF header over HTTP +(``/vsicurl/``): ~0.02 s per scene across a thread pool, so a whole day's ~60 +fires are screened in under two seconds without downloading a pixel. + +The companion BAER **soil** burn severity archive (``baer-data/``, ``product= +"sbs"``) is read the same way -- see the warning about it below. + +Gotchas +------- +* **dNBR is a *vegetation* index, not soil burn severity.** The USGS post-fire + debris-flow models are calibrated on **soil** burn severity (BAER SBS), which + is dNBR *adjusted by field crews* for soil hydrophobicity, ground cover and + duff consumption. BRISK is explicitly an **interim** product: use it to act + early, then supersede it with SBS when the BAER assessment lands. +* **The scenes are EPSG:3857 at "60 m", which is not 60 m of ground.** Web + Mercator metres shrink with latitude, so a 60 m cell is ~46.5 m at 39 deg N and + ~42 m at 45 deg N. Fields are kept on their native 3857 grid (resampling the + science data to make a rounder number would only lose fidelity), and every + scene tested lands on an exact 60 m multiple, so mosaicking neighbouring fires + is a paste rather than a warp. +* **NaN marks the area outside the burn, but the files do not tag a nodata + value.** Reading with ``masked=True`` therefore masks nothing; test + ``np.isfinite`` instead. Typically only ~10-15% of a scene's pixels are valid. +* **dNBR is unscaled here** (about -0.3 to 1.0), not the x1000 integer form the + MTBS/USGS severity thresholds are usually quoted in. :data:`SEVERITY_SCHEMES` + holds the same breaks divided by 1000. +* A fire's footprint **grows day to day**, so the same fire's scenes differ in + extent; :func:`find_scenes` takes the latest per fire by default. + +Outputs mirror :mod:`stormscape.mrms`: a result ``dict`` (``fields / transform / +crs / profile / meta``) ready for :func:`stormscape.mrms.save_fields` and +:func:`stormscape.plot.drape_i15`. No new dependencies. +""" + +from __future__ import annotations + +import datetime as dt +import json +import os +import re +import time +import urllib.request +from concurrent.futures import ThreadPoolExecutor +from typing import Optional, Sequence + +import numpy as np +import pandas as pd + +from .aoi import bbox_polygon, load_aoi + +#: per-fire daily dNBR composites (CIMSS BRISK) +BRISK_BASE = "https://bin.ssec.wisc.edu/pub/realearth/brisk" +#: BAER soil burn severity, when an assessment exists (sparse) +BAER_BASE = "https://bin.ssec.wisc.edu/pub/realearth/baer-data" + +PRODUCTS = { + "dnbr": dict(base=BRISK_BASE, token="dNBR", kind="continuous", + label="dNBR", scale=1.0, first_year=2025), + "sbs": dict(base=BAER_BASE, token="sbs", kind="class", + label="soil burn severity class", scale=1.0, first_year=2025), + # The BAER teams' own dNBR, published beside their soil-burn-severity + # rasters. int16 **x1000** (the BARC convention -- NOAA's own + # ``BARC256 = dNBR*5 - 275`` identity confirms the scaling), so it is + # divided back to a plain index on read. 2025 only so far. + "baer_dnbr": dict(base=BAER_BASE, token="dNBR", kind="continuous", + label="BAER dNBR", scale=1000.0, first_year=2025), +} + +WORKERS = 12 # parallel header reads / downloads +CACHE_DIR = "brisk_cache" # sits at the event root, like nexrad_cache/ +INDEX_TTL_H = 6.0 # re-list the *current* year this often +TIMEOUT = 60.0 + +#: dNBR severity breaks. The USGS/MTBS thresholds are published on the x1000 +#: integer scale (100 / 270 / 440 / 660); BRISK's own portal colour scale steps +#: at 0.10 / 0.40 / 0.70. Class 0 is always "unburned or below the low break". +SEVERITY_SCHEMES = { + "usgs": dict( + breaks=(0.10, 0.27, 0.44, 0.66), + labels=("unburned", "low", "moderate-low", "moderate-high", "high")), + "brisk": dict( + breaks=(0.10, 0.40, 0.70), + labels=("unburned", "low", "moderate", "high")), +} + +#: Days a BRISK composite needs before its *magnitude* is trustworthy. +#: Measured against the BAER teams' own dNBR on the 39 fires that carry both +#: (``examples/brisk_vs_baer.py``): BRISK's pattern is right immediately, but a +#: composite scored on the BAER assessment date reads **low** until it has +#: ingested a clear post-fire Landsat/Sentinel overpass. Every poor performer +#: recovered to r = 0.80-0.96 given a scene 5-21 days later, and a fixed +14 d +#: rule took fires at r >= 0.90 from 14/39 to 24/39. Advisory by default; a hard +#: filter with ``--min-age``. +MATURITY_DAYS = 14 + +#: BAER soil burn severity is delivered already classified. +SBS_LABELS = {1: "unburned/very low", 2: "low", 3: "moderate", 4: "high"} +#: Valid SBS class range. The rasters carry an embedded palette that colours 1-4 +#: with BRISK's own four severity colours and paints **0 and 5+ the same black**, +#: i.e. the product itself treats anything outside 1-4 as not-a-severity (a +#: water/inholding/unmapped mask, ~3% of the one NV scene checked). So values +#: outside this range are read as missing rather than charted as a class. +SBS_VALID = (1, 4) + +#: The **BAER burn-severity class colours**, verified from the products +#: themselves: all 77 of the 2025 soil-burn-severity rasters -- written by many +#: different BAER teams, in ERDAS Imagine -- carry this identical embedded +#: palette for classes 1-4 (only the class-5 mask colour varies). BRISK's own +#: ``qgis_BRISK_dNBR_colorscale_v2.txt`` uses exactly these four colours, so the +#: portal and the BAER deliverables already share one scheme; this is it. +BAER_CLASS_COLORS = ((0, 128, 128), # 1 unburned / very low teal + (82, 204, 204), # 2 low cyan + (255, 232, 32), # 3 moderate yellow + (168, 0, 0)) # 4 high dark red + +#: The same palette as a continuous dNBR ramp, as (dNBR, RGB) anchors -- the +#: interpolated form BRISK publishes, with its hard step at the 0.10 unburned +#: break. Registered as the ``"baer"`` colormap over [0, 1]. +BAER_ANCHORS = ((0.00, (0, 128, 128)), (0.10, (0, 128, 128)), + (0.10, (82, 204, 204)), (0.40, (255, 232, 32)), + (0.70, (168, 0, 0)), (1.00, (114, 0, 0))) +BRISK_ANCHORS = BAER_ANCHORS # back-compat alias; one scheme + +_HREF = re.compile(r'href="([^"?/][^"]*\.tif)"') + + +# --------------------------------------------------------------------------- # +# catalog: parse the archive's directory index +# --------------------------------------------------------------------------- # +def _product(product: str) -> dict: + try: + return PRODUCTS[product] + except KeyError: + raise ValueError(f"product must be one of {sorted(PRODUCTS)}, " + f"got {product!r}") from None + + +def parse_name(filename: str, token: str = "dNBR") -> Optional[dict]: + """Split ``Hidden-Valley-NV-dNBR_20260814_235959.tif`` into its parts. + + Returns ``{fire, state, date, filename}`` or ``None`` if the name does not + match. The state is the trailing two-letter code when present -- a few + scenes are named by incident number alone (``0231-OR``) or carry a + non-state code, so it is informational, never a filter. + """ + m = re.match(rf"^(?P.+?)-{re.escape(token)}_" + r"(?P\d{8})_(?P