name: yellowstone-burn-recovery
version: "0.2"
category: change-detection
themes: [vegetation]
family: fire-burn-severity-dnbr-recovery
description: >
  Two-epoch NBR over the 1988 Yellowstone fires, which burned about 3,200 km²
  of the park in one season. The baseline epoch is the pre-fire forest and the
  recent epoch is today's canopy, so the difference reads as RECOVERY rather
  than severity: the stand-replacing lodgepole burn has been regrowing for more
  than three decades.

# WHAT THIS MEASURES, AND WHAT IT DOES NOT. READ BEFORE CITING A NUMBER.
#   NBR is a CHAR-AND-MOISTURE PROXY. A recovered NBR says the canopy is green
#   again; it does not say the forest has recovered its structure, species mix
#   or carbon. The answer below is a count of pixels still below the recovery
#   threshold at 30 m (0.09 ha each), not a mapped unrecovered burn area.
#
#   This is NOT a dNBR severity product. Classic dNBR is pre-fire minus
#   IMMEDIATELY post-fire; the recent epoch here is 30+ years later, so the
#   difference confounds severity with three decades of regrowth. Read it as a
#   recovery deficit, and nothing else.
#
#   No per-year axis. Each epoch is ONE median composite over its window, so the
#   answer is a two-epoch comparison, not an annual series. A genuine annual
#   series needs a per-year grouping reduction, and no shipped op provides one
#   (recipe-library tracker, defect 63).
#
#   UPSTREAM BUG (recipe-library tracker, defects 61/62). Some Landsat 7 (LE07)
#   assets on the Planetary Computer answer an HTML error page where a COG
#   should be, and the read dies as "Unknown byte order marker: [3C, 21]" (`<!`).
#   Confirmed 2026-08-17 on LE07_L2SP_036035_20221027_02_T1. There is no
#   platform filter to exclude L7 with: the source params carry none, and the
#   `query.platform` block the previous version of this recipe declared was
#   silently ignored. The only lever an author has is the calendar. Windows that
#   overlap 1999-2023 can pull L7 and can therefore hit this.
#
#   No month filter. The engine's source params carry no month selector, so the
#   seasonal windows below are calendar spans and off-season scenes do enter the
#   median. A median over many scenes is robust to that; a mean would not be.

settings:
  nbr_recovery_threshold: -0.1   # dNBR below this → still short of the pre-fire canopy
  cloud_cover_max_pct: 20        # eo:cloud_cover upper bound on the STAC search
  pixel_area_m2: 900             # 30 m Landsat pixel — multiply `count` by this for area


inputs:
  aoi:
    bbox: [-111.2, 44.3, -109.8, 45.1]
    crs_target: EPSG:32612
    notes: >
      Northwestern Yellowstone, covering the North Fork and Clover-Mist burn
      perimeters. Summer windows only: snow cover destroys NBR.


  baseline_epoch:
    start: "1985-07-01"
    end: "1988-06-30"
    notes: Landsat 5 TM, three pre-fire summers ending before ignition.

  recent_epoch:
    start: "2021-07-01"
    end: "2025-09-30"
    notes: Landsat 8/9 OLI, the present canopy.

  collection:
    uri: "@usgs/landsat/c2-l2"

layers:

  # ── SOURCE ──────────────────────────────────────────────────────────────────
  #
  # Bands are CANONICAL concept codes from registry/datasets/usgs/landsat/c2-l2/
  # _dataset.yaml (N = nir, S2 = swir-2, QA = qa-pixel). Each distribution maps them onto its own
  # producer asset keys, so ONE band list is correct across TM, ETM+ and OLI —
  # which the previous SR_Bx spelling was not (SR_B5 is SWIR-1 on TM/ETM+ and
  # NIR on OLI, so the old files asked for a different physical band per sensor).
  #
  # The FUSED per-scene form (`scenes: true` + `expression:` + `mask_band:`) is
  # what makes the index correct: the distribution's declared scale/offset is
  # applied per band before the expression, fill DNs become NaN, and QA_PIXEL is
  # the platform mask band, excluded from the b1..bN binding. `temporal_reduce`
  # then takes the median across the surviving scenes.
  #
  # Do NOT collapse this to a bare `bands: [A, B, QA]` source and difference the
  # result. A source read emits the DECLARED BANDS (measured 2026-08-17: three
  # bands out for a three-band request), so `b2 - b1` there is raw band math,
  # not a normalized-difference change. It plans green and answers plausibly.

  source/nbr-baseline-scenes:
    type: raster
    description: Per-scene NBR, QA-masked, baseline epoch.
    uri: "@usgs/landsat/c2-l2"
    params:
      scenes: true
      bands: [N, S2, QA]
      mask_band: QA
      expression: "(b1 - b2) / (b1 + b2)"
      bbox: [-111.2, 44.3, -109.8, 45.1]
      temporal: ["1985-07-01", "1988-06-30"]
      cloud_lt: 20
      resolution: 30
      limit: 60

  source/nbr-recent-scenes:
    type: raster
    description: Per-scene NBR, QA-masked, recent epoch.
    uri: "@usgs/landsat/c2-l2"
    params:
      scenes: true
      bands: [N, S2, QA]
      mask_band: QA
      expression: "(b1 - b2) / (b1 + b2)"
      bbox: [-111.2, 44.3, -109.8, 45.1]
      temporal: ["2021-07-01", "2025-09-30"]
      cloud_lt: 20
      resolution: 30
      limit: 60

  # ── COMPUTE ─────────────────────────────────────────────────────────────────

  compute/nbr-baseline:
    type: raster
    description: Median NBR over the baseline epoch's scenes.
    compute:
      op: temporal_reduce
      inputs:
        time_series: { layer: source/nbr-baseline-scenes }
      params:
        method: median

  compute/nbr-recent:
    type: raster
    description: Median NBR over the recent epoch's scenes.
    compute:
      op: temporal_reduce
      inputs:
        time_series: { layer: source/nbr-recent-scenes }
      params:
        method: median

  compute/nbr-pair:
    type: raster
    description: The two NBR composites stacked into one 2-band raster (b1 = baseline, b2 = recent).
    compute:
      op: stack_concat
      inputs:
        layers:
          - { layer: compute/nbr-baseline }
          - { layer: compute/nbr-recent }

  compute/nbr-change:
    type: raster
    description: >
      dNBR = recent - baseline. Negative where today's canopy is still short of
      the pre-fire stand; near zero where it has caught up.
    compute:
      op: raster_calc
      inputs:
        raster: { layer: compute/nbr-pair }
      params:
        expression: "b2 - b1"
        output_band: dnbr

  compute/deficit-baseline:
    type: raster
    description: >
      Recovery-deficit footprint. Pixels still below the pre-fire canopy keep
      the value 1; recovered pixels become NaN (`b1 / 0.0`), which the stats
      kernel skips — so the answer layer's `count` is the deficit-pixel count.
    compute:
      op: raster_calc
      inputs:
        raster: { layer: compute/nbr-baseline }
      params:
        expression: "where(b1 < ${settings.nbr_recovery_threshold}, 1.0, b1 / 0.0)"
        output_band: deficit

  compute/deficit-recent:
    type: raster
    description: >
      Recovery-deficit footprint. Pixels still below the pre-fire canopy keep
      the value 1; recovered pixels become NaN (`b1 / 0.0`), which the stats
      kernel skips — so the answer layer's `count` is the deficit-pixel count.
    compute:
      op: raster_calc
      inputs:
        raster: { layer: compute/nbr-recent }
      params:
        expression: "where(b1 < ${settings.nbr_recovery_threshold}, 1.0, b1 / 0.0)"
        output_band: deficit

  # ── ANSWERS ─────────────────────────────────────────────────────────────────

  answer/deficit-baseline:
    type: answer
    description: >
      Recovery-deficit pixel count against the baseline composite (the baseline
      arm is a self-check and should be near zero). Area = count x 900 m².
    compute:
      op: summary
      inputs:
        raster: { layer: compute/deficit-baseline }
      params:
        stats: [count]

  answer/deficit-recent:
    type: answer
    description: >
      Recovery-deficit pixel count against the recent composite (the baseline
      arm is a self-check and should be near zero). Area = count x 900 m².
    compute:
      op: summary
      inputs:
        raster: { layer: compute/deficit-recent }
      params:
        stats: [count]

outputs:
  deficit_recent:
    layer: answer/deficit-recent
    value_kind: answer
    description: Pixels whose canopy is still short of the pre-fire stand. THE answer.
  deficit_baseline:
    layer: answer/deficit-baseline
    value_kind: answer
    description: The same test applied to the baseline against itself — a near-zero self-check.
  nbr_change:
    layer: compute/nbr-change
    value_kind: raster
    description: dNBR raster, recent minus pre-fire (show-your-work).


license:
  recipe: Apache-2.0
  output_dataset: CC-BY-4.0
  upstream_data:
    - { source: "@usgs/landsat/c2-l2", license: usgs-public-domain }