name: crop-anomaly version: "0.1" status: draft # compute graph authored; NOT yet run through # the engine. Validate with `folia run` on the # example AOI. See PROVENANCE. category: time-series-trends themes: [agriculture, imagery] family: crop-phenology-monitoring cites_papers: - verbesselt-2010 # BFAST: detecting and characterizing changes in trend and seasonal components - zhu-2014 # Continuous change detection and classification (CCDC) description: > Per-field crop anomaly detection from a Sentinel-2 NDVI time series. For an AOI covered by Fields of the World boundaries, fit a harmonic seasonal model to each pixel's NDVI history, compute the residual between the latest observed NDVI and the model-expected value, and aggregate that anomaly residual to each field. Fields running persistently below their own fitted season are flagged as anomalous. The method follows the BFAST season-anomaly framing (Verbesselt et al. 2010): the harmonic fit captures the expected phenological trajectory, and the signed residual (observed minus fitted) is the anomaly signal. CCDC-style persistence logic (Zhu & Woodcock 2014) informs the threshold and flagging step. The per-field anomaly table is the primary product; the NDVI stack, fitted surface, and residual raster are the scrutable "show your work" layers. Swap the NDVI expression for EVI/NDRE to shift the vigor proxy, or tighten the anomaly threshold to control alert sensitivity. # ───────────────────────────────────────────────────────────────────────────── # PROVENANCE (why this recipe exists) # Authored 2026-07-14 from the crop white-space suite (see crop-monitoring). Where # crop-monitoring delivers the fitted season itself (amplitude/phase/trend), this # recipe takes the fit as a *baseline* and computes the signed departure from it — # the BFAST season-anomaly loop (verbesselt-2010). The CCDC persistence model # (zhu-2014) informs the flag threshold: a field is anomalous when its mean # residual falls below –sigma, matching CCDC's break-evidence accumulation logic. # # This is the alert-facing counterpart to crop-monitoring: monitoring tells you # what the season looks like; anomaly tells you which fields are off their own # baseline right now. # # STATUS: draft. The op composition follows each op's documented signature. The # recipe has NOT been run through the engine. Validate with `folia run` on the # example AOI, reconcile any param/shape mismatches, add a validations/ block, # then drop status:draft. # ───────────────────────────────────────────────────────────────────────────── # ───────────────────────────────────────────────────────────────────────────── # PARAMETERIZATION # settings.index: the greenness index (NDVI) — formula + bands. Swap for # EVI/NDRE to change the vigor proxy. # settings.n_harmonics / fit_trend: the phenology model. 2 harmonics captures a # main season + secondary (double-cropping). # settings.anomaly_threshold: z-like residual cutoff for per-field flagging. # Negative: below-baseline anomaly (stress, failure). Tune # toward 0 to widen the alert, toward –1 to tighten it. # inputs.fields: @ftw/global-field-boundaries — the per-field zones. The # anomaly residual is aggregated into these polygons; fields # carry a per-field confidence attribute (ftw_conf) usable as # a weight or post-hoc filter. # inputs.aoi / date_range: the bound instance. # ───────────────────────────────────────────────────────────────────────────── settings: index: slug: ndvi expression: "((B08 - reflectance_offset) - (B04 - reflectance_offset)) / ((B08 - reflectance_offset) + (B04 - reflectance_offset))" bands: [B04, B08] # both native 10m, no grid alignment needed range: [-1, 1] colormap: rdylgn # diverging: red = stress, green = vigorous n_harmonics: 2 # annual + semi-annual (double-cropping) fit_trend: true # remove multi-year drift before computing residual period: 365.25 # day-of-year time axis anomaly_threshold: -0.5 # residual cutoff (z-like); fields below = flagged field_stats: [mean, std, min, count] # per-field aggregation of the anomaly residual reflectance_offset: 1000 # >=04.00 BOA offset (DN - offset); scale cancels in ND max_bad_fraction: 0.30 # AOI-level cloud gate bad_scl_classes: [0, 1, 3, 8, 9, 10] inputs: aoi: name: San Joaquin Valley cropland block bbox: [-120.05, 36.55, -119.85, 36.72] crs_target: EPSG:32611 mgrs_tile: "11SKB" fields: # Per-field zones sourced from Fields of the World (ftw). Each polygon carries # ftw_conf (field-boundary confidence, 0–1) — use as a post-hoc filter or weight # when ranking anomaly flags. Binding "@ftw/global-field-boundaries" pulls the # FTW vector tile covering the AOI bbox at query time. uri: "@ftw/global-field-boundaries" description: > Fields of the World global field-boundary dataset (CC-BY-4.0). The per-field confidence attribute (ftw_conf) indicates boundary quality; prefer fields with ftw_conf >= 0.7 for production alert runs. date_range: start: "2023-01-01" # >= 2 full growing seasons for a stable harmonic fit end: "2026-06-12" collection: uri: "@esa/sentinel-2/l2a" assets: [B04, B08, SCL] # NDVI bands + SCL cloud gate layers: # ── SOURCE ────────────────────────────────────────────────────────────────── source/s2: type: raster_timeseries uri: "@esa/sentinel-2/l2a" description: > Every clear S2 L2A pass over the AOI, one item per acquisition date, gated at the AOI-polygon level from its SCL window before the NDVI bands are read. The SCL gate removes cloud, cirrus, and snow pixels before any index is computed, so the NDVI stack reflects only valid surface observations. params: aoi: { bbox: "${inputs.aoi.bbox}" } datetime: "${inputs.date_range.start}/${inputs.date_range.end}" assets: "${inputs.collection.assets}" dedup_by: acquisition_date aoi_cloud_gate: scl_asset: SCL bad_classes: "${settings.bad_scl_classes}" max_bad_fraction: "${settings.max_bad_fraction}" # ── COMPUTE ────────────────────────────────────────────────────────────────── compute/ndvi: type: raster_timeseries description: > NDVI per scene (raster_calc temporal mode over B04/B08). This (t, y, x) stack is the scrutable greenness layer and the sole input to the harmonic fit. Using raster_calc rather than normalized_diff keeps the BOA reflectance offset handling explicit and consistent with the crop-monitoring recipe. compute: op: raster_calc inputs: data: { layer: source/s2 } params: expression: "${settings.index.expression}" reflectance_offset: "${settings.reflectance_offset}" output_band: ndvi compute/expected: type: raster description: > Per-pixel harmonic seasonal model fit to the NDVI stack (verbesselt-2010). harmonic_fit returns three output bands: fitted (the model-predicted NDVI at every timestep), amplitude (annual greenness swing), and phase (day-of-year of peak greenness). fit_trend removes the multi-year greening/browning signal before the residual is computed so the anomaly reflects within-season departure, not a structural trend. The fitted band is the direct input to the anomaly differencing step. compute: op: harmonic_fit inputs: collection: { layer: compute/ndvi } params: n_harmonics: "${settings.n_harmonics}" fit_trend: "${settings.fit_trend}" period: "${settings.period}" compute/anomaly: type: raster_timeseries description: > Signed residual: observed NDVI minus harmonic-fitted NDVI, at each timestep (BFAST season-anomaly framing, verbesselt-2010). Negative values = below- baseline stress; positive values = above-baseline vigor. raster_calc in two-input mode subtracts the fitted surface from the observed stack. The latest timestep slice is the primary anomaly raster shown on the map; the full (t,y,x) residual cube is retained for show_your_work. # NOTE: the expression below references the two input bands by their output_band # names. If the engine requires positional b1/b2 syntax instead of named bands, # replace "(ndvi - fitted)" with "(b1 - b2)" and map inputs accordingly. compute: op: raster_calc inputs: data: { layer: compute/ndvi } fitted: { layer: compute/expected, band: fitted } params: expression: "(ndvi - fitted)" output_band: anomaly_residual compute/field_flags: type: table description: > Per-field anomaly summary: zonal_stats of the anomaly residual raster over @ftw/global-field-boundaries, returning mean/std/min/count of the residual inside each polygon. A derived column flags fields whose mean residual falls below settings.anomaly_threshold (CCDC-style persistence logic, zhu-2014): these are the alert rows surfaced to the dashboard. The ftw_conf attribute carried by the @ftw zones propagates to the output table for post-hoc quality filtering. # NOTE: the flag column (mean_residual < threshold -> is_anomaly bool) is a # table-level derived column. If the engine does not support inline derived # columns in zonal_stats params, add a reclassify layer over the output table # or compute it in the dashboard pattern. No dedicated "table_calc" op exists. compute: op: zonal_stats inputs: raster: { layer: compute/anomaly } zones: { feature: "${inputs.fields.uri}" } params: stats: "${settings.field_stats}" # Carry through the ftw_conf attribute from the @ftw zones feature so # downstream consumers can filter by field-boundary confidence. carry_zone_attributes: [ftw_conf] flag_column: name: is_anomaly condition: "mean_anomaly_residual < ${settings.anomaly_threshold}" outputs: primary: - compute/field_flags # geojson/table — per-field anomaly flags (THE product) show_your_work: - compute/anomaly # zarr (t,y,x) — full residual cube; latest slice on map - compute/expected # cog (fitted, amplitude, phase) — the harmonic baseline - compute/ndvi # zarr (t,y,x) — the raw NDVI stack behind the fit provenance: - source_scene_ids # STAC item id per timestamp # ── DASHBOARD ─────────────────────────────────────────────────────────────── dashboard: patterns: - slug: anomaly-vs-climatology # anomaly residual vs the fitted season curve slot: temporal layer: compute/anomaly - slug: zonal-rank # fields ranked by mean anomaly; table + choropleth slot: aggregation layer: compute/field_flags - slug: derived-alert-subscribe # subscribe to field-flag alerts via threshold slot: subscription layer: compute/field_flags license: recipe: Apache-2.0 output_dataset: CC-BY-4.0 upstream_data: - { source: sentinel-2-l2a, license: usgs-public-domain } - { source: ftw-global-field-boundaries, license: CC-BY-4.0 } # Related: crop-monitoring (the fitted-season product this anomaly recipe uses as # its baseline); index-timeseries (the source->index skeleton both extend). # Paper provenance: verbesselt-2010 (BFAST season-anomaly method); # zhu-2014 (CCDC persistence/flagging logic). # White-space provenance: crop white-space suite (crop-anomaly / crop-stress / # crop-yield-forecast), authored as the alert-facing complement to crop-monitoring.