diff --git a/README.md b/README.md index 57fc703..c5d8090 100755 --- a/README.md +++ b/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 | diff --git a/docs/methods.md b/docs/methods.md index a344e5a..ef7398d 100644 --- a/docs/methods.md +++ b/docs/methods.md @@ -2,7 +2,7 @@ ## 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). @@ -10,6 +10,33 @@ Separate `actual_evaporation`, `actual_transpiration`, `unmet_evaporation` and ` 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. diff --git a/docs/migration.md b/docs/migration.md index 980ee34..8a6d676 100644 --- a/docs/migration.md +++ b/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 diff --git a/docs/native-grids.md b/docs/native-grids.md index d4c3a17..6ed3248 100644 --- a/docs/native-grids.md +++ b/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. diff --git a/docs/workflows.md b/docs/workflows.md index dd5b4ab..b6d3eca 100644 --- a/docs/workflows.md +++ b/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 /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. diff --git a/examples/custom-soil.toml b/examples/custom-soil.toml index 47e3425..d34c511 100644 --- a/examples/custom-soil.toml +++ b/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" diff --git a/examples/workflow.toml b/examples/workflow.toml index 148fa70..9fa0aa3 100644 --- a/examples/workflow.toml +++ b/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" diff --git a/pysweb/io/acquisition.py b/pysweb/io/acquisition.py index 834d818..7733adc 100644 --- a/pysweb/io/acquisition.py +++ b/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')) diff --git a/pysweb/ssebop/api.py b/pysweb/ssebop/api.py index 4b2d6b5..0af1015 100644 --- a/pysweb/ssebop/api.py +++ b/pysweb/ssebop/api.py @@ -4,7 +4,7 @@ Objective: Provide package-owned SSEBop input-preparation and model-run APIs. Author: Yi Yu (with assistance from Codex) Created: 2026-04-17 -Last updated: 2026-10-17 +Last updated: 2026-10-22 Inputs: API parameters, optional Landsat config templates, local Landsat GeoTIFFs, meteorology NetCDFs, and DEM rasters. Outputs: Prepared inputs plus SSEBop ET GeoTIFF and NetCDF products in the requested output directory. Usage: import pysweb.ssebop.api @@ -103,6 +103,7 @@ def prepare_inputs( dem_source: str = "nasadem", gee_config_template: str | None = None, sensors: tuple[str, ...] | list[str] = ("LC08", "LC09"), + landcover_file: str | None = None, ) -> None: gee_project = gee_project.strip() if not gee_project: @@ -126,6 +127,13 @@ def prepare_inputs( sensors = sensors, ) + from pysweb.io.acquisition import acquire_worldcover, worldcover_extent + from pysweb.ssebop.landcover import WORLDCOVER_FILENAME + if landcover_file is None or not Path(landcover_file).exists(): + coverage = worldcover_extent(Path(landsat_dir).glob('*.tif'), extent) + acquire_worldcover(landcover_file or str(Path(dem_dir)/WORLDCOVER_FILENAME), + coverage, gee_project=gee_project) + month_start, month_end = month_bounds(start_date, end_date) era5land_download.download_era5land_daily( start_date = month_start.strftime("%Y-%m-%d"), @@ -513,6 +521,8 @@ def _process_landsat_scene_worker(tif_path: str) -> Tuple[np.datetime64, str, st water_mask = _SCENE_WORKER_CONTEXT["water_mask"], # type: ignore[arg-type] dt_clim = _SCENE_WORKER_CONTEXT["dt_clim"], # type: ignore[arg-type] tcold_config = _SCENE_WORKER_CONTEXT["tcold_config"], # type: ignore[arg-type] + tcold_method = str(_SCENE_WORKER_CONTEXT["tcold_method"]), + landcover = _SCENE_WORKER_CONTEXT["landcover"], template_crs = _SCENE_WORKER_CONTEXT["template_crs"], # type: ignore[arg-type] etf_dir = str(_SCENE_WORKER_CONTEXT["etf_dir"]), ndvi_dir = str(_SCENE_WORKER_CONTEXT["ndvi_dir"]), @@ -539,7 +549,11 @@ def process_landsat_scene( output_template=None, meteorology_context=None, padding_protocol=None, + tcold_method: str = "worldcover_fano", + landcover: xr.DataArray | str | None = None, ) -> Tuple[np.datetime64, str, str]: + from pysweb.ssebop.fano import validate_method, worldcover_fano, FANO_VERSION + validate_method(tcold_method, landcover) date_str = parse_date_from_filename(tif_path) doy = datetime.fromisoformat(date_str).timetuple().tm_yday bands = read_geotiff_bands(tif_path) @@ -605,12 +619,21 @@ def process_landsat_scene( dt = reproject_match(dt_clim.sel(dayofyear=doy), lst, resampling="bilinear") if tcold_config is None: tcold_config = LocalFanoConfig() - tcold = tcold_fano_simple_xr(lst, ndvi, dt, config=tcold_config) + fano = None + if tcold_method == 'worldcover_fano': + if isinstance(landcover, (str, Path)): + landcover = load_worldcover_landcover(landcover, template=lst) + fano = worldcover_fano(lst, ndvi, dt, landcover, config=tcold_config) + tcold = fano.tcold_K + else: + tcold = tcold_fano_simple_xr(lst, ndvi, dt, config=tcold_config) etf = et_fraction_xr(lst, tcold, dt).rename("etf") from pysweb.io.provenance import atomic_netcdf fields = {"lst_K": lst, "tcold_K": tcold, "dt_K": dt, "ndvi_observed": ndvi, "etf_observed": etf, "scene_valid": np.isfinite(etf).astype('uint8')} + if fano is not None: + fields.update({name: value for name, value in fano.data_vars.items() if name != 'tcold_K'}) for name in ("QA_PIXEL", "QA_RADSAT"): if name in bands: fields[name] = bands[name] @@ -630,11 +653,12 @@ def process_landsat_scene( if output_template is None: raise ValueError('Padding screening requires the actual preparation template') from pysweb.ssebop.scene_protocol import screened_fraction - accepted, diagnostics, padding_gate = screened_fraction(lst,ndvi,dt,output_template,tcold_config,padding_protocol) + accepted, diagnostics, padding_gate = screened_fraction(lst,ndvi,dt,output_template,tcold_config,padding_protocol, + landcover=landcover, tcold_method=tcold_method) fields['etf_observed'] = accepted fields.update(diagnostics) if output_template is not None: - flags = {'QA_PIXEL','QA_RADSAT','scan_gap_candidate'} + flags = {'QA_PIXEL','QA_RADSAT','scan_gap_candidate','worldcover_class','fano_reference_eligible'} fields = {name:align_to_grid(field,output_template,resampling='nearest' if name in flags or name.endswith( ('_provenance','_method','_window','_common_count')) else 'bilinear') for name,field in fields.items()} etf = fields['etf_observed'] @@ -653,6 +677,14 @@ def process_landsat_scene( auxiliary = xr.Dataset(fields) auxiliary.attrs.update(source_scene=str(Path(tif_path).resolve()), source_crs=str(lst.rio.crs), source_transform=list(lst.rio.transform())[:6]) + auxiliary.attrs.update(tcold_method=tcold_method, + tcold_algorithm=FANO_VERSION if fano is not None else 'legacy-local-fano') + if fano is not None: + auxiliary.attrs.update(fano.attrs) + if padding_protocol: + # The retained temperature uses the primary window, so a scalar + # computed on the full download would misdescribe its fallback. + auxiliary.attrs.pop('tcold_scene_fallback_K', None) if padding_gate is not None: import json auxiliary.attrs['padding_gate_json'] = json.dumps(padding_gate) @@ -781,6 +813,9 @@ def run_ssebop_workflow( ea_var = _cfg_value(cfg, args, "ea_var") dem_path = _cfg_value(cfg, args, "dem") landcover_path = _cfg_value(cfg, args, "landcover") + tcold_method = _cfg_value(cfg, args, "tcold_method", default="worldcover_fano") + from pysweb.ssebop.fano import validate_method, FANO_VERSION + validate_method(tcold_method, landcover_path) output_dir = _cfg_value(cfg, args, "output_dir") max_gap_days = int(_cfg_value(cfg, args, "max_gap_days", default=32)) apply_water_mask = bool(_cfg_value(cfg, args, "apply_water_mask", default=False)) @@ -816,6 +851,8 @@ def run_ssebop_workflow( fine_scale_m=tcold_fine_scale_m, coarse_scale_m=tcold_coarse_scale_m, smooth_scale_m=tcold_smooth_scale_m, + smooth_resolution_m=float(_cfg_value(cfg, args, "tcold_smooth_resolution_m", default=90.0)), + smooth_radius_pixels=int(_cfg_value(cfg, args, "tcold_smooth_radius_pixels", default=5)), ) if workers < 1: @@ -862,8 +899,8 @@ def scene_rank(path): template = next(iter(sample_bands.values())) buffer = 2 * float(template.rio.resolution()[0]) - if landcover_path: - landcover = load_worldcover_landcover(landcover_path) + if landcover_path and apply_water_mask: + landcover = load_worldcover_landcover(landcover_path, template=template) landcover = reproject_match_crop_first(landcover, template, resampling="nearest", buffer=buffer) _, _, water_mask = worldcover_masks(landcover) else: @@ -973,6 +1010,8 @@ def scene_rank(path): water_mask = water_mask, dt_clim = dt_clim, tcold_config = tcold_config, + tcold_method = tcold_method, + landcover = landcover_path, template_crs = template_crs, etf_dir = etf_dir, ndvi_dir = ndvi_dir, @@ -990,6 +1029,8 @@ def scene_rank(path): "water_mask": water_mask, "dt_clim": dt_clim, "tcold_config": tcold_config, + "tcold_method": tcold_method, + "landcover": landcover_path, "template_crs": template_crs, "etf_dir": etf_dir, "ndvi_dir": ndvi_dir, @@ -1115,6 +1156,8 @@ def scene_rank(path): "f_t": transpiration_fraction_daily, } ) + result.attrs.update(tcold_method=tcold_method, + tcold_algorithm=FANO_VERSION if tcold_method == 'worldcover_fano' else 'legacy-local-fano') if reconstruction_flags is not None: for name, value in reconstruction_flags.data_vars.items(): result[name] = value @@ -1127,7 +1170,7 @@ def scene_rank(path): if not auxiliary_path.exists(): raise ValueError(f"Missing SSEBop scene auxiliaries: {auxiliary_path}") with xr.open_dataset(auxiliary_path, decode_coords='all') as source: - fields = {name: align_to_grid(da, template, resampling='nearest' if name in {'QA_PIXEL','QA_RADSAT','scene_valid','original_support','scan_gap_candidate'} or name.endswith(('_provenance','_method','_window','_common_count')) else 'bilinear') + fields = {name: align_to_grid(da, template, resampling='nearest' if name in {'QA_PIXEL','QA_RADSAT','scene_valid','original_support','scan_gap_candidate','worldcover_class','fano_reference_eligible'} or name.endswith(('_provenance','_method','_window','_common_count')) else 'bilinear') for name, da in source.data_vars.items()} scene_records.append(xr.Dataset(fields).expand_dims(scene_time=[ts])) if scene_records: diff --git a/pysweb/ssebop/cli_prepare.py b/pysweb/ssebop/cli_prepare.py index 2aa04ef..2d5ea44 100644 --- a/pysweb/ssebop/cli_prepare.py +++ b/pysweb/ssebop/cli_prepare.py @@ -1,12 +1,12 @@ #!/usr/bin/env python3 """ Script: cli_prepare.py -Objective: Prepare Landsat, ERA5-Land, DEM, and stacked meteorology inputs for the first SSEBop step. +Objective: Prepare Landsat, WorldCover, ERA5-Land, DEM and stacked meteorology for SSEBop. Author: Yi Yu (with assistance from Codex) Created: 2026-02-17 -Last updated: 2026-10-17 +Last updated: 2026-10-22 Inputs: CLI arguments (--date-range, --extent, --met-source, --gee-project, --out-dir). -Outputs: Package-managed Landsat, NASADEM, and meteorology inputs for the requested SSEBop run. +Outputs: Package-managed Landsat, WorldCover, NASADEM and meteorology for the requested SSEBop run. Usage: pysweb ssebop-prepare --help Dependencies: Python standard library; pysweb """ @@ -24,7 +24,7 @@ def build_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser( description=( - "Prepare Landsat, NASADEM, ERA5-Land daily downloads, and " + "Prepare Landsat, WorldCover v200, NASADEM, ERA5-Land daily downloads, and " "stacked ERA5-Land meteorology inputs for SSEBop." ) ) diff --git a/pysweb/ssebop/cli_run.py b/pysweb/ssebop/cli_run.py index e24f68b..4b61626 100644 --- a/pysweb/ssebop/cli_run.py +++ b/pysweb/ssebop/cli_run.py @@ -4,7 +4,7 @@ Objective: Provide a thin CLI wrapper around the package-owned SSEBop run API. Author: Yi Yu (with assistance from Codex) Created: 2026-02-17 -Last updated: 2026-10-17 +Last updated: 2026-10-22 Inputs: YAML config/CLI options, Landsat GeoTIFFs, meteorology NetCDF files, DEM, landcover raster. Outputs: Daily SSEBop ET NetCDF outputs and optional gap-filled ETf diagnostics in output directory. Usage: pysweb ssebop-run --help @@ -50,7 +50,8 @@ def build_parser() -> argparse.ArgumentParser: parser.add_argument("--ea", default=None) parser.add_argument("--ea-var", default=None) parser.add_argument("--dem", default=None) - parser.add_argument("--landcover", default=None) + parser.add_argument("--landcover", default=None, help="WorldCover v200 GeoTIFF covering the padded scene; required for worldcover_fano.") + parser.add_argument("--tcold-method", choices=["worldcover_fano", "legacy"], default="worldcover_fano") parser.add_argument("--output-dir", default=None) parser.add_argument("--max-gap-days", type=int, default=32) parser.add_argument("--development-end", help="Last development day; temporal reconstruction never crosses this boundary.") @@ -77,7 +78,10 @@ def build_parser() -> argparse.ArgumentParser: parser.add_argument("--tcold-anchor-ndvi-threshold", type=float, default=0.4) parser.add_argument("--tcold-fine-scale-m", type=float, default=240.0) parser.add_argument("--tcold-coarse-scale-m", type=float, default=4800.0) - parser.add_argument("--tcold-smooth-scale-m", type=float, default=240.0) + parser.add_argument("--tcold-smooth-scale-m", type=float, default=240.0, help="Legacy method only.") + parser.add_argument("--tcold-smooth-resolution-m", type=float, default=90.0) + parser.add_argument("--tcold-smooth-radius-pixels", type=int, default=5, + help="WorldCover FANO square radius: 5 means an 11 by 11 window.") parser.add_argument("--workers", type=int, default=1) return parser diff --git a/pysweb/ssebop/core.py b/pysweb/ssebop/core.py index 3a8d686..0cd45a8 100644 --- a/pysweb/ssebop/core.py +++ b/pysweb/ssebop/core.py @@ -4,7 +4,7 @@ Objective: Provide package-native SSEBop physics helpers for dT, ET fraction, Tcold estimation, climatology, and daily ET calculations. Author: Yi Yu (with assistance from Codex) Created: 2026-04-17 -Last updated: 2026-10-17 +Last updated: 2026-10-22 Inputs: xarray temperature and geospatial arrays plus pandas ET series. Outputs: SSEBop helper DataArray and Series results. Usage: import pysweb.ssebop.core @@ -85,6 +85,10 @@ class LocalFanoConfig: fine_scale_m: float = 240.0 coarse_scale_m: float = 4800.0 smooth_scale_m: float = 240.0 + # WorldCover FANO uses a separate 90 m grid and square radius-five mean. + # smooth_scale_m is retained solely for legacy retrievals. + smooth_resolution_m: float = 90.0 + smooth_radius_pixels: int = 5 @dataclass(init=False) diff --git a/pysweb/ssebop/fano.py b/pysweb/ssebop/fano.py new file mode 100644 index 0000000..2b5cca9 --- /dev/null +++ b/pysweb/ssebop/fano.py @@ -0,0 +1,167 @@ +""" +Script: fano.py +Objective: Estimate land-cover-aware FANO cold boundaries on fixed projected grids. +Author: Yi Yu (with assistance from Codex) +Created: 2026-10-22 +Last updated: 2026-10-22 +Inputs: Aligned Landsat temperature, NDVI, dT and ESA WorldCover classifications. +Outputs: Cold temperature, scene-fallback contribution and reference diagnostics. +Usage: from pysweb.ssebop.fano import worldcover_fano +Dependencies: numpy, scipy, xarray, rasterio, rioxarray, pyproj, pysweb +""" +import math + +import numpy as np +import xarray as xr +from pyproj import CRS +from rasterio import Affine +from rasterio.enums import Resampling +from rasterio.warp import reproject +from scipy.ndimage import uniform_filter + +from pysweb.ssebop.core import LocalFanoConfig +from pysweb.ssebop.landcover import WORLDCOVER_CLASSES, worldcover_masks + + +FANO_VERSION = "worldcover-fano-v1" +TCOLD_METHODS = ("worldcover_fano", "legacy") + + +def validate_method(method, landcover): + if method not in TCOLD_METHODS: + raise ValueError(f"tcold_method must be one of {TCOLD_METHODS}") + if method == "worldcover_fano" and landcover is None: + raise ValueError("worldcover_fano requires ESA WorldCover landcover; run ssebop-prepare " + "or provide landcover/forcing.retrieval.landcover_file. " + "Use tcold_method='legacy' only to reproduce historical retrievals.") + + +def config_from_options(options): + """Use the same tuning contract in direct and configured scene workflows.""" + from dataclasses import fields + return LocalFanoConfig(**{field.name: options['tcold_'+field.name] + for field in fields(LocalFanoConfig) + if 'tcold_'+field.name in options}) + + +def _grid(template, scale): + """Anchor grids globally, so cropping a scene cannot shift the coarse cells.""" + left, bottom, right, top = template.rio.bounds() + left = 15 + math.floor((left-15)/scale)*scale + top = 15 + math.ceil((top-15)/scale)*scale + width, height = math.ceil((right-left)/scale), math.ceil((top-bottom)/scale) + transform = Affine(scale, 0, left, 0, -scale, top) + return xr.DataArray(np.full((height, width), np.nan), dims=('y', 'x'), + coords={'x': left+(np.arange(width)+.5)*scale, + 'y': top-(np.arange(height)+.5)*scale}).rio.write_crs( + template.rio.crs).rio.write_transform(transform) + + +def _warp(source, target, method='nearest'): + output = np.full(target.shape, np.nan, dtype='float64') + reproject(np.asarray(source, dtype='float64'), output, + src_transform=source.rio.transform(), src_crs=source.rio.crs, src_nodata=np.nan, + dst_transform=target.rio.transform(), dst_crs=target.rio.crs, dst_nodata=np.nan, + resampling=Resampling[method], num_threads=1) + return target.copy(data=output) + + +def _mean(data, radius): + """Square focal mean, ignoring missing cells without extending edge values.""" + values = np.asarray(data, dtype='float64') + finite = np.isfinite(values) + size = 2*radius+1 + numerator = uniform_filter(np.where(finite, values, 0), size=size, mode='constant') + denominator = uniform_filter(finite.astype(float), size=size, mode='constant') + result = np.full(values.shape, np.nan) + np.divide(numerator, denominator, out=result, where=denominator > 1e-12) + return data.copy(data=result) + + +def worldcover_fano(lst, ndvi, dt, landcover, *, config=None): + """OpenET-style FANO with a documented WorldCover crosswalk. + + Fine/coarse area means use all available raster pixels, rather than Earth + Engine's bestEffort pyramid reductions. Categorical masks use nearest + sampling. See docs/methods.md for crosswalk, water and scene-extent limits. + """ + validate_method('worldcover_fano', landcover) + config = config or LocalFanoConfig() + scales = (config.fine_scale_m, config.coarse_scale_m, config.smooth_resolution_m) + if (not all(np.isfinite(s) and s > 0 for s in scales) + or config.coarse_scale_m < config.fine_scale_m): + raise ValueError('FANO scales must be positive with coarse >= fine') + radius = config.smooth_radius_pixels + if not isinstance(radius, (int, np.integer)) or radius < 0: + raise ValueError('smooth_radius_pixels must be a nonnegative integer') + if (not np.isfinite(config.dt_coeff) or config.dt_coeff <= 0 + or not 0 <= config.anchor_ndvi_threshold <= config.high_ndvi_threshold <= 1): + raise ValueError('Invalid FANO coefficient or NDVI thresholds') + for field in (lst, ndvi, dt): + if field.dims != ('y', 'x') or field.rio.crs != lst.rio.crs: + raise ValueError('FANO requires aligned two-dimensional y/x inputs') + xr.align(lst, field, join='exact') + crs = CRS(lst.rio.crs) + transform = lst.rio.transform() + if (not crs.is_projected or any(axis.unit_conversion_factor != 1 for axis in crs.axis_info) + or transform.b or transform.d or transform.a <= 0 or transform.e >= 0): + raise ValueError('FANO requires a north-up projected grid in metres') + if landcover.dims != ('y', 'x') or landcover.rio.crs is None: + raise ValueError('WorldCover requires a georeferenced two-dimensional raster') + cover = _warp(landcover, lst) + valid = np.isfinite(lst) & np.isfinite(ndvi) & np.isfinite(dt) & (dt > 0) + valid = valid & (ndvi >= -1) & (ndvi <= 1) + # Unknown or absent cover is a missing input, never an implicit non-ag class. + if bool((valid & ~cover.isin(WORLDCOVER_CLASSES)).any()): + raise ValueError('WorldCover does not cover all valid scene pixels with recognised classes') + ag, anomalous, water = worldcover_masks(cover) + land = valid & ~water.astype(bool) + fine, coarse, smooth = (_grid(lst, scale) for scale in scales) + fine_ndvi = _warp(ndvi.where(land), fine, 'average') + fine_lst = _warp(lst.where(land), fine, 'average') + fine_dt = _warp(dt.where(land), fine, 'average') + fine_tc = fine_lst-config.dt_coeff*fine_dt*(config.high_ndvi_threshold-fine_ndvi)*10 + + # OpenET masks fine-resolution fields with native green agricultural + # pixels before the coarse reduction. Nearest sampling retains class IDs. + anchor = land & ag.astype(bool) & (ndvi >= config.anchor_ndvi_threshold) + fine_anchor = _warp(anchor.astype(float), fine) == 1 + coarse_ndvi = _warp(fine_ndvi.where(fine_anchor), coarse, 'average') + coarse_lst = _warp(fine_lst.where(fine_anchor), coarse, 'average') + coarse_dt = _warp(dt.where(land), coarse, 'average') + coarse_tc = coarse_lst-config.dt_coeff*coarse_dt*(config.high_ndvi_threshold-coarse_ndvi)*10 + regional = _mean(coarse_tc, 1) + available = np.asarray(fine_tc)[np.isfinite(fine_tc)] + scene_tc = float(available.mean()) if available.size else np.nan + missing = regional.isnull() + regional = _mean(regional.fillna(scene_tc), 1) + # Propagate the fraction of the temperature estimate supplied by the + # scene scalar through exactly the same second and final means. + regional_fallback = _mean(missing.astype(float).where(np.isfinite(regional)), 1) + ag_smooth = _warp(ag.astype(float), smooth) == 1 + green_smooth = _warp(fine_ndvi, smooth) >= config.anchor_ndvi_threshold + use_regional = ag_smooth & ~green_smooth + local = _warp(fine_tc, smooth) + layered = xr.where(use_regional, _warp(regional, smooth), local) + weight = xr.where(use_regional, _warp(regional_fallback, smooth), 0).where(np.isfinite(layered)) + cold = _warp(_mean(layered, radius), lst) + fallback = _warp(_mean(weight, radius), lst).clip(0, 1) + # Mirror FANO's negative-NDVI special case for shrub/bare/built classes. + # Water is excluded from terrestrial ET, including when optional masking + # is disabled by the caller; it cannot enter reference statistics. + cold = xr.where(ndvi < 0, xr.where(anomalous.astype(bool), 250., lst), cold).where(land) + fallback = xr.where(ndvi < 0, 0., fallback).where(land) + result = xr.Dataset(dict(tcold_K=cold, cold_boundary_fallback_fraction=fallback, + fano_reference_eligible=anchor.astype('uint8'), worldcover_class=cover)) + result.attrs.update(tcold_method='worldcover_fano', tcold_algorithm=FANO_VERSION, + landcover_crosswalk='reference=30,40,90; anomalous=20,50,60; water=80', + landcover_epoch='2021 (ESA WorldCover v200); static for all scene dates', + tcold_scene_fallback_K=scene_tc) + result.tcold_K.attrs.update(units='K', long_name='FANO cold boundary temperature') + result.cold_boundary_fallback_fraction.attrs.update(units='1', + long_name='Scene-scalar contribution to cold boundary after branch selection and smoothing') + result.fano_reference_eligible.attrs.update(units='1', + long_name='Valid native pixel with eligible WorldCover class and NDVI above anchor threshold') + result.worldcover_class.attrs.update(source='ESA/WorldCover/v200', epoch=2021, + spatial_resampling='nearest', long_name='ESA WorldCover class') + return result.rio.write_crs(lst.rio.crs).rio.write_transform(transform) diff --git a/pysweb/ssebop/landcover.py b/pysweb/ssebop/landcover.py index 156cdef..030f5d9 100644 --- a/pysweb/ssebop/landcover.py +++ b/pysweb/ssebop/landcover.py @@ -4,7 +4,7 @@ Objective: Provide package-native SSEBop landcover loaders and masks for ESA WorldCover inputs. Author: Yi Yu (with assistance from Codex) Created: 2026-04-17 -Last updated: 2026-04-17 +Last updated: 2026-10-22 Inputs: WorldCover raster paths and xarray landcover grids. Outputs: Landcover DataArray values and class masks. Usage: import pysweb.ssebop.landcover @@ -20,6 +20,10 @@ __all__ = ["load_worldcover_landcover", "worldcover_masks"] +WORLDCOVER_CLASSES = (10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100) +WORLDCOVER_ASSET = 'ESA/WorldCover/v200' +WORLDCOVER_FILENAME = 'worldcover_v200_2021.tif' + def _ensure_spatial_dims(data_array: xr.DataArray) -> xr.DataArray: """Ensure spatial dims are named consistently for package helpers.""" @@ -33,18 +37,32 @@ def _ensure_spatial_dims(data_array: xr.DataArray) -> xr.DataArray: def load_worldcover_landcover( path: Optional[str] = None, masked: bool = True, + template: xr.DataArray | None = None, ) -> xr.DataArray: """Load ESA WorldCover (v200) GeoTIFF and return landcover classes.""" if path is None: raise ValueError("path is required for package-level WorldCover loading") - lc = rioxarray.open_rasterio(path, masked=masked).squeeze("band", drop=True) + with rioxarray.open_rasterio(path, masked=masked) as source: + if template is not None: + from rasterio.warp import transform_bounds + bounds = transform_bounds(template.rio.crs, source.rio.crs, + *template.rio.bounds(), densify_pts=21) + rx, ry = map(abs, source.rio.resolution()) + source = source.rio.clip_box(bounds[0]-rx, bounds[1]-ry, + bounds[2]+rx, bounds[3]+ry) + lc = source.squeeze("band", drop=True).load() return _ensure_spatial_dims(lc) def worldcover_masks( landcover: xr.DataArray, ) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray]: - """Build ag, anomalous, and water masks from WorldCover classes.""" + """Build reference, anomalous and water masks (not an exact NLCD crosswalk). + + Tree cover and mangroves are excluded from agricultural references. + WorldCover does not distinguish developed open space from dense built-up + land, so all built-up pixels use the anomalous/non-agricultural branch. + """ ag_mask = landcover.isin([30, 40, 90]).astype("uint8") anomalous_mask = landcover.isin([20, 50, 60]).astype("uint8") water_mask = landcover.isin([80]).astype("uint8") diff --git a/pysweb/ssebop/scene_protocol.py b/pysweb/ssebop/scene_protocol.py index eb4b9f4..aab61a7 100644 --- a/pysweb/ssebop/scene_protocol.py +++ b/pysweb/ssebop/scene_protocol.py @@ -3,7 +3,7 @@ Objective: Retrieve padded native scenes using bounded development climatology and acceptance gates. Author: Yi Yu (with assistance from Codex) Created: 2026-10-17 -Last updated: 2026-10-17 +Last updated: 2026-10-22 Inputs: Native meteorology, DEM, Landsat fields and configured padding/climatology rules. Outputs: Calendar-day dT, local-anchor fractions and independently screened retrievals. Usage: from pysweb.ssebop.scene_protocol import calendar_day_dt, screened_fraction @@ -70,16 +70,23 @@ def _window(data, bounds, crs): return selected -def screened_fraction(lst, ndvi, dt, template, config, protocol): +def screened_fraction(lst, ndvi, dt, template, config, protocol, *, landcover=None, + tcold_method='worldcover_fano'): """Compare fixed padding windows after excluding cold-anchor fallback.""" + from pysweb.ssebop.fano import validate_method, worldcover_fano + validate_method(tcold_method, landcover) fractions, fallback, cold_fields = {}, {}, {} for label in ("primary", "comparison"): temperature = _window(lst, protocol[label+"_bounds"], protocol["crs"]) vegetation = ndvi.sel(x=temperature.x, y=temperature.y) delta = dt.sel(x=temperature.x, y=temperature.y) - cold = tcold_fano_local_xr(temperature, vegetation, delta, config=config) + if tcold_method == 'worldcover_fano': + diagnostic = worldcover_fano(temperature, vegetation, delta, landcover, config=config) + cold, weight = diagnostic.tcold_K, diagnostic.cold_boundary_fallback_fraction + else: + cold = tcold_fano_local_xr(temperature, vegetation, delta, config=config) + weight = fallback_fraction(vegetation, np.isfinite(temperature) & np.isfinite(vegetation) & np.isfinite(delta), config) fraction = et_fraction_xr(temperature, cold, delta) - weight = fallback_fraction(vegetation, np.isfinite(temperature) & np.isfinite(vegetation) & np.isfinite(delta), config) fallback[label] = align_to_grid(weight, template, resampling="bilinear") fractions[label] = align_to_grid(fraction, template, resampling="bilinear").where( fallback[label] <= protocol["maximum_fallback_fraction"]) @@ -116,10 +123,12 @@ def _retrieve_scene(path): for sub in ('etf','ndvi'): (folder/sub).mkdir(parents=True,exist_ok=True) paths = [path, context['meteorology']['dem'], *context['met_paths']] + if context['landcover']: + paths.append(context['landcover']) gapfill = context['gapfill_products'].get(str(Path(path).resolve())) if gapfill: paths.append(gapfill) - config = dict(padding=context['padding'], climatology_start=context['meteorology']['start'], + config = dict(tcold=context['tcold_options'], padding=context['padding'], climatology_start=context['meteorology']['start'], climatology_end=context['meteorology']['end'], template=context['template_identity']) manifest = input_manifest(config, paths) outputs = [folder/'etf'/f'auxiliary_{date}.nc', folder/'etf'/f'etf_{date}.tif', folder/'ndvi'/f'ndvi_{date}.tif'] @@ -128,10 +137,12 @@ def _retrieve_scene(path): if previous.get('fingerprint') == manifest['fingerprint'] and all( p.exists() and previous.get('outputs', {}).get(str(p)) == sha256(p) for p in outputs): return np.datetime64(date, 'ns'), str(outputs[1]), str(outputs[2]) - from pysweb.ssebop.core import LocalFanoConfig + from pysweb.ssebop.fano import config_from_options result = process_landsat_scene(str(path), 'auto', 'ndvi', 'SR_B4', 'SR_B5', False, None, None, context['template'].rio.crs, str(folder/'etf'), str(folder/'ndvi'), - tcold_config=LocalFanoConfig(), gapfill_product=gapfill, output_template=context['template'], + tcold_config=config_from_options(context['tcold_options']), + tcold_method=context['tcold_options']['tcold_method'], landcover=context['landcover'], + gapfill_product=gapfill, output_template=context['template'], meteorology_context=context['meteorology'], padding_protocol=context['padding']) atomic_json(receipt, {**manifest, 'outputs': {str(p): sha256(p) for p in outputs}}) return result @@ -156,7 +167,7 @@ def scene_support_calendar(fractions, support_aoi=None): def run_scene_workflow(*, date_range, scene_files, gapfill_products, dem, - output_dir, extent, protocol, reconstruction, workers=1): + output_dir, extent, protocol, reconstruction, workers=1, tcold_options=None): """Configured bounded-memory retrieval with persistent shared scene caching. The meteorology files explicitly cover the climatology and actual model @@ -171,6 +182,10 @@ def run_scene_workflow(*, date_range, scene_files, gapfill_products, dem, from pysweb.io.provenance import input_manifest, atomic_netcdf, save_output_manifest, atomic_json from pysweb.forcing import reconstruct from pysweb.ssebop.core import et_from_fraction + from pysweb.ssebop.fano import validate_method, FANO_VERSION + tcold_options = {'tcold_method': 'worldcover_fano', **(tcold_options or {})} + landcover = tcold_options.get('landcover') or protocol.get('landcover_file') + validate_method(tcold_options['tcold_method'], landcover) if workers < 1 or not scene_files: raise ValueError('Positive workers and selected scenes required') start, end = parse_date_range(date_range) @@ -192,6 +207,7 @@ def run_scene_workflow(*, date_range, scene_files, gapfill_products, dem, context = dict(template=template, template_identity=dict(crs=str(template.rio.crs), transform=list(template.rio.transform())[:6], shape=list(template.shape)), cache=str(cache), met_paths=met_paths, gapfill_products=gapfill_products, + tcold_options=tcold_options, landcover=landcover, meteorology=dict(fields={k:meteorology[k] for k in ('tmax','tmin','rs','ea')}, dem=str(dem), start=protocol['climatology_start'], end=protocol['climatology_end']), padding=protocol['padding']) @@ -246,6 +262,8 @@ def run_scene_workflow(*, date_range, scene_files, gapfill_products, dem, from pysweb.components import COMPONENT_SCHEMA result = result.rio.write_crs(template.rio.crs).rio.write_transform(template.rio.transform()) result.attrs.update(component_schema=COMPONENT_SCHEMA, selected_scenes='\n'.join(map(str,scene_files)), + tcold_method=tcold_options['tcold_method'], + tcold_algorithm=FANO_VERSION if tcold_options['tcold_method'] == 'worldcover_fano' else 'legacy-local-fano', climatology_start=protocol['climatology_start'], climatology_end=protocol['climatology_end'], climatology_calendar='month/day; Feb29 uses Feb28/Mar1', development_end=boundary, same_day_selection='configured immutable catalogue: cloud cover then scene id', @@ -268,6 +286,6 @@ def run_scene_workflow(*, date_range, scene_files, gapfill_products, dem, field.attrs.pop('_FillValue', None) encoding[name]['_FillValue'] = None atomic_netcdf(result,path,encoding=encoding) - save_output_manifest(path,input_manifest(dict(protocol=protocol,reconstruction=reconstruction), - [*scene_files,*met_paths,dem,*gapfill_products.values()])) + save_output_manifest(path,input_manifest(dict(protocol=protocol,reconstruction=reconstruction,tcold=tcold_options), + [*scene_files,*met_paths,dem,*gapfill_products.values(), *([landcover] if landcover else [])])) return path diff --git a/pysweb/workflow.py b/pysweb/workflow.py index 1045d7e..f37f44e 100644 --- a/pysweb/workflow.py +++ b/pysweb/workflow.py @@ -3,7 +3,7 @@ Objective: Run explicitly configured stages with shared paths and model settings. Author: Yi Yu (with assistance from Codex) Created: 2026-09-12 -Last updated: 2026-10-21 +Last updated: 2026-10-22 Inputs: Study TOML configuration, requested stages and optional dry-run flag. Outputs: Resolved configuration and products from the selected model stages. Usage: pysweb workflow --config examples/workflow.toml --dry-run @@ -173,7 +173,6 @@ def load_config(path): "tmin", "rs", "ea", - "landcover", } for section, allowed_keys in [ ("model", allowed_model), @@ -194,7 +193,8 @@ def load_config(path): unknown = set(config.get(section, {})) - keys if unknown: raise ValueError(f"Unknown {section} settings: {sorted(unknown)}") - for section, keys in {"landsat": ["manifest_file"], "evaluation": ["observations_file", "regimes_file"]}.items(): + for section, keys in {"landsat": ["manifest_file"], "evaluation": ["observations_file", "regimes_file"], + "ssebop": ["landcover"]}.items(): for key in keys: if key in config.get(section, {}): value = Path(config[section][key]).expanduser() @@ -203,7 +203,7 @@ def load_config(path): config["forcing"]["peer_runs"] = [str(Path(p) if Path(p).is_absolute() else (path.parent/p).resolve()) for p in config["forcing"]["peer_runs"]] retrieval = config.get('forcing', {}).get('retrieval', {}) - for key in ('meteorology_raw_dir','meteorology_stack_dir','dem_file','scene_cache_dir','template_file'): + for key in ('meteorology_raw_dir','meteorology_stack_dir','dem_file','scene_cache_dir','template_file','landcover_file'): if key in retrieval: value = Path(retrieval[key]).expanduser() retrieval[key] = str(value if value.is_absolute() else (path.parent/value).resolve()) @@ -273,6 +273,7 @@ def run_config(config, stages=STAGES, *, dry_run=False): met_raw_dir=str(prepared / "met/era5land/raw"), met_stack_dir=str(prepared / "met/era5land/stack"), dem_dir=str(prepared / "dem"), + landcover_file=config.get('ssebop', {}).get('landcover'), ) if "ssebop" in stages: options = { @@ -282,6 +283,8 @@ def run_config(config, stages=STAGES, *, dry_run=False): "ndvi_band": "ndvi", } options.update(config.get("ssebop", {})) + if options.get('tcold_method', 'worldcover_fano') == 'worldcover_fano': + options.setdefault('landcover', str(prepared/'dem/worldcover_v200_2021.tif')) pysweb.ssebop.run( date_range=period, landsat_dir=str(prepared / "landsat"), diff --git a/pysweb/workflow_native.py b/pysweb/workflow_native.py index 7c3afb0..63f2970 100644 --- a/pysweb/workflow_native.py +++ b/pysweb/workflow_native.py @@ -3,7 +3,7 @@ Objective: Orchestrate native downloads, GSSM calibration and shared Landsat-grid products. Author: Yi Yu (with assistance from Codex) Created: 2026-10-17 -Last updated: 2026-10-21 +Last updated: 2026-10-22 Inputs: Validated workflow configuration and selected execution stages. Outputs: Two final NetCDFs, coarse calibration inputs, parameters and grid checks. Usage: pysweb workflow --config examples/workflow.toml @@ -75,7 +75,7 @@ def native_plan(config, stages): raise ValueError('SSEBop climatology cannot include dates after the development boundary') if boundary is not None and boundary < start: raise ValueError('The run must include development meteorology before the holdout boundary') - if run.get('processing_buffer_m',6000) < 0: + if run.get('processing_buffer_m',14400) < 0: raise ValueError('processing_buffer_m must be non-negative') if calibration.get('sm_res') is not None: raise ValueError('Native calibration uses the GSSM affine grid; remove calibration.sm_res') @@ -179,7 +179,7 @@ def run_native(config, stages, dry_run=False): write_dataarray(data,reference) reference_grid = subset_template(open_template(reference,soil.get('reference_var','reference_ssm')),extent) extent = list(transform_bounds(reference_grid.rio.crs,'EPSG:4326',*reference_grid.rio.bounds(),densify_pts=21)) - download_extent = buffered_extent(extent,run.get('processing_buffer_m',6000)) + download_extent = buffered_extent(extent,run.get('processing_buffer_m',14400)) met = prepared/'met/era5land/stack' period = f'{start} to {end}' if 'prepare' in stages: @@ -190,13 +190,17 @@ def run_native(config, stages, dry_run=False): job_directory = os.environ.get('PBS_JOBFS') or os.environ.get('TMPDIR') if not job_directory: raise ValueError('Set TMPDIR for bounded temporary meteorology expansion') + retrieval = dict(retrieval) + if config.get('ssebop', {}).get('landcover'): + retrieval['landcover_file'] = config['ssebop']['landcover'] prepare_native_sources(config['landsat']['manifest_file'],retrieval, gee_project=run.get('gee_project',''),job_directory=job_directory) else: pysweb.ssebop.prepare_inputs(date_range=period,extent=download_extent, met_source='era5land',gee_project=run.get('gee_project',''), landsat_dir=str(prepared/'landsat'),met_raw_dir=str(prepared/'met/era5land/raw'), - met_stack_dir=str(met),dem_dir=str(prepared/'dem'),sensors=run['sensors']) + met_stack_dir=str(met),dem_dir=str(prepared/'dem'),sensors=run['sensors'], + landcover_file=config.get('ssebop',{}).get('landcover')) seb_internal = internal/f'et_daily_ssebop_{start}_{end}.nc' seb_final = products/'SSEBop_ET.nc' if 'ssebop' in stages: @@ -204,6 +208,10 @@ def run_native(config, stages, dry_run=False): 'reference_et':str(met/f'et_tall_crop_daily_{start}_{end}.nc'), 'reference_et_var':'et_tall_crop','et_reference_type':'alfalfa'} options.update(config.get('ssebop',{})) + if options.get('tcold_method', 'worldcover_fano') == 'worldcover_fano': + retrieval = config.get('forcing', {}).get('retrieval', {}) + options.setdefault('landcover', retrieval.get('landcover_file', + str(Path(retrieval.get('dem_file', prepared/'dem/nasadem.tif')).parent/'worldcover_v200_2021.tif'))) if config.get('landsat'): from pysweb.ssebop.preprocessing import select_scenes, prepare_l7 selection = config['landsat'] @@ -242,6 +250,7 @@ def run_native(config, stages, dry_run=False): run_scene_workflow(date_range=period,scene_files=options['scene_files'], gapfill_products=options.get('gapfill_products',{}),dem=retrieval['dem_file'], output_dir=internal,extent=extent,protocol=scene_protocol, + tcold_options={key:value for key,value in options.items() if key.startswith('tcold_') or key == 'landcover'}, reconstruction=config['forcing']['reconstruction'],workers=run['workers']) else: pysweb.ssebop.run(date_range=period,landsat_dir=str(prepared/'landsat'),met_dir=str(met), diff --git a/tests/ssebop/test_api_prepare_inputs.py b/tests/ssebop/test_api_prepare_inputs.py index 8c9ab4e..50bf844 100644 --- a/tests/ssebop/test_api_prepare_inputs.py +++ b/tests/ssebop/test_api_prepare_inputs.py @@ -4,7 +4,7 @@ Objective: Verify the SSEBop package API orchestrates Landsat and meteorological input preparation steps. Author: Yi Yu (with assistance from Codex) Created: 2026-04-17 -Last updated: 2026-10-17 +Last updated: 2026-10-22 Inputs: Temporary paths and monkeypatched package functions supplied by pytest. Outputs: Test assertions. Usage: python -m pytest tests/ssebop/test_api_prepare_inputs.py @@ -13,6 +13,7 @@ from pathlib import Path import json import sys +import pytest ROOT = Path(__file__).resolve().parents[2] @@ -23,6 +24,14 @@ from pysweb.ssebop.landsat import update_gee_config +@pytest.fixture(autouse=True) +def worldcover_acquisition(monkeypatch): + calls = [] + monkeypatch.setattr('pysweb.io.acquisition.acquire_worldcover', + lambda *args, **kwargs: calls.append((args, kwargs))) + return calls + + def _load_config_payload(path: Path) -> dict: payload = path.read_text(encoding="utf-8") try: @@ -32,7 +41,7 @@ def _load_config_payload(path: Path) -> dict: return yaml.safe_load(payload) or {} -def test_prepare_inputs_calls_landsat_era5land_dem_and_stack_steps(monkeypatch, tmp_path: Path): +def test_prepare_inputs_calls_landsat_era5land_dem_and_stack_steps(monkeypatch, tmp_path: Path, worldcover_acquisition): recorded = [] monkeypatch.setattr( "pysweb.met.era5land.stack.native_download_extent", @@ -69,6 +78,8 @@ def test_prepare_inputs_calls_landsat_era5land_dem_and_stack_steps(monkeypatch, gee_project="workflow-project", ) + assert worldcover_acquisition == [((str(tmp_path/'dem/worldcover_v200_2021.tif'), + [147.2, -35.1, 147.3, -35.0]), {'gee_project': 'workflow-project'})] assert recorded == [ ( "landsat", diff --git a/tests/ssebop/test_api_run.py b/tests/ssebop/test_api_run.py index 734bc04..d22ed41 100644 --- a/tests/ssebop/test_api_run.py +++ b/tests/ssebop/test_api_run.py @@ -4,7 +4,7 @@ Objective: Verify SSEBop dispatch, process fallbacks and preservation of observations during interpolation. Author: Yi Yu (with assistance from Codex) Created: 2026-04-17 -Last updated: 2026-10-17 +Last updated: 2026-10-22 Inputs: Package API calls, temporary files, and monkeypatched package functions supplied by pytest. Outputs: Test assertions. Usage: python -m pytest tests/ssebop/test_api_run.py @@ -223,6 +223,7 @@ def fail_strict_local(*args, **kwargs): monkeypatch.setattr(ssebop_api, "tcold_fano_local_xr", fail_strict_local) ssebop_api.process_landsat_scene( + tcold_method="legacy", tif_path = "Landsat_2019-01-04.tif", lst_band = "ST_B10", ndvi_band = "ndvi", @@ -331,6 +332,7 @@ def fake_process_landsat_scene(**kwargs): monkeypatch.setattr(provenance, "input_manifest", lambda config, paths: { "fingerprint": "executor-test", "software": {"code_sha256": "test"}}) ssebop_api.run_ssebop_workflow( + tcold_method="legacy", date_range="2024-01-01 to 2024-01-02", landsat_dir=str(tmp_path), dem=str(tmp_path / "dem.tif"), diff --git a/tests/ssebop/test_fano.py b/tests/ssebop/test_fano.py new file mode 100644 index 0000000..ab2998c --- /dev/null +++ b/tests/ssebop/test_fano.py @@ -0,0 +1,162 @@ +""" +Script: test_fano.py +Objective: Verify FANO reference selection, physical branches, missing inputs and grid invariance. +Author: Yi Yu (with assistance from Codex) +Created: 2026-10-22 +Last updated: 2026-10-22 +Inputs: Synthetic projected fields with independently specified temperature expectations. +Outputs: Numerical and workflow regression assertions. +Usage: python -m pytest tests/ssebop/test_fano.py +Dependencies: numpy, pytest, xarray, rasterio, pysweb +""" +from dataclasses import replace +import numpy as np +import pytest +import xarray as xr +from rasterio import Affine +from tests.io.test_native_grids import raster +from pysweb.ssebop.core import LocalFanoConfig +from pysweb.ssebop.fano import worldcover_fano, _grid, _mean, validate_method + + +def scene(ndvi=.8, cover=40, size=96): + grid = raster(np.ones((size, size)), Affine(30, 0, 15, 0, -30, 15+30*size)) + return grid*300, grid*ndvi, grid*8, grid*cover + + +@pytest.mark.parametrize('cover', [30, 40, 90]) +def test_green_reference_matches_analytic_cold_temperature(cover): + fields = scene(cover=cover) + output = worldcover_fano(*fields) + # 300 - .125 * 8 * (.9-.8) * 10 = 299 K. + np.testing.assert_allclose(output.tcold_K, 299.) + np.testing.assert_allclose(output.cold_boundary_fallback_fraction, 0.) + assert output.fano_reference_eligible.all() + + +@pytest.mark.parametrize('cover', [10, 20, 50, 60, 70, 95, 100]) +def test_nonag_uses_fine_temperature_without_claiming_reference_support(cover): + output = worldcover_fano(*scene(ndvi=.2, cover=cover)) + np.testing.assert_allclose(output.tcold_K, 293.) + assert not output.fano_reference_eligible.any() + np.testing.assert_allclose(output.cold_boundary_fallback_fraction, 0.) + + +def test_dry_crop_without_green_references_reports_scene_fallback(): + output = worldcover_fano(*scene(ndvi=.2)) + np.testing.assert_allclose(output.tcold_K, 293.) + np.testing.assert_allclose(output.cold_boundary_fallback_fraction, 1.) + + +def test_green_forest_cannot_replace_green_crop_reference(): + lst, ndvi, dt, lc = scene(size=160) + lc.values[:, 128:] = 10 + ndvi.values[:, 128:] = .9 + lst.values[:, 128:] = 270. + # Dry field surrounded by the 299 K crop reference. + ndvi.values[64:80, 64:80] = .2 + lst.values[64:80, 64:80] = 315. + cfg = LocalFanoConfig(coarse_scale_m=480, smooth_radius_pixels=0) + first = worldcover_fano(lst, ndvi, dt, lc, config=cfg) + assert first.tcold_K.values[70, 70] == pytest.approx(299.) + assert first.cold_boundary_fallback_fraction.values[70, 70] == 0 + assert not first.fano_reference_eligible.values[:, 128:].any() + lst.values[:, 128:] = 320. + second = worldcover_fano(lst, ndvi, dt, lc, config=cfg) + assert second.tcold_K.values[70, 70] == pytest.approx(299.) + + +def test_water_excluded_and_negative_ndvi_anomalous_branch(): + lst, ndvi, dt, lc = scene(ndvi=-.2, cover=60) + lc.values[:, :32] = 80 + output = worldcover_fano(lst, ndvi, dt, lc) + assert output.tcold_K.values[:, :32].shape == (96, 32) + assert np.isnan(output.tcold_K.values[:, :32]).all() + np.testing.assert_allclose(output.tcold_K.values[:, 32:], 250.) + assert not output.fano_reference_eligible.any() + + +@pytest.mark.parametrize('missing', [np.nan, 0, 42]) +def test_unknown_landcover_is_an_error(missing): + lst, ndvi, dt, lc = scene() + lc.values[20, 20] = missing + with pytest.raises(ValueError, match='recognised classes'): + worldcover_fano(lst, ndvi, dt, lc) + + +def test_missing_observations_remain_missing(): + lst, ndvi, dt, lc = scene() + lst.values[10:30, 10:30] = np.nan + result = worldcover_fano(lst, ndvi, dt, lc) + assert np.isnan(result.tcold_K.values[10:30, 10:30]).all() + dt.values[:] = 0 + assert worldcover_fano(lst, ndvi, dt, lc).tcold_K.isnull().all() + + +def test_fixed_grid_survives_shifted_crop(): + lst, ndvi, dt, lc = scene(size=160) + # Spatial gradient detects accidental anchoring at each array's top-left. + lst.values += np.arange(160)[None, :]*.05 + fields = (lst, ndvi, dt, lc) + full = worldcover_fano(*fields) + cropped = worldcover_fano(*(v.isel(x=slice(7, 151), y=slice(5, 150)) for v in fields)) + target = dict(x=slice(45, 100), y=slice(45, 100)) + expected = full.tcold_K.isel(**target) + np.testing.assert_allclose(cropped.tcold_K.sel(x=expected.x, y=expected.y), expected) + grid = _grid(cropped.tcold_K, 4800) + assert (grid.rio.transform().c-15) % 4800 == 0 + assert (grid.rio.transform().f-15) % 4800 == 0 + + +def test_radius_five_is_eleven_cells_and_ignores_missing(): + data = xr.DataArray(np.zeros((15, 15)), dims=('y', 'x')) + data.values[7, 7] = 121 + result = _mean(data, 5) + assert result.values[7, 7] == pytest.approx(1.) + assert result.values[7, 2] > 0 + assert result.values[7, 1] == 0 + + +def test_default_rejects_missing_cover_and_legacy_is_explicit(): + with pytest.raises(ValueError, match='requires ESA WorldCover'): + validate_method('worldcover_fano', None) + validate_method('legacy', None) + with pytest.raises(ValueError, match='tcold_method'): + validate_method('typo', None) + + +@pytest.mark.parametrize('changes', [dict(fine_scale_m=0), dict(coarse_scale_m=-1), + dict(smooth_radius_pixels=2.5), dict(anchor_ndvi_threshold=2), dict(dt_coeff=np.nan)]) +def test_invalid_fano_parameters(changes): + with pytest.raises(ValueError): + worldcover_fano(*scene(), config=replace(LocalFanoConfig(), **changes)) + + +def test_worldcover_download_covers_exported_cells(tmp_path): + from rasterio.warp import transform_bounds + from pysweb.io.acquisition import worldcover_extent + grid = scene()[0] + path = tmp_path/'scene.tif' + grid.rio.to_raster(path) + bounds = transform_bounds(grid.rio.crs, 'EPSG:4326', *grid.rio.bounds(), densify_pts=65) + # Deliberately request less than the exported extent. + requested = [bounds[0]+.001, bounds[1]+.001, bounds[2]-.001, bounds[3]-.001] + actual = worldcover_extent([path], requested) + assert actual[0] < bounds[0] and actual[1] < bounds[1] + assert actual[2] > bounds[2] and actual[3] > bounds[3] + + +def test_categorical_landcover_alignment_and_windowed_loading(tmp_path): + from pysweb.ssebop.landcover import load_worldcover_landcover + lst, ndvi, dt, _ = scene(size=96) + cover = raster(np.tile(np.repeat([10., 40.], 144), (288, 1)), + Affine(10, 0, 15, 0, -10, 2895)) + path = tmp_path/'worldcover.tif' + cover.rio.to_raster(path) + chip = lst.isel(x=slice(20, 80), y=slice(20, 80)) + loaded = load_worldcover_landcover(path, template=chip) + assert loaded.size < cover.size + result = worldcover_fano(chip, ndvi.sel(x=chip.x, y=chip.y), + dt.sel(x=chip.x, y=chip.y), loaded) + assert set(np.unique(result.worldcover_class)) == {10., 40.} + assert not result.fano_reference_eligible.where(result.worldcover_class == 10, 0).any() diff --git a/tests/ssebop/test_multisensor_quality.py b/tests/ssebop/test_multisensor_quality.py index e269ddc..d3a6afd 100644 --- a/tests/ssebop/test_multisensor_quality.py +++ b/tests/ssebop/test_multisensor_quality.py @@ -4,7 +4,7 @@ Objective: Check sensor-specific band mappings, unused cirrus bits and Landsat 7 dropped-pixel handling. Author: Yi Yu (with assistance from Codex) Created: 2026-10-17 -Last updated: 2026-10-17 +Last updated: 2026-10-22 Inputs: Independently specified Collection 2 bit-field fixtures. Outputs: pytest assertions. Usage: pytest tests/ssebop/test_multisensor_quality.py @@ -44,7 +44,7 @@ def test_scene_auto_bands_preserve_the_same_physical_signal(sensor, tmp_path, mo monkeypatch.setattr(api, "reproject_match", lambda source, match, **kw: source) monkeypatch.setattr(api, "tcold_fano_simple_xr", lambda lst, ndvi, dt, **kw: lst - 4.) _, etf, ndvi = api.process_landsat_scene(str(path), "auto", "ndvi", "auto", "auto", False, None, dt, - template.rio.crs, str(tmp_path), str(tmp_path)) + template.rio.crs, str(tmp_path), str(tmp_path), tcold_method="legacy") with rasterio.open(etf) as src: np.testing.assert_allclose(src.read(1), .5) with rasterio.open(ndvi) as src: np.testing.assert_allclose(src.read(1), .055/.205, rtol=1e-6) diff --git a/tests/ssebop/test_quality.py b/tests/ssebop/test_quality.py index 6d2e1fd..b4a4dc1 100644 --- a/tests/ssebop/test_quality.py +++ b/tests/ssebop/test_quality.py @@ -4,7 +4,7 @@ Objective: Verify Landsat QA retention, cloud rejection and scene processing with tiny GeoTIFFs. Author: Yi Yu (with assistance from Codex) Created: 2026-10-16 -Last updated: 2026-10-16 +Last updated: 2026-10-22 Inputs: Synthetic Collection 2 flags and temporary raw Landsat rasters. Outputs: Regression assertions for masks, zero-valued QA and NDVI exports. Usage: python -m pytest tests/ssebop/test_quality.py @@ -65,7 +65,7 @@ def run_scene(path, monkeypatch, tmp_path, **kwargs): monkeypatch.setattr(api, "tcold_fano_simple_xr", lambda lst, ndvi, dt, **kw: lst) return api.process_landsat_scene( str(path), "ST_B10", "ndvi", "SR_B4", "SR_B5", False, None, dt, - template.rio.crs, str(tmp_path), str(tmp_path), **kwargs) + template.rio.crs, str(tmp_path), str(tmp_path), tcold_method="legacy", **kwargs) def test_scene_keeps_zero_saturation_and_masks_contaminated_ndvi(tmp_path, monkeypatch): diff --git a/tests/test_protocols.py b/tests/test_protocols.py index 9b737ac..e75bccb 100644 --- a/tests/test_protocols.py +++ b/tests/test_protocols.py @@ -3,7 +3,7 @@ Objective: Verify holdout isolation, donor restrictions, spin-up and evaluation contracts. Author: Yi Yu (with assistance from Codex) Created: 2026-10-17 -Last updated: 2026-10-21 +Last updated: 2026-10-22 Inputs: Small deterministic synthetic scenes, forcing and observation arrays. Outputs: Scientific protocol regression assertions. Usage: python -m pytest tests/test_protocols.py @@ -305,7 +305,7 @@ def test_l7_file_gate_and_receipt_reject_mutation(tmp_path): (tmp_path/'etf').mkdir() (tmp_path/'ndvi').mkdir() process_landsat_scene(records[0]['file'],'auto','ndvi','SR_B3','SR_B4',False,None, - delta,template.rio.crs,str(tmp_path/'etf'),str(tmp_path/'ndvi'),gapfill_product=outputs[records[0]['file']]) + delta,template.rio.crs,str(tmp_path/'etf'),str(tmp_path/'ndvi'),gapfill_product=outputs[records[0]['file']],tcold_method='legacy') with xr.open_dataset(tmp_path/'etf/auxiliary_2014-02-01.nc') as ds: assert ds.lst_K_provenance.values[20,20] == 2 assert ds.original_support.values[20,20] == 0 @@ -334,6 +334,8 @@ def test_calendar_climatology_requires_full_source_and_ignores_holdout(): def test_screened_scene_workflow_offline(tmp_path): + import json + from pathlib import Path import rasterio from rasterio import Affine from rasterio.warp import transform_bounds @@ -354,6 +356,8 @@ def test_screened_scene_workflow_offline(tmp_path): scenes.append(str(path)) dem=tmp_path/'dem.tif' (source_grid*100).rio.to_raster(dem) + landcover=tmp_path/'worldcover.tif' + (source_grid*40).rio.to_raster(landcover) geographic=raster(np.ones((6,6)),Affine(.1,0,8.7,0,-.1,55.3),'EPSG:4326') dates=pd.date_range('2014-01-01','2014-12-31') files={} @@ -364,7 +368,7 @@ def test_screened_scene_workflow_offline(tmp_path): padding=dict(primary_bounds=[497120,6097120,502880,6102880],comparison_bounds=list(source_grid.rio.bounds()), crs='EPSG:32632',maximum_fallback_fraction=1e-8,minimum_agreement=.99,maximum_p99_difference=.01) protocol=dict(development_end='2014-12-31',climatology_start='2014-01-01',climatology_end='2014-12-31', - meteorology_files=files,scene_cache_dir=str(tmp_path/'cache'),padding=padding) + meteorology_files=files,scene_cache_dir=str(tmp_path/'cache'),padding=padding,landcover_file=str(landcover)) output=run_scene_workflow(date_range='2014-06-01 to 2014-06-05',scene_files=scenes,gapfill_products={}, dem=str(dem),output_dir=tmp_path/'seb',extent=extent,protocol=protocol, reconstruction=dict(acquisition_max_gap_days=64,boundary_carry=True),workers=1) @@ -373,6 +377,19 @@ def test_screened_scene_workflow_offline(tmp_path): assert ds.ET.notnull().any() assert {'original_supported','reconstructed','dt_K','cold_boundary_fallback_fraction'} <= set(ds) assert ds.rio.crs == source_grid.rio.crs + assert ds.attrs['tcold_method'] == 'worldcover_fano' + assert ds.worldcover_class.min() == 40 + # Cached retrievals must invalidate when categorical input content changes, + # even if the replacement class has the same reference eligibility. + from pysweb.ssebop.scene_protocol import _retrieve_scene + receipt=tmp_path/'cache'/Path(scenes[0]).stem/'2014-06-01.json' + fingerprint=json.loads(receipt.read_text())['fingerprint'] + original_mtime=receipt.stat().st_mtime_ns + _retrieve_scene(scenes[0]) + assert receipt.stat().st_mtime_ns == original_mtime + (source_grid*30).rio.to_raster(landcover) + _retrieve_scene(scenes[0]) + assert json.loads(receipt.read_text())['fingerprint'] != fingerprint with pytest.raises(ValueError,match='before first supported scene'): run_scene_workflow(date_range='2014-05-31 to 2014-06-05',scene_files=scenes,gapfill_products={}, dem=str(dem),output_dir=tmp_path/'early',extent=extent,