Use WorldCover-aware FANO for SSEBop cold references

README.md

@@ -55,8 +55,8 @@ The five stages share one implementation in `pysweb.workflow`:
 
 | Stage | Role | Inputs / requirements |
 |---|---|---|
-| `prepare` | Download Landsat, ERA5-Land and NASADEM; stack meteorology | Earth Engine authentication and your project |
-| `ssebop` | Estimate daily ET and partition E_s/E_t | Prepared satellite, meteorology and DEM files |
+| `prepare` | Download Landsat, ERA5-Land, NASADEM and WorldCover; stack meteorology | Earth Engine authentication and your project |
+| `ssebop` | Estimate daily ET with land-cover-aware FANO and partition E_s/E_t | Prepared satellite, meteorology, DEM and WorldCover files |
 | `preprocess` | Align forcing and soil properties; calculate effective rain | Full calendar-month rain; local or online soil inputs |
 | `calibrate` | Fit two domain parameters to reference surface soil moisture | Enabled explicitly; local reference or configured Earth Engine asset |
 | `run` | Simulate independent soil columns | Prepared forcing and hydraulic properties |

docs/methods.md

@@ -2,14 +2,41 @@
 
 ## Process chain
 
-SSEBop uses Landsat surface temperature, vegetation information and meteorological reference evapotranspiration to estimate ET fraction and daily ET. The maintained local cold-temperature/FANO correction uses the existing configurable NDVI thresholds and spatial scales in `pysweb.ssebop.core`. Daily NDVI partitions ET using `f_t = clip(1.26 * NDVI - 0.18, 0, 1)`, `E_t = f_t * ET_a`, and `E_s = ET_a - E_t`. The forcing identity is `ET_a = ET_f * ET_0 = E_s + E_t`; NetCDF total ET retains its `ET`/`et` names. This empirical partition requires evaluation for the vegetation and season being studied.
+SSEBop uses Landsat surface temperature, vegetation information and meteorological reference evapotranspiration to estimate ET fraction and daily ET. The default cold boundary is the WorldCover FANO adaptation described below; the previous local correction is available with `tcold_method="legacy"`. Daily NDVI partitions ET using `f_t = clip(1.26 * NDVI - 0.18, 0, 1)`, `E_t = f_t * ET_a`, and `E_s = ET_a - E_t`. The forcing identity is `ET_a = ET_f * ET_0 = E_s + E_t`; NetCDF total ET retains its `ET`/`et` names. This empirical partition requires evaluation for the vegetation and season being studied.
 
 The soil model uses conservative finite volumes with physical storage bounds, beta-distributed root uptake and a bottom drainage boundary. Evaporation draws from the surface layer and can continue below WP; root extraction stops at WP. A fixed three-day supply closure limits E/T, and joint donor constraints prevent transfers and sinks from spending the same water. Physical field capacity initialises the profile. The exact equations, frozen GSSM reference derivation, daily update and migration contract are recorded in [the formulation specification](swb-formulation.md).
 
 Separate `actual_evaporation`, `actual_transpiration`, `unmet_evaporation` and `unmet_transpiration` budgets distinguish forcing demand from water-limited extraction. `actual_et` is their realised E/T sum. The `E_s` and `E_t` SEB variables retain their forcing meanings. No interception store or flux is introduced.
 
 OpenLandMap clay and sand predictors are percentages. Organic carbon uses the asset scale of 5 g kg⁻¹ per encoded unit. The existing mapping samples bands at 0, 10, 30, 60 and 100 cm and assigns them to model bottoms at 50, 150, 300, 600 and 1000 mm. These source point-depth predictions and model layer averages have different vertical support. Missing organic carbon uses the existing 5 g kg⁻¹ fallback, recorded in preprocessing configuration; clay/sand gaps remain missing. The local MLConstraints organic-matter convention is documented separately in [soil inputs](soil-inputs.md).
 
+## Land-cover-aware FANO cold boundary
+
+`pysweb.ssebop.fano.worldcover_fano` follows the fine/coarse and mixed-landscape logic in [OpenET SSEBop](https://github.com/Open-ET/openet-ssebop/blob/main/openet/ssebop/image.py), inspected 22 October 2026 (Git blob `5aef6a7a3833679a2a1598ef33b2e7adbdc25cdd`). The implementation is an Australian land-cover adaptation, not an exact reproduction of OpenET ET products.
+
+The input classification is [ESA WorldCover v200](https://developers.google.com/earth-engine/datasets/catalog/ESA_WorldCover_v200), band `Map`, a global 10 m **2021** map. It is static for every Landsat date, including earlier years. Land-cover change and classification error therefore remain sources of uncertainty. Credit: Zanaga et al. (2022), [doi:10.5281/zenodo.7254221](https://doi.org/10.5281/zenodo.7254221), CC BY 4.0.
+
+| WorldCover classes | FANO use |
+| --- | --- |
+| 30 grassland, 40 cropland, 90 herbaceous wetland | Agricultural branch; eligible for regional references when native NDVI ≥ 0.4 |
+| 20 shrubland, 50 built-up, 60 bare/sparse vegetation | Local fine branch; anomalous negative-NDVI treatment |
+| 10 tree cover, 70 snow/ice, 95 mangroves, 100 moss/lichen | Local fine branch; excluded from the agricultural reference mask |
+| 80 permanent water | Excluded from terrestrial ET and all cold-reference means |
+| Missing or unknown class | Error wherever the scene otherwise has valid LST, NDVI and positive dT |
+
+This positive reference mask excludes more than forests. OpenET's NLCD mapping also includes developed open/low-intensity areas and woody wetlands. WorldCover cannot distinguish those developed subclasses; PySWEB treats all built-up land as non-agricultural and excludes mangroves from the reference group. This is a deliberate, conservative crosswalk.
+
+1. Average valid, non-water LST, NDVI and dT onto a 240 m grid. Compute `Tc_fine = LST_fine - 0.125 × dT_fine × (0.9 − NDVI_fine) × 10`.
+2. Sample the native green agricultural mask onto the fine grid by nearest neighbour. Aggregate the eligible fine LST/NDVI cells to 4800 m and apply the same expression with coarse dT. As in the OpenET structure, fine means can mix neighbouring land covers at class boundaries; exclusion from the reference mask does not make mixed fine cells pure crop observations.
+3. Apply a 3 × 3 coarse mean, fill remaining holes with the mean of all valid fine Tc cells over the supplied scene extent, then apply another 3 × 3 coarse mean. The scene fallback can include non-agricultural land and is explicitly tracked.
+4. At 90 m, use local fine Tc for green agricultural cells (fine NDVI ≥ 0.4) and for non-agricultural land. Dry agricultural cells use the regional estimate. Smooth the combined field with a square **radius of 5 pixels (11 × 11)**, then sample back to native Landsat pixels. Negative NDVI uses native LST, except anomalous shrub/bare/built classes use 250 K. Water remains missing.
+
+The projected grid origin is `(15, 15)` at every scale, matching OpenET's transforms; changing a crop window does not reset the grid origin. CRS must be north-up and in metres. Continuous reductions use full area-weighted raster means, not Earth Engine's `bestEffort` pyramid reductions. Categorical fields and final grid sampling use nearest neighbour. PySWEB averages dT at each working scale and retains its existing meteorology, Landsat temperature scaling and QA. It does not add OpenET's separate emissivity correction or its NDWI/JRC surface-water logic. These differences prevent claims of numerical parity with OpenET.
+
+Reference extent is a spatial-support decision, not a multiple of AOI area. Two coarse means can depend on a 5 × 5 set of original coarse cells (24 km across), plus final smoothing; the fallback depends on the whole supplied scene. Native workflows default to 14.4 km acquisition padding, but this is not a proof of sufficient coverage for every scene. Configured retrieval protocols continue to compare two explicit padding windows and reject fallback-dependent or unstable estimates using their stated thresholds. `cold_boundary_fallback_fraction` now follows the actual branch selection and smoothing weights. `worldcover_class` and `fano_reference_eligible` accompany scene temperatures and ET fractions in diagnostics.
+
+Removing forest references does not guarantee positive ET: the selected reference may still place a hot, dry field above the hot boundary. Independent ET evaluation and padding checks remain necessary before interpreting a changed fraction as an improvement.
+
 ## Calendars, storage and closure
 
 Effective rain uses the repository's monthly Smith calculation. Preprocessing requires every daily timestamp in each intersecting calendar month, computes monthly effective totals and distributes them proportionally to daily precipitation, then selects the requested simulation dates. Any missing pixel-day makes that pixel's effective precipitation missing for the month. This avoids interpreting an incomplete month as a complete rainfall total.

docs/migration.md

@@ -1,5 +1,13 @@
 # Migration to the maintained package
 
+## WorldCover FANO default (22 October 2026)
+
+New SSEBop retrievals use `tcold_method="worldcover_fano"`, with a required ESA WorldCover v200 land-cover raster. Preparation downloads `dem/worldcover_v200_2021.tif`; configured native runs pass it automatically. Standalone runs need `--landcover`. Existing prepared inputs can acquire the additional raster through `pysweb.io.acquisition.acquire_worldcover`, using the complete scene/reference extent and the existing Earth Engine project.
+
+This changes reference selection, fine/coarse branch selection and smoothing, so it intentionally changes temperature and ET values. Recompute upstream forcing and recalibrate for new experiments. Keep previous outputs and parameter fits together. To reproduce the previous calculation, select `tcold_method="legacy"` explicitly. The existing `tcold_fano_local_xr` and `tcold_fano_simple_xr` numerical helpers retain their legacy behavior.
+
+The new diagnostic `cold_boundary_fallback_fraction` measures scene-scalar contribution after branch selection and smoothing, rather than the former coarse-hole indicator. Scene cache keys include the selected method/settings and land-cover contents. Native acquisition padding defaults to 14400 m (previously 6000 m); existing explicit padding remains honored and must be checked. [Methods](methods.md#land-cover-aware-fano-cold-boundary) documents the WorldCover crosswalk, 2021 epoch and differences from OpenET.
+
 ## Explicit repeated-spin-up switch (21 October 2026)
 
 Repeated forcing is disabled by default. Set `[spinup] enabled = true` to enable

docs/native-grids.md

@@ -45,7 +45,7 @@ This follows Earth Engine's [export alignment guidance](https://developers.googl
 
 **Calibration:** select every complete native GSSM cell intersecting the requested chip. Inputs cover those entire footprints, including portions outside the chip. Continuous fields use area-overlap averaging; missing source support remains missing. GSSM observations are exact subsets, without interpolation. A nominal 5.1 km × 5.1 km chip contains roughly a few dozen 1 km cells; counts depend on latitude and alignment.
 
-**Simulation:** select complete pixels from the chosen native Landsat UTM scene. The exact 30 m grid is shared by SEB and SWB products. Meteorology and soil are explicitly interpolated to this grid. The processing extent is padded by `processing_buffer_m` (default 6000 m) for FANO and whole GSSM footprints, then the final SEB file is cropped without resampling. Increase the buffer if spatial FANO scales require it. Native UTM grids avoid degree-based pixel dimensions; source coverage limitations still apply.
+**Simulation:** select complete pixels from the chosen native Landsat UTM scene. The exact 30 m grid is shared by SEB and SWB products. Meteorology and soil are explicitly interpolated to this grid. The processing extent is padded by `processing_buffer_m` (default 14400 m) for FANO and whole GSSM footprints, then the final SEB file is cropped without resampling. WorldCover covers that acquisition extent and uses categorical nearest-neighbour alignment. Check padding sensitivity and fallback contributions; a fixed buffer does not guarantee adequate reference support for every scene. Native UTM grids avoid degree-based pixel dimensions; source coverage limitations still apply. See [WorldCover FANO](methods.md#land-cover-aware-fano-cold-boundary).
 
 The code rejects a source template that is not a 30 m UTM grid. This deliberate validation also rejects earlier geographic Landsat downloads. Reacquire those inputs into a new run directory; changing only the CRS label cannot correct them.
 

docs/workflows.md

@@ -51,6 +51,14 @@ Run `pysweb demo` first to generate those synthetic source files. `sm_res` is in
 
 ## Calibration and persistence
 
+### WorldCover FANO inputs
+
+SSEBop defaults to `tcold_method="worldcover_fano"`. `ssebop-prepare` downloads a native WorldCover v200 chip to `dem/worldcover_v200_2021.tif`, covering the same padded acquisition extent as Landsat. The native workflow passes this path automatically. For standalone execution, supply `--landcover <prepared>/dem/worldcover_v200_2021.tif`. Existing inputs without land cover require this additional preparation before retrieval.
+
+An optional `[ssebop] landcover = "path/to/worldcover_v200_2021.tif"` selects an existing raster (resolved relative to the TOML). In the configured scene workflow, `forcing.retrieval.landcover_file` is also supported; `[ssebop].landcover` takes precedence. Without either override, the prepared file is placed beside `forcing.retrieval.dem_file`. Preparation preserves explicitly supplied existing land-cover rasters. They must contain WorldCover v200 class codes and cover all valid pixels of every downloaded scene; the loader reads only the required spatial window. Acquisition follows the existing bounded direct-download contract; requests exceeding Earth Engine's download limits must be split into smaller study regions.
+
+Use `tcold_method="legacy"` explicitly when reproducing previous land-cover-unaware retrievals. `tcold_smooth_scale_m` applies only to that method. The new method retains the existing coefficient, NDVI threshold, fine-scale and coarse-scale options, and adds `tcold_smooth_resolution_m=90` and `tcold_smooth_radius_pixels=5`. Method/settings, code and land-cover contents participate in retrieval cache identity. Recompute SSEBop forcing and recalibrate after changing methods; do not treat old calibrated parameters as a fit to the new forcing. See [the algorithm and crosswalk](methods.md#land-cover-aware-fano-cold-boundary).
+
 Calibration is opt-in. Set its dates, random seed and spin-up exclusion in `[calibration]`. A non-converged optimizer fails unless `allow_unconverged=true` is explicit; the saved manifest records convergence, objective, sample count and optimizer status. The native workflow uses GSSM cell-day errors. Individual legacy calibration calls retain optional `sm_res` block coarsening and `objective="domain_mean"` for deliberate reproduction; these options are not the native workflow defaults.
 
 The workflow automatically supplies the calibrated CSV to the run. Its checksum is included in run provenance. The calibrated parameters are `diff_factor` and `root_beta`. Bounds and units are available in CLI help. Defaults are research starting values, and independent validation remains necessary.

examples/custom-soil.toml

@@ -7,7 +7,7 @@ end = "2024-01-31"
 extent = [148.0, -35.02, 148.02, -35.0]
 grid_mode = "native"
 sensors = ["LC08", "LC09"] # Add "LT05", "LE07" explicitly when needed.
-processing_buffer_m = 6000
+processing_buffer_m = 14400
 workers = 1
 gee_project = "your-gee-project"
 state_timing = "start"

examples/workflow.toml

@@ -12,7 +12,7 @@ center = [147.25, -35.05] # Default 5100 m square: 170 x 170 native 30 m pixels.
 # Repeated spin-up is disabled by default; enable it with [spinup] enabled = true.
 grid_mode = "native"
 sensors = ["LC08", "LC09"] # Add "LT05", "LE07" explicitly when needed.
-processing_buffer_m = 6000
+processing_buffer_m = 14400 # FANO regional neighbourhood; verify padding sensitivity.
 workers = 1
 gee_project = "your-gee-project"
 state_timing = "start"
@@ -33,6 +33,8 @@ diff_factor = 1000.0
 root_beta = 0.96
 
 [ssebop]
+tcold_method = "worldcover_fano"
+# landcover = "inputs/dem/worldcover_v200_2021.tif" # Optional override; prepare downloads it.
 landsat_qa = "required"
 lst_band = "auto"
 red_band = "SR_B4"

pysweb/io/acquisition.py

@@ -3,7 +3,7 @@
 Objective: Acquire exact native Landsat scenes and compact native meteorology with restart receipts.
 Author: Yi Yu (with assistance from Codex)
 Created: 2026-10-17
-Last updated: 2026-10-17
+Last updated: 2026-10-22
 Inputs: Audited scene manifest, geographic envelope and an existing Earth Engine project.
 Outputs: Original-grid rasters, explicit L7 optical masks and source metadata ledgers.
 Usage: from pysweb.io.acquisition import acquire_landsat_manifest, acquire_monthly_meteorology
@@ -162,6 +162,11 @@ def prepare_native_sources(manifest_file, protocol, *, gee_project, job_director
     catalogue=acquire_landsat_manifest(manifest_file,protocol['acquisition_extent'],gee_project=gee_project,sample=sample)
     if sample:
         return catalogue
+    landcover_file = protocol.get('landcover_file',
+        str(Path(protocol['dem_file']).parent/'worldcover_v200_2021.tif'))
+    if 'landcover_file' not in protocol or not Path(landcover_file).exists():
+        extent = worldcover_extent(catalogue['files'], protocol['acquisition_extent'])
+        acquire_worldcover(landcover_file, extent, gee_project=gee_project)
     raw_met=Path(protocol['meteorology_raw_dir'])
     acquire_monthly_meteorology(raw_met,protocol['meteorology_extent'],protocol['climatology_start'],
         protocol['meteorology_end'],gee_project=gee_project)
@@ -177,4 +182,40 @@ def prepare_native_sources(manifest_file, protocol, *, gee_project, job_director
     unpack_monthly_meteorology(raw_met,daily)
     stack_era5land_daily_inputs(daily,protocol['dem_file'],protocol['climatology_start'],
         protocol['meteorology_end'],protocol['meteorology_stack_dir'])
-    return dict(landsat=catalogue,meteorology=monthly,dem=protocol['dem_file'])
+    return dict(landsat=catalogue,meteorology=monthly,dem=protocol['dem_file'],landcover=landcover_file)
+
+
+def worldcover_extent(scene_files, requested_extent):
+    """Cover actual exported scene cells, including UTM bounding-box corners.
+
+    Native exports snap the requested geographic box to a projected rectangle;
+    its corners can extend outside that original box. Inspect metadata only.
+    """
+    import rasterio
+    from rasterio.warp import transform_bounds
+    left, bottom, right, top = requested_extent
+    for path in scene_files:
+        with rasterio.open(path) as source:
+            bounds = source.bounds
+            rx, ry = map(abs, source.res)
+            lo, bo, hi, to = transform_bounds(source.crs, 'EPSG:4326',
+                bounds.left-rx, bounds.bottom-ry, bounds.right+rx, bounds.top+ry, densify_pts=65)
+        left, bottom, right, top = min(left, lo), min(bottom, bo), max(right, hi), max(top, to)
+    return [left, bottom, right, top]
+
+
+def acquire_worldcover(output_file, extent, *, gee_project):
+    """Acquire a native categorical WorldCover v200 chip with a hash receipt.
+
+    This is a static 2021 classification, not a scene-year land-cover product.
+    The extent must cover the complete Landsat reference region, not just AOI.
+    """
+    import ee
+    from pysweb.ssebop.landcover import WORLDCOVER_ASSET
+    ee.Initialize(project=gee_project)
+    ee.data.setDeadline(180000)
+    Path(output_file).parent.mkdir(parents=True, exist_ok=True)
+    image = ee.ImageCollection(WORLDCOVER_ASSET).first().select('Map')
+    return _download(image.toUint8(), output_file, ['Map'], extent,
+        dict(collection=WORLDCOVER_ASSET, year=2021, band='Map',
+             citation='Zanaga et al. (2022), doi:10.5281/zenodo.7254221', license='CC-BY-4.0'))