# Five canonical spectral indices over an AOI, each a cloud-masked median composite.
#
# ORDER OF OPERATIONS, because it is the whole scientific content of this file:
# every index here is computed PER SCENE and only then reduced across time. It is
# median-of-the-index, never index-of-the-median-bands. That is what the engine's
# composite actually does (`composite.rs::scene_index`: reflectance -> index -> cloud
# mask, per scene, then `reduce_scenes`), and it is the defensible ordering.
#
# The distinction is not pedantic once you change the reducer. For a MAXIMUM VALUE
# COMPOSITE (MVC) or any percentile, the two orderings give different answers, and
# only this one is the statistic you asked for: the max of an index is NOT the index
# of the max-composited bands, because the band that maximises a ratio is not the
# maximum of either band. A recipe that composites reflectance first and then divides
# is silently answering a different question. If you are writing an MVC, keep this
# per-scene shape and change `method:` — do not restructure it.
#
# THIS FILE WAS BROKEN AND IS KEPT AS THE FIXED REFERENCE. Until 2026-08-10 it declared
# a `raster_calc` layer with no `expression:` (the op's one required param), fed it
# `scl_mask:` and `target_resolution_m:` params that no op reads, and passed
# `aoi:`/`datetime:`/`assets:` to a source layer that wanted `bbox:`/`temporal:`/`bands:`.
# It planned green the whole time: `folia plan` is structural and checks neither op
# params nor param names. It would have failed on the first real run.
#
# WHAT DOES THE WORK, so you can see there is no magic left in it:
#   - reflectance:  the DATASET owns the DN -> reflectance transform (S2 L2A is
#                   x*0.0001 - 0.1, in registry/datasets/esa/sentinel-2/l2a/_dataset.yaml),
#                   applied per band BEFORE the index. No recipe-level offset arithmetic;
#                   the earlier hand-rolled `(DN - 1000)/10000` was a second, divergent
#                   copy of a number the registry already holds.
#   - cloud mask:   declaring SCL as the LAST band makes it the platform QA band, and the
#                   Sentinel-2 kernel NaN-s SCL 3 / 8 / 9 / 10 (shadow, cloud medium, cloud
#                   high, cirrus) per scene. The class set is compiled in and is not
#                   overridable here; if you need a different one, build the mask
#                   explicitly (docs/RECIPE_AUTHORING.md §6, and
#                   registry/recipes/openeo-scl-cloud-mask).
#   - reduction:    NaN-aware, so a pixel cloudy in every pass answers NaN rather than a
#                   fabricated value.
#
# Four of the five indices are plain normalized differences and need only a band pair, so
# they are one layer each. BSI is not, so it uses the per-scene expression form and an
# explicit `temporal_reduce` — the same shape as registry/recipes/eefa-f40-evi.
#
# Canonical band codes (registry/datasets/esa/sentinel-2/l2a/_dataset.yaml), NOT asset keys:
#   B = blue (B02)   G = green (B03)   R = red (B04)   N = nir (B08)
#   S1 = swir-1 (B11)   S2 = swir-2 (B12)   SCL = scene classification
# Note `S2` is swir-2, not "Sentinel-2".

name: spectral-indices
version: "0.1"
category: time-series-trends
themes: [vegetation, hydrology, imagery]
family: spectral-index-map-timeseries
description: >
  Five canonical spectral indices over a named AOI from Sentinel-2 L2A: NDVI, NDSI,
  NDWI, NBR and BSI. Each is computed per scene, cloud-masked from the scene
  classification layer, and reduced to a median over the date window. One snapshot per
  run, not a time series; rerunning on a schedule updates the snapshot.

  NDVI = (N - R)   / (N + R)      vegetation greenness (NIR vs red)
  NDSI = (G - S1)  / (G + S1)     snow / cloud discrimination (green vs SWIR-1)
  NDWI = (G - N)   / (G + N)      surface water (McFeeters; green vs NIR)
  NBR  = (N - S2)  / (N + S2)     burn severity / health (NIR vs SWIR-2)
  BSI  = ((S1 + R) - (N + B)) / ((S1 + R) + (N + B))     bare soil (SWIR-1 + red vs NIR + blue)

  On NDSI: this recipe treats snow as a real surface signal rather than something to
  remove, so nothing here masks the snow class. The platform cloud mask drops cloud,
  shadow and cirrus only. For snow-cover fraction as a question in its own right, use
  snow-cover-fsc, which drives per-scene NDSI through a sigmoid calibrated to FSCTOC.

settings:
  # Scene-level gate on STAC `eo:cloud_cover`: skip an acquisition before reading it.
  # Coarse and cheap, and independent of the per-pixel SCL mask, which still runs on
  # every scene that passes. Partly-cloudy scenes are worth keeping: they contribute
  # wherever they ARE clear, and the median ignores the NaNs.
  cloud_max: 60

  # Common output grid, in metres. B02/B03/B04/B08 are natively 10 m and B11/B12 are
  # 20 m, so 20 m is the finest grid every band in this file can answer at honestly.
  resolution_m: 20

  # Cap on scenes read per index. Raise for a longer window.
  scene_limit: 120

inputs:
  aoi:
    # Salt Lake region, Utah.
    bbox: [-113.0, 40.0, -111.0, 41.5]
  date_range:
    # Absolute on purpose. A relative window ("last 90 days") has four different
    # spellings across the tree and no canonical one yet, so a recipe that uses one
    # is not portable.
    start: "2024-06-01"
    end: "2024-08-31"

layers:

  # ── NORMALIZED-DIFFERENCE INDICES ────────────────────────────────────────────────
  # One layer each. `bands: [a, b, SCL]` means: per scene compute (a - b) / (a + b) on
  # reflectance, NaN-out the cloud classes using SCL, then take the median across scenes.
  # The band ORDER is the index definition — swapping the pair negates the result.

  compute/ndvi:
    type: raster
    uri: "@esa/sentinel-2/l2a"
    description: NDVI = (N - R) / (N + R), per-scene, median over the window. Range [-1, 1].
    params:
      bands: [N, R, SCL]
      bbox: "${inputs.aoi.bbox}"
      temporal: ["${inputs.date_range.start}", "${inputs.date_range.end}"]
      method: median
      resolution: "${settings.resolution_m}"
      cloud_lt: "${settings.cloud_max}"
      limit: "${settings.scene_limit}"
    style:
      palette: rdylgn
      range: [-0.2, 0.9]

  compute/ndsi:
    type: raster
    uri: "@esa/sentinel-2/l2a"
    description: NDSI = (G - S1) / (G + S1), per-scene, median over the window. Range [-1, 1].
    params:
      bands: [G, S1, SCL]
      bbox: "${inputs.aoi.bbox}"
      temporal: ["${inputs.date_range.start}", "${inputs.date_range.end}"]
      method: median
      resolution: "${settings.resolution_m}"
      cloud_lt: "${settings.cloud_max}"
      limit: "${settings.scene_limit}"

  compute/ndwi:
    type: raster
    uri: "@esa/sentinel-2/l2a"
    description: NDWI = (G - N) / (G + N), McFeeters, per-scene, median. Range [-1, 1].
    params:
      bands: [G, N, SCL]
      bbox: "${inputs.aoi.bbox}"
      temporal: ["${inputs.date_range.start}", "${inputs.date_range.end}"]
      method: median
      resolution: "${settings.resolution_m}"
      cloud_lt: "${settings.cloud_max}"
      limit: "${settings.scene_limit}"

  compute/nbr:
    type: raster
    uri: "@esa/sentinel-2/l2a"
    description: NBR = (N - S2) / (N + S2), per-scene, median. Range [-1, 1].
    params:
      bands: [N, S2, SCL]
      bbox: "${inputs.aoi.bbox}"
      temporal: ["${inputs.date_range.start}", "${inputs.date_range.end}"]
      method: median
      resolution: "${settings.resolution_m}"
      cloud_lt: "${settings.cloud_max}"
      limit: "${settings.scene_limit}"

  # ── BSI ──────────────────────────────────────────────────────────────────────────
  # Not a normalized-difference pair, so it takes the general per-scene expression form:
  # `scenes: true` for a real time axis, a fused `expression:` evaluated per acquisition
  # as the bytes decode, and `mask_band:` naming SCL as the QA band (which excludes it
  # from the b1..bN binding). The reduce is then explicit and separate, which is what
  # makes the intermediate stack inspectable.
  #
  # Bands bind POSITIONALLY in declared order: b1 = S1, b2 = R, b3 = N, b4 = B.

  compute/bsi-scenes:
    type: raster
    uri: "@esa/sentinel-2/l2a"
    description: >
      Per-scene BSI, one band per acquisition, cloud-masked. The show-your-work layer:
      scrub it to see which passes each pixel of the composite was built from.
    params:
      scenes: true
      bands: [S1, R, N, B, SCL]
      mask_band: SCL
      expression: "((b1 + b2) - (b3 + b4)) / ((b1 + b2) + (b3 + b4))"
      bbox: "${inputs.aoi.bbox}"
      temporal: ["${inputs.date_range.start}", "${inputs.date_range.end}"]
      resolution: "${settings.resolution_m}"
      cloud_lt: "${settings.cloud_max}"
      limit: "${settings.scene_limit}"

  compute/bsi:
    type: raster
    description: BSI median over the per-scene stack. Range [-1, 1].
    compute:
      op: temporal_reduce
      inputs: { data: { layer: compute/bsi-scenes } }
      params:
        method: median

outputs:
  ndvi:
    layer: compute/ndvi
    value_kind: raster
    description: Median NDVI composite.
  ndsi:
    layer: compute/ndsi
    value_kind: raster
    description: Median NDSI composite.
  ndwi:
    layer: compute/ndwi
    value_kind: raster
    description: Median NDWI composite.
  nbr:
    layer: compute/nbr
    value_kind: raster
    description: Median NBR composite.
  bsi:
    layer: compute/bsi
    value_kind: raster
    description: Median BSI composite.
  bsi_stack:
    layer: compute/bsi-scenes
    value_kind: raster
    description: Show-your-work — the per-scene BSI stack behind the BSI composite.

license:
  recipe: Apache-2.0
  output_dataset: CC-BY-4.0
  upstream_data:
    - { source: "@esa/sentinel-2/l2a", license: CC-BY-SA-3.0-IGO }