From 950949929d5abb4f1f5f21b525d4f0fcabfb382e Mon Sep 17 00:00:00 2001 From: Manjila Singh Date: Tue, 28 Jul 2026 08:44:51 -0500 Subject: [PATCH 1/4] add fimbox.nextgen which gets nextgen hydrofabric catchments and datastream discharge for a given aoi --- src/fimbox/nextgen/README.md | 135 ++++++++++ src/fimbox/nextgen/__init__.py | 47 ++++ src/fimbox/nextgen/__main__.py | 17 ++ src/fimbox/nextgen/_common.py | 113 ++++++++ src/fimbox/nextgen/data/vpu_bbox.json | 28 ++ src/fimbox/nextgen/datastream.py | 312 ++++++++++++++++++++++ src/fimbox/nextgen/hydrofabric.py | 358 ++++++++++++++++++++++++++ src/fimbox/nextgen/pipeline.py | 321 +++++++++++++++++++++++ 8 files changed, 1331 insertions(+) create mode 100644 src/fimbox/nextgen/README.md create mode 100644 src/fimbox/nextgen/__init__.py create mode 100644 src/fimbox/nextgen/__main__.py create mode 100644 src/fimbox/nextgen/_common.py create mode 100644 src/fimbox/nextgen/data/vpu_bbox.json create mode 100644 src/fimbox/nextgen/datastream.py create mode 100644 src/fimbox/nextgen/hydrofabric.py create mode 100644 src/fimbox/nextgen/pipeline.py diff --git a/src/fimbox/nextgen/README.md b/src/fimbox/nextgen/README.md new file mode 100644 index 0000000..6037ada --- /dev/null +++ b/src/fimbox/nextgen/README.md @@ -0,0 +1,135 @@ +### NextGen-in-a-Box Integration + +This part of fimbox connects your **area of interest** (a map outline of the place you care about) to NextGen streamflow data, so that flood maps can be produced for that area. + +In simple terms, you hand it a shapefile or geopackage of an area, and it figures out which rivers and drainage areas fall inside that shape, then goes and fetches how much water is flowing in each of them. + +**What you give it** + +- An **area of interest** — a boundary file (a `.shp` shapefile or a `.gpkg` GeoPackage) outlining the region you want. For example, a city, a watershed, or a study site. + +**What you get back**, saved in a folder named after your area: + +- A **catchments** map file — the drainage areas that fall inside your boundary. +- A **discharge** table — how much water is flowing in each of those areas (in cubic metres per second). +- The **stream IDs** — the identifiers that link each drainage area to its river reach. + +Everything is downloaded from a [free public s3 bucket](https://ciroh-community-ngen-datastream.s3.amazonaws.com/index.html). + +### How to run it + +The simplest way, from a terminal: + +```bash +python -m fimbox.nextgen my_area.shp --out-dir out +``` + +Replace `my_area.shp` with the path to your boundary file. Results land in `out/my_area/`. + +Or from Python, in one line: + +```python +import fimbox + +result = fimbox.getNextGenAOI("my_area.shp", out_dir="out") + +print(result.catchments_path) # the catchments map file +print(result.discharge_csvs) # the water-flow table(s) +print(result.feature_ids) # the stream IDs +``` + +By default it grabs the **most recent** available forecast. + +**Output folder structure** + +``` +out/my_area/ + feature_id.csv the list of stream IDs + hydrofabric/ + aoi_catchments.gpkg the catchments (drainage areas) map + aoi_flowpaths.gpkg the rivers/streams inside your area + network_crosswalk.csv a lookup table linking IDs together + discharge-inputs/ + nextgen_..._maximum.csv the water-flow numbers for each stream +``` + +### A few useful options + +You can add these to the command above: + +- `--no-discharge` — just find the catchments and IDs, skip downloading the water-flow data (faster). +- `--date 20260721 --cycle 09` — use a specific day/time instead of the latest. +- `--forecast medium_range` — use a longer-range forecast (default is `short_range`). +- `--sortby maximum` — for each stream, keep the **peak** flow over the forecast (the default, best for flood extent). Other choices: `minimum`, `mean`, or `none` to keep every hour. + + +--- + +### Developer reference + +`fimbox.nextgen` maps an AOI to the NOAA/OWP **NextGen v2.2 hydrofabric** and the ngen/t-route **discharge** published on the public [CIROH community NextGen DataStream](https://ciroh-community-ngen-datastream.s3.amazonaws.com) S3 bucket. + +**How the crosswalk works.** Each NextGen divide (`divide_id` = `cat-`) drains to exactly one flowpath (`id` = `wb-`). The integer `` is the `feature_id` used in the ngen/t-route outputs, so catchment → flowpath → discharge is a direct join. The `network` layer also carries `hf_id` (the NWM COMID), available as an optional crosswalk (`nwm_crosswalk=True`, slower because it scans that table over S3). + +```mermaid +flowchart LR + A[AOI shp/gpkg] --> B[Resolve VPU
cached bbox index] + B --> C[Read intersecting
divides from S3 gpkg] + C --> D[cat-* -> wb-* -> feature_id
+ optional NWM hf_id] + D --> E[Locate ngen run
model/forecast/date/cycle/VPU] + E --> F[Read t-route flow
parquet or tar.gz/netCDF] + F --> G[FIM-ready CSVs
feature_id, discharge_cms] + D --> H[aoi_catchments.gpkg
network_crosswalk.csv] +``` + +Newer runs store discharge as `ngen-run/outputs/troute/troute_output_*.parquet`; older runs as an `ngen-run.tar.gz` containing `outputs/troute/*.nc`. Both are handled transparently. + +**Module contents** + +| File | What it contains | +|---|---| +| `hydrofabric.py` | `NextGenHydrofabric` / `AOIHydrofabric`: AOI → VPU resolution, intersecting `divides`, and the `cat`/`wb`/`feature_id`(+`hf_id`) crosswalk. `build_vpu_index()` regenerates the cached VPU bbox index. | +| `datastream.py` | `NextGenDatastream` / `DischargeRun`: locate an ngen run on the bucket (default `cfe_nom` / `short_range` / latest), read t-route `flow` for the feature ids, and write FIM-ready CSVs. | +| `pipeline.py` | `NextGenAOI` and `getNextGenAOI`: AOI → catchments + discharge in one call, written into the standard fimbox AOI layout. Also the module CLI. | +| `_common.py` | Bucket constants, anonymous S3 handle, AOI-layout helpers, the cached VPU bbox index, and `wb`→`feature_id` parsing. | +| `data/vpu_bbox.json` | Per-VPU bounding boxes (EPSG:5070) so AOI→VPU resolution needs no network round-trip. | + +**Full Python API** + +```python +import fimbox + +res = fimbox.getNextGenAOI( + "my_aoi.shp", # AOI shapefile / GeoPackage / GeoJSON + out_dir="out", # AOI folder created at out/ + # aoi_layer=None, # layer name for a multi-layer gpkg + # predicate="intersects", # or "within" (catchments fully inside the AOI) + # nwm_crosswalk=False, # also resolve NWM hf_id/COMID (slower S3 read) + # model="cfe_nom", # ngen model output set + # forecast="short_range", # short_range | medium_range | analysis_assim_extend + # date=None, cycle=None, # YYYYMMDD + cycle hour; default latest available + # sortby="maximum", # horizon aggregation: maximum|minimum|mean|None + # at_time=None, # single timestamp instead of aggregating + # fetch_discharge=True, # False -> resolve catchments/ids only +) + +res.catchments_path # /hydrofabric/aoi_catchments.gpkg +res.feature_ids # NextGen feature_ids == streamflow network ids +res.network_ids # wb-* flowpath ids +res.discharge # long DataFrame [feature_id, time, flow] (flow in cms) +res.discharge_csvs # FIM-ready CSVs in /discharge-inputs/ +res.run # e.g. "cfe_nom/short_range/ngen.20260721/09/VPU_03N" + +# Step-by-step (class-based) +from fimbox.nextgen import NextGenHydrofabric, NextGenDatastream + +hf = NextGenHydrofabric("my_aoi.shp").resolve() # AOIHydrofabric +ds = NextGenDatastream(hf.vpus[0]) # latest cfe_nom short_range run +flow = ds.read_discharge(hf.feature_ids) # [feature_id, time, flow] +csvs = ds.to_fim_inputs("out/my_aoi", hf.feature_ids, sortby="maximum") +``` + +**Notes** + +- The bucket is read anonymously; no AWS credentials are required. +- `feature_id.csv` and `discharge-inputs/` match the `fimbox.streamflow` layout, so NextGen discharge is a drop-in alternative to the NWM/GEOGLOWS sources for FIM generation. diff --git a/src/fimbox/nextgen/__init__.py b/src/fimbox/nextgen/__init__.py new file mode 100644 index 0000000..a2bd820 --- /dev/null +++ b/src/fimbox/nextgen/__init__.py @@ -0,0 +1,47 @@ +""" +Author: Manjila Singh +Date Updated: July 2026 + +NextGen-in-a-Box integration for fimbox. + +Map an AOI to the NOAA/OWP NextGen v2.2 hydrofabric and the +ngen/t-route discharge published on the public CIROH community NextGen +DataStream bucket (https://ciroh-community-ngen-datastream.s3.amazonaws.com), +so the resulting streamflow can drive FIM generation. + +Resolution +---------- +NextGenHydrofabric AOI (shp/gpkg) -> intersecting catchments + network-id crosswalk +NextGenDatastream NextGen network ids -> ngen/t-route discharge (S3) + +Orchestration +------------- +NextGenAOI / getNextGenAOI AOI -> catchments + discharge in one call, written + into the standard fimbox AOI layout. + +Example +------- + import fimbox + + res = fimbox.getNextGenAOI("my_aoi.shp", out_dir="out") + res.catchments_path # /hydrofabric/aoi_catchments.gpkg + res.feature_ids # NextGen feature_ids (== streamflow network ids) + res.discharge_csvs # FIM-ready /discharge-inputs/*.csv +""" + +from __future__ import annotations + +from .datastream import DischargeRun, NextGenDatastream +from .hydrofabric import AOIHydrofabric, NextGenHydrofabric, build_vpu_index +from .pipeline import NextGenAOI, NextGenResult, getNextGenAOI + +__all__ = [ + "NextGenHydrofabric", + "AOIHydrofabric", + "NextGenDatastream", + "DischargeRun", + "NextGenAOI", + "NextGenResult", + "getNextGenAOI", + "build_vpu_index", +] diff --git a/src/fimbox/nextgen/__main__.py b/src/fimbox/nextgen/__main__.py new file mode 100644 index 0000000..f5b0c18 --- /dev/null +++ b/src/fimbox/nextgen/__main__.py @@ -0,0 +1,17 @@ +""" +Author: Manjila Singh +Date Created: July 2026 + +Package entry point so ``python -m fimbox.nextgen`` runs the pipeline CLI. + +Kept separate from ``pipeline.py`` (and not imported by ``__init__``) so that +``python -m fimbox.nextgen`` does not trip the runpy "found in sys.modules" +double-import warning. +""" + +from __future__ import annotations + +from .pipeline import _main + +if __name__ == "__main__": + _main() diff --git a/src/fimbox/nextgen/_common.py b/src/fimbox/nextgen/_common.py new file mode 100644 index 0000000..8deeb79 --- /dev/null +++ b/src/fimbox/nextgen/_common.py @@ -0,0 +1,113 @@ +""" +Shared helpers for the ``fimbox.nextgen`` subpackage: constants for the public +CIROH community NextGen DataStream bucket, an anonymous s3fs handle, AOI-layout +resolution, and the cached per-VPU bounding-box index used to resolve an area +of interest to its hydrofabric VPU(s) without opening every GeoPackage. + +The bucket is read anonymously (``--no-sign-request`` equivalent) so it works +for any caller regardless of whether AWS credentials are configured. +""" + +from __future__ import annotations + +import importlib +import json +import logging +from functools import lru_cache +from pathlib import Path +from typing import Optional, Union + +from ..logging_utils import aoi_root, attach_case_log + +log = logging.getLogger(__name__) + +PathLike = Union[str, Path] + +BUCKET = "ciroh-community-ngen-datastream" +S3_BUCKET_URL = f"https://{BUCKET}.s3.amazonaws.com" + +HF_VERSION = "v2.2" +HF_GEOPACKAGES_PREFIX = f"{BUCKET}/resources/{HF_VERSION}_hydrofabric/geopackages" + +# outputs//_hydrofabric/ngen.///VPU_/ +OUTPUTS_PREFIX = f"{BUCKET}/outputs" + +HF_EPSG = 5070 # CONUS Albers — all spatial ops happen in this CRS +DISCHARGE_COL = "discharge_cms" # column name the FIM Inundator expects + +HYDROFABRIC_DIR = "hydrofabric" +DISCHARGE_INPUTS_DIR = "discharge-inputs" + + +def require(module: str): + """Import an optional dependency, raising a clear install hint if missing.""" + try: + return importlib.import_module(module) + except ImportError as exc: # pragma: no cover - trivial + raise ImportError( + f"'{module}' is required for fimbox.nextgen. " + f"Install it with: pip install {module}" + ) from exc + + +@lru_cache(maxsize=1) +def s3() -> "object": + """Anonymous s3fs filesystem for the public bucket (cached).""" + s3fs = require("s3fs") + return s3fs.S3FileSystem(anon=True) + + +def resolve_aoi(aoi_dir: PathLike) -> Path: + """Return the AOI root for any directory the caller passes (the root itself + or its ``watershed-data`` subfolder).""" + return aoi_root(Path(aoi_dir)) + + +def hydrofabric_dir(aoi_dir: PathLike) -> Path: + """``/hydrofabric`` — subset catchments / flowpaths / crosswalk land + here. Created on demand.""" + d = resolve_aoi(aoi_dir) / HYDROFABRIC_DIR + d.mkdir(parents=True, exist_ok=True) + return d + + +def discharge_inputs_dir(aoi_dir: PathLike) -> Path: + """``/discharge-inputs`` — the FIM-ready discharge CSVs the generator + iterates (shared with the streamflow subpackage).""" + d = resolve_aoi(aoi_dir) / DISCHARGE_INPUTS_DIR + d.mkdir(parents=True, exist_ok=True) + return d + + +def attach_log(aoi_dir: PathLike) -> None: + """Route nextgen log records into the AOI's combined processing.log.""" + attach_case_log(aoi_dir) + + +_DATA_DIR = Path(__file__).parent / "data" +_VPU_INDEX_PATH = _DATA_DIR / "vpu_bbox.json" + + +@lru_cache(maxsize=1) +def vpu_bbox_index() -> dict[str, list[float]]: + """Load the cached ``{VPU_id: [minx, miny, maxx, maxy]}`` index (EPSG:5070). + + Shipped as package data so AOI->VPU resolution needs no network round-trip. + """ + payload = json.loads(_VPU_INDEX_PATH.read_text()) + return {k: list(map(float, v)) for k, v in payload["bboxes"].items()} + + +def vpu_gpkg_uri(vpu: str) -> str: + """``/vsis3`` URI for a VPU's NextGen GeoPackage (read anonymously).""" + return f"/vsis3/{HF_GEOPACKAGES_PREFIX}/{vpu}/nextgen_{vpu}.gpkg" + + +def feature_id_from_wb(wb_id: str) -> Optional[int]: + """``"wb-2855078"`` -> ``2855078``. Returns None for non-``wb`` ids.""" + if not isinstance(wb_id, str) or "-" not in wb_id: + return None + prefix, _, num = wb_id.partition("-") + if prefix != "wb" or not num.isdigit(): + return None + return int(num) diff --git a/src/fimbox/nextgen/data/vpu_bbox.json b/src/fimbox/nextgen/data/vpu_bbox.json new file mode 100644 index 0000000..8a70251 --- /dev/null +++ b/src/fimbox/nextgen/data/vpu_bbox.json @@ -0,0 +1,28 @@ +{ + "_comment": "NextGen v2.2 hydrofabric per-VPU bounding boxes in EPSG:5070 [minx, miny, maxx, maxy]. Derived from the divides layer total_bounds of each nextgen_VPU_.gpkg in s3://ciroh-community-ngen-datastream/resources/v2.2_hydrofabric/geopackages/. Used to resolve an AOI to its VPU(s) without opening every gpkg. Regenerate with fimbox.nextgen.hydrofabric.build_vpu_index().", + "crs": "EPSG:5070", + "hydrofabric_version": "v2.2", + "bboxes": { + "VPU_01": [1824824.9997, 2215544.9958, 2258234.9955, 3086925.0003], + "VPU_02": [1350645.0003, 1687635.0, 1991054.9979, 2720864.9997], + "VPU_03N": [1115894.9997, 1023914.9979, 1833435.0045, 1740465.0], + "VPU_03S": [1061444.9997, 340785.0, 1598235.0003, 1331925.0003], + "VPU_03W": [578895.0003, 788805.0, 1161915.0003, 1402994.9997], + "VPU_04": [213255.0, 1990604.9997, 1744665.0003, 2824454.9997], + "VPU_05": [588854.9997, 1413285.0003, 1488224.9997, 2292855.0003], + "VPU_06": [645014.9997, 1272255.0003, 1294065.0, 1672544.9997], + "VPU_07": [-106125.0003, 1565505.0, 837735.0003, 2754015.0003], + "VPU_08": [150975.0, 673394.9976, 689775.0003, 1661985.0], + "VPU_09": [-646274.9997, 2498744.9997, 455654.9997, 3083384.9997], + "VPU_10L": [-1074915.0, 1556235.0, 505784.9997, 2300715.0], + "VPU_10U": [-1389615.0003, 2117324.9997, 31545.0, 3039464.9997], + "VPU_11": [-907575.0003, 908984.9997, 518355.0, 1866134.9997], + "VPU_12": [-729254.9997, 312795.0, 271244.9997, 1319204.9997], + "VPU_13": [-1206315.0, 209715.0003, -115845.0003, 1759995.0], + "VPU_14": [-1425105.0, 1456994.9997, -811784.9997, 2356695.0], + "VPU_15": [-1772655.0003, 975944.9997, -1081784.9997, 1971974.9997], + "VPU_16": [-2063565.0, 1528364.9997, -1192725.0, 2368545.0003], + "VPU_17": [-2294145.0, 2191574.9997, -1102155.0003, 3506235.0003], + "VPU_18": [-2355824.9997, 1202325.0003, -1696155.0003, 2538315.0] + } +} diff --git a/src/fimbox/nextgen/datastream.py b/src/fimbox/nextgen/datastream.py new file mode 100644 index 0000000..24b55bc --- /dev/null +++ b/src/fimbox/nextgen/datastream.py @@ -0,0 +1,312 @@ +""" +Retrieve NextGen-in-a-Box discharge from the public CIROH community NextGen +DataStream bucket for a set of hydrofabric ``feature_id``s. + +Output layout on the bucket:: + + outputs//v2.2_hydrofabric/ngen.///VPU_/ + ngen-run/outputs/troute/troute_output_.parquet (newer runs) + ngen-run.tar.gz -> ngen-run/outputs/troute/*.nc (older runs) + +t-route output is the routed streamflow that drives FIM. It is keyed by +``feature_id`` (the integer of the ``wb-`` flowpath id) with an hourly +``flow`` series in cubic metres per second. This module locates the run +(defaulting to model ``cfe_nom``, ``short_range``, and the latest available +date/cycle that has the requested VPU), reads the ``flow`` series for the +requested feature ids, and emits FIM-ready CSVs (``feature_id, discharge_cms``) +into ``/discharge-inputs/``. +""" + +from __future__ import annotations + +import io +import logging +import tarfile +from dataclasses import dataclass +from pathlib import Path +from typing import Optional, Union + +import pandas as pd + +from . import _common as C + +log = logging.getLogger(__name__) + +PathLike = Union[str, Path] + +DEFAULT_MODEL = "cfe_nom" +DEFAULT_FORECAST = "short_range" +_HF_SEG = f"{C.HF_VERSION}_hydrofabric" + + +@dataclass +class DischargeRun: + """A located ngen datastream run.""" + + model: str + forecast: str + date: str # YYYYMMDD + cycle: str # e.g. "00" + vpu: str + + @property + def prefix(self) -> str: + return ( + f"{C.OUTPUTS_PREFIX}/{self.model}/{_HF_SEG}/ngen.{self.date}/" + f"{self.forecast}/{self.cycle}/{self.vpu}" + ) + + def __str__(self) -> str: + return f"{self.model}/{self.forecast}/ngen.{self.date}/{self.cycle}/{self.vpu}" + + +class NextGenDatastream: + """Locate and read NextGen datastream discharge for one VPU. + + Parameters + ---------- + vpu : str + Hydrofabric VPU id, e.g. ``"VPU_03N"``. + model : str + ngen model output set. Default ``"cfe_nom"``. + forecast : str + Forecast product: ``"short_range"`` (default), ``"medium_range"``, + ``"analysis_assim_extend"``, ... + date : str, optional + Run date ``YYYYMMDD`` (or ``YYYY-MM-DD``). Default: latest available. + cycle : str, optional + Cycle hour, e.g. ``"06"``. Default: latest cycle that has ``vpu``. + """ + + def __init__( + self, + vpu: str, + *, + model: str = DEFAULT_MODEL, + forecast: str = DEFAULT_FORECAST, + date: Optional[str] = None, + cycle: Optional[str] = None, + ): + self.vpu = vpu + self.model = model + self.forecast = forecast + self._fs = C.s3() + self.run = self._locate(date, cycle) + + def _model_root(self) -> str: + return f"{C.OUTPUTS_PREFIX}/{self.model}/{_HF_SEG}" + + def _ls_names(self, prefix: str) -> list[str]: + try: + return [p.split("/")[-1] for p in self._fs.ls(prefix)] + except FileNotFoundError: + return [] + + def _available_dates(self) -> list[str]: + dates = [ + n.split("ngen.")[-1] + for n in self._ls_names(self._model_root()) + if n.startswith("ngen.") + ] + return sorted(dates, reverse=True) + + def _locate(self, date: Optional[str], cycle: Optional[str]) -> DischargeRun: + date = date.replace("-", "") if date else None + + def cycle_has_vpu(d: str, c: str) -> bool: + base = f"{self._model_root()}/ngen.{d}/{self.forecast}/{c}/{self.vpu}" + return self._fs.exists(base) + + if date and cycle: + if not cycle_has_vpu(date, cycle): + raise FileNotFoundError( + f"No {self.vpu} output for {self.model}/{self.forecast} " + f"ngen.{date} cycle {cycle}." + ) + return DischargeRun(self.model, self.forecast, date, cycle, self.vpu) + + candidate_dates = [date] if date else self._available_dates() + if not candidate_dates: + raise FileNotFoundError( + f"No dated runs under {self._model_root()} (model={self.model})." + ) + for d in candidate_dates: + cycles = self._ls_names(f"{self._model_root()}/ngen.{d}/{self.forecast}") + wanted = [cycle] if cycle else sorted(cycles, reverse=True) + for c in wanted: + if c in cycles and cycle_has_vpu(d, c): + log.info( + "Located discharge run: %s/%s ngen.%s cycle %s %s", + self.model, + self.forecast, + d, + c, + self.vpu, + ) + return DischargeRun(self.model, self.forecast, d, c, self.vpu) + raise FileNotFoundError( + f"No {self.forecast} run with {self.vpu} found for model " + f"{self.model} (searched {len(candidate_dates)} date(s))." + ) + + def _troute_dir(self) -> str: + return f"{self.run.prefix}/ngen-run/outputs/troute" + + def read_discharge(self, feature_ids: Optional[list[int]] = None) -> pd.DataFrame: + """Return a long DataFrame ``[feature_id, time, flow]`` (flow in cms), + filtered to ``feature_ids`` when given. + + Handles both the extracted parquet/netCDF layout (newer runs) and the + ``ngen-run.tar.gz`` archive layout (older runs). + """ + want = set(int(f) for f in feature_ids) if feature_ids is not None else None + troute_dir = self._troute_dir() + + if self._fs.exists(troute_dir): + files = self._fs.ls(troute_dir) + df = self._read_extracted(files) + elif self._fs.exists(f"{self.run.prefix}/ngen-run.tar.gz"): + df = self._read_tarball(f"{self.run.prefix}/ngen-run.tar.gz") + else: + raise FileNotFoundError( + f"No troute output (extracted or tarball) under {self.run.prefix}" + ) + + if want is not None: + df = df[df["feature_id"].isin(want)].copy() + df = df.sort_values(["feature_id", "time"]).reset_index(drop=True) + log.info( + "Discharge: %d reaches x %d timesteps from %s", + df["feature_id"].nunique(), + df["time"].nunique(), + self.run, + ) + return df + + def _read_extracted(self, files: list[str]) -> pd.DataFrame: + parquet = [f for f in files if f.endswith(".parquet")] + netcdf = [f for f in files if f.endswith(".nc")] + if parquet: + frames = [pd.read_parquet(self._fs.open(f)) for f in parquet] + return self._normalize(pd.concat(frames, ignore_index=True)) + if netcdf: + frames = [self._read_netcdf(self._fs.open(f)) for f in netcdf] + return self._normalize(pd.concat(frames, ignore_index=True)) + raise FileNotFoundError("troute dir has no .parquet or .nc output") + + def _read_tarball(self, key: str) -> pd.DataFrame: + log.info("Downloading archive %s ...", key.split("/")[-1]) + raw = self._fs.cat_file(key) + frames = [] + with tarfile.open(fileobj=io.BytesIO(raw), mode="r:gz") as tf: + members = [ + m + for m in tf.getmembers() + if "outputs/troute" in m.name and m.name.endswith((".nc", ".parquet")) + ] + for m in members: + fh = tf.extractfile(m) + buf = io.BytesIO(fh.read()) + if m.name.endswith(".parquet"): + frames.append(pd.read_parquet(buf)) + else: + frames.append(self._read_netcdf(buf)) + if not frames: + raise FileNotFoundError(f"{key} has no troute output members") + return self._normalize(pd.concat(frames, ignore_index=True)) + + @staticmethod + def _read_netcdf(fileobj) -> pd.DataFrame: + xr = C.require("xarray") + # netCDF4 can't read a file-like object; materialise to a temp file + import tempfile + + data = fileobj.read() if hasattr(fileobj, "read") else fileobj + with tempfile.NamedTemporaryFile(suffix=".nc") as tmp: + tmp.write(data) + tmp.flush() + with xr.open_dataset(tmp.name) as ds: + keep = [v for v in ("flow", "velocity", "depth") if v in ds] + df = ds[keep].to_dataframe().reset_index() + return df + + @staticmethod + def _normalize(df: pd.DataFrame) -> pd.DataFrame: + """Coerce whatever troute layout we read into ``[feature_id, time, flow]``.""" + cols = {c.lower(): c for c in df.columns} + fid = cols.get("feature_id") or cols.get("featureid") + flow = cols.get("flow") or cols.get("streamflow") or cols.get("q") + time = cols.get("time") or cols.get("current_time") + if fid is None or flow is None: + raise ValueError(f"Unrecognized troute columns: {list(df.columns)}") + # keep flowpaths (wb), drop nexus rows when a 'type' column is present + if "type" in cols: + df = df[df[cols["type"]].astype(str).str.lower() == "wb"] + out = df.rename(columns={fid: "feature_id", flow: "flow"}) + if time is not None: + out = out.rename(columns={time: "time"}) + else: + out["time"] = pd.NaT + out["feature_id"] = out["feature_id"].astype("int64") + out["time"] = pd.to_datetime(out["time"]) + return out[["feature_id", "time", "flow"]] + + def to_fim_inputs( + self, + aoi_dir: PathLike, + feature_ids: list[int], + *, + sortby: Optional[str] = "maximum", + at_time: Optional[str] = None, + ) -> list[Path]: + """Write FIM-ready discharge CSV(s) to ``/discharge-inputs/``. + + ``sortby`` in {``"maximum"``, ``"minimum"``, ``"mean"``} collapses the + forecast horizon into one CSV of that statistic per reach (default + ``"maximum"`` — the peak, most relevant for flood extent). ``at_time`` + selects a single timestamp instead. ``sortby=None`` and no ``at_time`` + writes one CSV per timestep. + """ + C.attach_log(aoi_dir) + df = self.read_discharge(feature_ids) + if df.empty: + log.warning("No discharge matched the requested feature_ids.") + return [] + out_dir = C.discharge_inputs_dir(aoi_dir) + tag = ( + f"nextgen_{self.model}_{self.run.forecast}_{self.run.date}{self.run.cycle}" + ) + written: list[Path] = [] + + if at_time is not None: + ts = pd.to_datetime(at_time) + sub = df[df["time"] == ts] + if sub.empty: + raise ValueError( + f"No timestep at {at_time}; available: " + f"{df['time'].min()}..{df['time'].max()}" + ) + written.append( + self._write_csv(sub, out_dir / f"{tag}_{C.DISCHARGE_COL}.csv") + ) + elif sortby: + agg = {"maximum": "max", "minimum": "min", "mean": "mean"} + if sortby not in agg: + raise ValueError("sortby must be maximum, minimum, mean, or None") + g = df.groupby("feature_id")["flow"].agg(agg[sortby]).reset_index() + written.append(self._write_csv(g, out_dir / f"{tag}_{sortby}.csv")) + else: + for ts, sub in df.groupby("time"): + stamp = pd.to_datetime(ts).strftime("%Y%m%dT%H%M") + written.append(self._write_csv(sub, out_dir / f"{tag}_{stamp}.csv")) + return written + + @staticmethod + def _write_csv(df: pd.DataFrame, path: Path) -> Path: + out = df.rename(columns={"flow": C.DISCHARGE_COL})[ + ["feature_id", C.DISCHARGE_COL] + ].copy() + out["feature_id"] = out["feature_id"].astype("int64") + out.to_csv(path, index=False) + log.info("FIM-ready discharge (%d reaches) --> %s", len(out), path.name) + return path diff --git a/src/fimbox/nextgen/hydrofabric.py b/src/fimbox/nextgen/hydrofabric.py new file mode 100644 index 0000000..30cac10 --- /dev/null +++ b/src/fimbox/nextgen/hydrofabric.py @@ -0,0 +1,358 @@ +""" +Resolve an area of interest (AOI) to the NextGen v2.2 hydrofabric. + +Given an AOI (shapefile / GeoPackage / any vector fiona can read), this module +finds which hydrofabric VPU(s) the AOI falls in, reads the catchment +(``divides``) polygons that intersect the AOI directly from the public +per-VPU GeoPackage on the CIROH community NextGen DataStream bucket, and +crosswalks them to the NextGen network ids used by the ngen/t-route outputs. + +Each NextGen divide (``divide_id`` = ``cat-``) drains to exactly one +flowpath (``id`` = ``wb-``); the integer ```` is the ``feature_id`` used +in the t-route discharge outputs. The ``network`` layer additionally carries +``hf_id`` — the NOAA reference-fabric / NWM COMID — kept as a crosswalk to the +NWM feature ids the rest of fimbox uses. + +Outputs (written under ``/hydrofabric/`` when ``save=True``): + aoi_catchments.gpkg -- intersecting divides (catchment polygons) + aoi_flowpaths.gpkg -- the flowpaths draining those catchments + network_crosswalk.csv -- divide_id, wb_id, feature_id, hf_id (NWM), vpu + feature_id.csv -- unique feature_ids (drop-in for the streamflow step) +""" + +from __future__ import annotations + +import logging +import os +from contextlib import contextmanager +from dataclasses import dataclass, field +from pathlib import Path +from typing import Optional, Union + +import geopandas as gpd +import pandas as pd +from shapely.geometry import box + +from . import _common as C + +log = logging.getLogger(__name__) + +PathLike = Union[str, Path] + +# GDAL/pyogrio env needed to read the public gpkg anonymously over /vsis3. +_VSIS3_ENV = { + "AWS_NO_SIGN_REQUEST": "YES", + "AWS_S3_ENDPOINT": "s3.amazonaws.com", + "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": ".gpkg", + "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR", + "VSI_CACHE": "TRUE", + "GDAL_HTTP_MAX_RETRY": "3", + "GDAL_HTTP_RETRY_DELAY": "1", +} + + +@contextmanager +def _vsis3_env(): + """Temporarily set the GDAL env vars required for anonymous /vsis3 reads.""" + old = {k: os.environ.get(k) for k in _VSIS3_ENV} + os.environ.update(_VSIS3_ENV) + try: + yield + finally: + for k, v in old.items(): + if v is None: + os.environ.pop(k, None) + else: + os.environ[k] = v + + +@dataclass +class AOIHydrofabric: + """Result of an AOI -> hydrofabric resolution. + + Attributes + ---------- + catchments : GeoDataFrame + Intersecting ``divides`` (catchment polygons), EPSG:5070. + flowpaths : GeoDataFrame + The flowpaths (``wb-*``) draining those catchments, EPSG:5070. + crosswalk : DataFrame + One row per catchment: divide_id, wb_id, feature_id, hf_id, vpu. + vpus : list[str] + VPU(s) the AOI intersects (usually one). + """ + + catchments: gpd.GeoDataFrame + flowpaths: gpd.GeoDataFrame + crosswalk: pd.DataFrame + vpus: list[str] = field(default_factory=list) + + @property + def feature_ids(self) -> list[int]: + """Unique t-route ``feature_id``s (integers of the ``wb-*`` ids).""" + return sorted( + {int(f) for f in self.crosswalk["feature_id"].dropna().astype("int64")} + ) + + @property + def catchment_ids(self) -> list[str]: + """``cat-*`` divide ids.""" + return self.crosswalk["divide_id"].dropna().astype(str).tolist() + + @property + def network_ids(self) -> list[str]: + """``wb-*`` NextGen flowpath (network) ids.""" + return self.crosswalk["wb_id"].dropna().astype(str).tolist() + + @property + def nwm_feature_ids(self) -> list[int]: + """NWM/NHD reference COMIDs (``hf_id``) crosswalked to the catchments.""" + return sorted( + {int(f) for f in self.crosswalk["hf_id"].dropna().astype("int64")} + ) + + +class NextGenHydrofabric: + """AOI -> NextGen v2.2 catchments + network-id crosswalk. + + Parameters + ---------- + aoi : str or Path or GeoDataFrame + Area of interest: a path to a shapefile / GeoPackage / GeoJSON, or an + already-loaded GeoDataFrame. + aoi_layer : str, optional + Layer name when ``aoi`` is a multi-layer GeoPackage. + predicate : str, optional + Spatial test used to select catchments. ``"intersects"`` (default) + keeps every catchment that touches the AOI; ``"within"`` keeps only + catchments fully inside the AOI. + """ + + def __init__( + self, + aoi: Union[PathLike, gpd.GeoDataFrame], + *, + aoi_layer: Optional[str] = None, + predicate: str = "intersects", + nwm_crosswalk: bool = False, + ): + if predicate not in ("intersects", "within"): + raise ValueError("predicate must be 'intersects' or 'within'") + self.predicate = predicate + # opt-in: reading the non-spatial `network` table scans the whole layer + self.nwm_crosswalk = nwm_crosswalk + self.aoi = self._load_aoi(aoi, aoi_layer).to_crs(epsg=C.HF_EPSG) + self._aoi_geom = self.aoi.geometry.union_all() + + @staticmethod + def _load_aoi( + aoi: Union[PathLike, gpd.GeoDataFrame], layer: Optional[str] + ) -> gpd.GeoDataFrame: + if isinstance(aoi, gpd.GeoDataFrame): + gdf = aoi.copy() + else: + gdf = gpd.read_file(aoi, layer=layer) if layer else gpd.read_file(aoi) + if gdf.empty: + raise ValueError("AOI is empty.") + if gdf.crs is None: + raise ValueError("AOI has no CRS; set one before resolving.") + return gdf + + def candidate_vpus(self) -> list[str]: + """VPU(s) whose bounding box overlaps the AOI (bbox pre-filter).""" + aoi_box = box(*self.aoi.total_bounds) + hits = [ + vpu + for vpu, bbox in C.vpu_bbox_index().items() + if box(*bbox).intersects(aoi_box) + ] + if not hits: + raise ValueError( + "AOI does not overlap any CONUS hydrofabric VPU bounding box. " + "Check the AOI CRS / extent (must be within CONUS)." + ) + return hits + + def _read_divides(self, vpu: str) -> gpd.GeoDataFrame: + """Read ``divides`` in the AOI bbox from a VPU gpkg over /vsis3, then + apply the precise spatial predicate against the AOI geometry.""" + uri = C.vpu_gpkg_uri(vpu) + bounds = tuple(self.aoi.total_bounds) + with _vsis3_env(): + divides = gpd.read_file(uri, layer="divides", bbox=bounds) + if divides.empty: + return divides + if divides.crs is None: + divides.set_crs(epsg=C.HF_EPSG, inplace=True) + else: + divides = divides.to_crs(epsg=C.HF_EPSG) + if self.predicate == "within": + mask = divides.within(self._aoi_geom) + else: + mask = divides.intersects(self._aoi_geom) + hit = divides[mask].copy() + hit["vpuid"] = hit.get("vpuid", vpu) + hit["vpu"] = vpu + return hit + + def _read_flowpaths(self, vpu: str, wb_ids: list[str]) -> gpd.GeoDataFrame: + """Read the flowpaths draining the selected catchments from a VPU gpkg. + + Uses the AOI bbox spatial filter (fast, index-backed) then keeps only + the ``wb`` ids we selected — an attribute ``WHERE id IN (...)`` over + /vsis3 would scan the whole layer instead.""" + if not wb_ids: + return gpd.GeoDataFrame(geometry=[], crs=f"EPSG:{C.HF_EPSG}") + uri = C.vpu_gpkg_uri(vpu) + bounds = tuple(self.aoi.total_bounds) + with _vsis3_env(): + fp = gpd.read_file(uri, layer="flowpaths", bbox=bounds) + if fp.empty: + return fp + if fp.crs is not None: + fp = fp.to_crs(epsg=C.HF_EPSG) + fp = fp[fp["id"].astype(str).isin(set(wb_ids))].copy() + return fp + + def _read_network_crosswalk(self, vpu: str, wb_ids: list[str]) -> pd.DataFrame: + """Read the NWM ``hf_id`` crosswalk for the selected ``wb`` ids from the + ``network`` layer (non-spatial). One hf_id per flowpath is kept.""" + if not wb_ids: + return pd.DataFrame(columns=["wb_id", "hf_id"]) + uri = C.vpu_gpkg_uri(vpu) + where = "id IN ({})".format(",".join(f"'{i}'" for i in wb_ids)) + try: + import pyogrio + + with _vsis3_env(): + net = pyogrio.read_dataframe( + uri, + layer="network", + where=where, + read_geometry=False, + columns=["id", "hf_id"], + ) + except Exception as exc: + log.warning("network crosswalk read failed for %s: %s", vpu, exc) + return pd.DataFrame(columns=["wb_id", "hf_id"]) + net = net.rename(columns={"id": "wb_id"}) + net = net.dropna(subset=["wb_id"]).drop_duplicates(subset=["wb_id"]) + return net[["wb_id", "hf_id"]] + + def resolve(self) -> AOIHydrofabric: + """Run the full AOI -> catchments + crosswalk resolution.""" + cats, fps, xwalks, vpus = [], [], [], [] + for vpu in self.candidate_vpus(): + hit = self._read_divides(vpu) + if hit.empty: + log.info("VPU %s: no catchments intersect the AOI", vpu) + continue + vpus.append(vpu) + wb_ids = hit["id"].dropna().astype(str).tolist() + log.info("VPU %s: %d catchments intersect the AOI", vpu, len(hit)) + + xwalk = pd.DataFrame( + { + "divide_id": hit["divide_id"].astype(str).values, + "wb_id": hit["id"].astype(str).values, + "vpu": vpu, + } + ) + xwalk["feature_id"] = xwalk["wb_id"].map(C.feature_id_from_wb) + if self.nwm_crosswalk: + nwm = self._read_network_crosswalk(vpu, wb_ids) + xwalk = xwalk.merge(nwm, on="wb_id", how="left") + else: + xwalk["hf_id"] = pd.NA + + cats.append(hit) + fps.append(self._read_flowpaths(vpu, wb_ids)) + xwalks.append(xwalk) + + if not cats: + raise ValueError( + "No hydrofabric catchments intersect the AOI in any candidate VPU." + ) + + catchments = gpd.GeoDataFrame( + pd.concat(cats, ignore_index=True), crs=f"EPSG:{C.HF_EPSG}" + ) + flowpaths = ( + gpd.GeoDataFrame( + pd.concat([f for f in fps if not f.empty], ignore_index=True), + crs=f"EPSG:{C.HF_EPSG}", + ) + if any(not f.empty for f in fps) + else gpd.GeoDataFrame(geometry=[], crs=f"EPSG:{C.HF_EPSG}") + ) + crosswalk = pd.concat(xwalks, ignore_index=True).drop_duplicates( + subset=["divide_id"] + ) + log.info( + "AOI hydrofabric: %d catchments across VPU(s) %s (%d unique feature_ids)", + len(catchments), + ", ".join(vpus), + crosswalk["feature_id"].nunique(), + ) + return AOIHydrofabric(catchments, flowpaths, crosswalk, vpus) + + def resolve_and_save(self, aoi_dir: PathLike) -> AOIHydrofabric: + """Resolve and write the catchment/flowpath/crosswalk outputs under + ``/hydrofabric/`` plus a ``feature_id.csv`` at the AOI root.""" + C.attach_log(aoi_dir) + result = self.resolve() + out = C.hydrofabric_dir(aoi_dir) + + cat_path = out / "aoi_catchments.gpkg" + result.catchments.to_file(cat_path, driver="GPKG") + log.info("catchments (%d) --> %s", len(result.catchments), cat_path.name) + + if not result.flowpaths.empty: + fp_path = out / "aoi_flowpaths.gpkg" + result.flowpaths.to_file(fp_path, driver="GPKG") + log.info("flowpaths (%d) --> %s", len(result.flowpaths), fp_path.name) + + xwalk_path = out / "network_crosswalk.csv" + result.crosswalk.to_csv(xwalk_path, index=False) + log.info("network crosswalk --> %s", xwalk_path.name) + + fid_path = C.resolve_aoi(aoi_dir) / "feature_id.csv" + pd.DataFrame({"feature_id": result.feature_ids}).to_csv(fid_path, index=False) + log.info("feature ids (%d) --> %s", len(result.feature_ids), fid_path.name) + return result + + +def build_vpu_index(out_path: Optional[PathLike] = None) -> dict[str, list[float]]: + """Regenerate the cached per-VPU bounding-box index from the public bucket. + + Reads each ``nextgen_VPU_.gpkg`` ``divides`` layer's ``total_bounds`` + (metadata only) over /vsis3 and writes them to ``data/vpu_bbox.json``. Run + this only if the hydrofabric VPU set or extents change upstream. + """ + import json + + import pyogrio + + fs = C.s3() + vpus = sorted( + p.split("/")[-1] + for p in fs.ls(C.HF_GEOPACKAGES_PREFIX) + if "VPU_" in p.split("/")[-1] + ) + bboxes: dict[str, list[float]] = {} + with _vsis3_env(): + for vpu in vpus: + info = pyogrio.read_info(C.vpu_gpkg_uri(vpu), layer="divides") + bboxes[vpu] = [float(x) for x in info["total_bounds"]] + log.info("VPU %s bounds %s", vpu, bboxes[vpu]) + payload = { + "_comment": "EPSG:5070 [minx, miny, maxx, maxy] per VPU; regenerate via " + "fimbox.nextgen.hydrofabric.build_vpu_index().", + "crs": "EPSG:5070", + "hydrofabric_version": C.HF_VERSION, + "bboxes": bboxes, + } + dest = Path(out_path) if out_path else C._VPU_INDEX_PATH + dest.write_text(json.dumps(payload, indent=2)) + C.vpu_bbox_index.cache_clear() + return bboxes diff --git a/src/fimbox/nextgen/pipeline.py b/src/fimbox/nextgen/pipeline.py new file mode 100644 index 0000000..45c78a4 --- /dev/null +++ b/src/fimbox/nextgen/pipeline.py @@ -0,0 +1,321 @@ +""" +NextGen-in-a-Box -> fimbox bridge. + +One entry point takes an area of interest (shapefile / GeoPackage) and returns +everything downstream FIM generation needs from the NextGen ecosystem: + + 1. the hydrofabric catchments (``divides``) that intersect the AOI, + 2. the NextGen network ids (``wb-*`` flowpaths / integer ``feature_id``s) and + their NWM ``hf_id`` crosswalk, and + 3. the ngen/t-route discharge for those reaches, pulled from the public + CIROH community NextGen DataStream S3 bucket and written as FIM-ready + ``feature_id, discharge_cms`` CSVs. + +Everything lands in the standard fimbox AOI layout:: + + / + feature_id.csv -- unique NextGen feature_ids + hydrofabric/ + aoi_catchments.gpkg -- intersecting catchment polygons + aoi_flowpaths.gpkg -- their flowpaths + network_crosswalk.csv -- divide_id, wb_id, feature_id, hf_id, vpu + discharge-inputs/ + nextgen____.csv +""" + +from __future__ import annotations + +import logging +from dataclasses import dataclass +from pathlib import Path +from typing import Optional, Union + +import geopandas as gpd +import pandas as pd + +from . import _common as C +from .datastream import DEFAULT_FORECAST, DEFAULT_MODEL, NextGenDatastream +from .hydrofabric import AOIHydrofabric, NextGenHydrofabric + +log = logging.getLogger(__name__) + +PathLike = Union[str, Path] + + +@dataclass +class NextGenResult: + """Everything the AOI resolution produced. + + Attributes + ---------- + aoi_dir : Path + AOI root the outputs were written under. + hydrofabric : AOIHydrofabric + Catchments, flowpaths, and the network crosswalk. + catchments_path : Path + Saved catchment GeoPackage (the hydrofabric "shapefile"). + crosswalk_path : Path + Saved ``network_crosswalk.csv``. + discharge : DataFrame + Long ``[feature_id, time, flow]`` discharge series (cms). + discharge_csvs : list[Path] + FIM-ready ``feature_id, discharge_cms`` CSVs. + run : str + The datastream run used (``model/forecast/ngen.date/cycle/VPU``). + """ + + aoi_dir: Path + hydrofabric: AOIHydrofabric + catchments_path: Path + crosswalk_path: Path + discharge: pd.DataFrame + discharge_csvs: list[Path] + run: str + + @property + def feature_ids(self) -> list[int]: + return self.hydrofabric.feature_ids + + @property + def network_ids(self) -> list[str]: + return self.hydrofabric.network_ids + + @property + def catchments(self) -> gpd.GeoDataFrame: + return self.hydrofabric.catchments + + +def getNextGenAOI( + aoi: Union[PathLike, gpd.GeoDataFrame], + out_dir: Optional[PathLike] = None, + *, + aoi_layer: Optional[str] = None, + predicate: str = "intersects", + nwm_crosswalk: bool = False, + model: str = DEFAULT_MODEL, + forecast: str = DEFAULT_FORECAST, + date: Optional[str] = None, + cycle: Optional[str] = None, + sortby: Optional[str] = "maximum", + at_time: Optional[str] = None, + fetch_discharge: bool = True, +) -> NextGenResult: + """Resolve an AOI to NextGen catchments + discharge in one call. + + Parameters + ---------- + aoi : path or GeoDataFrame + Area of interest (shapefile / GeoPackage / GeoJSON, or a GeoDataFrame). + out_dir : path, optional + Root output directory; an AOI folder ``/`` is created + underneath. Defaults to ``./out``. + aoi_layer : str, optional + Layer name when the AOI is a multi-layer GeoPackage. + predicate : str + ``"intersects"`` (default) or ``"within"`` for catchment selection. + nwm_crosswalk : bool + Also resolve each reach's NWM/NHD ``hf_id`` (COMID) from the hydrofabric + ``network`` table. Off by default because that read scans the whole + table over S3. + model, forecast, date, cycle + Discharge selection on the datastream bucket. Defaults: ``cfe_nom`` / + ``short_range`` / latest date / latest cycle that has the AOI's VPU. + sortby : str or None + Horizon aggregation for the FIM-ready CSV: ``"maximum"`` (default), + ``"minimum"``, ``"mean"``, or ``None`` for one CSV per timestep. + at_time : str, optional + Single timestamp to slice instead of aggregating. + fetch_discharge : bool + Set False to resolve catchments/ids only (no S3 discharge download). + + Returns + ------- + NextGenResult + """ + return NextGenAOI( + aoi, + out_dir, + aoi_layer=aoi_layer, + predicate=predicate, + nwm_crosswalk=nwm_crosswalk, + ).run( + model=model, + forecast=forecast, + date=date, + cycle=cycle, + sortby=sortby, + at_time=at_time, + fetch_discharge=fetch_discharge, + ) + + +class NextGenAOI: + """AOI -> NextGen hydrofabric catchments + datastream discharge.""" + + def __init__( + self, + aoi: Union[PathLike, gpd.GeoDataFrame], + out_dir: Optional[PathLike] = None, + *, + aoi_layer: Optional[str] = None, + predicate: str = "intersects", + nwm_crosswalk: bool = False, + ): + self.aoi = aoi + self.aoi_layer = aoi_layer + self.predicate = predicate + self.nwm_crosswalk = nwm_crosswalk + + stem = Path(aoi).stem if not isinstance(aoi, gpd.GeoDataFrame) else "aoi" + root = Path(out_dir) if out_dir else (Path.cwd() / "out") + self.aoi_dir = (root / stem).resolve() + self.aoi_dir.mkdir(parents=True, exist_ok=True) + C.attach_log(self.aoi_dir) + + def run( + self, + *, + model: str = DEFAULT_MODEL, + forecast: str = DEFAULT_FORECAST, + date: Optional[str] = None, + cycle: Optional[str] = None, + sortby: Optional[str] = "maximum", + at_time: Optional[str] = None, + fetch_discharge: bool = True, + ) -> NextGenResult: + log.info("=== NextGen AOI: %s ===", self.aoi_dir.name) + + log.info("--- Hydrofabric (AOI -> catchments) ---") + hf = NextGenHydrofabric( + self.aoi, + aoi_layer=self.aoi_layer, + predicate=self.predicate, + nwm_crosswalk=self.nwm_crosswalk, + ) + result = hf.resolve_and_save(self.aoi_dir) + cat_path = C.hydrofabric_dir(self.aoi_dir) / "aoi_catchments.gpkg" + xwalk_path = C.hydrofabric_dir(self.aoi_dir) / "network_crosswalk.csv" + + discharge = pd.DataFrame(columns=["feature_id", "time", "flow"]) + csvs: list[Path] = [] + run_label = "(discharge not fetched)" + + if fetch_discharge: + if len(result.vpus) > 1: + log.warning( + "AOI spans %d VPUs (%s); fetching discharge from each.", + len(result.vpus), + ", ".join(result.vpus), + ) + parts, labels = [], [] + for vpu in result.vpus: + log.info("--- Discharge (%s) ---", vpu) + fids = [ + int(f) + for f in result.crosswalk.loc[ + result.crosswalk["vpu"] == vpu, "feature_id" + ].dropna() + ] + ds = NextGenDatastream( + vpu, model=model, forecast=forecast, date=date, cycle=cycle + ) + parts.append(ds.read_discharge(fids)) + csvs += ds.to_fim_inputs( + self.aoi_dir, fids, sortby=sortby, at_time=at_time + ) + labels.append(str(ds.run)) + if parts: + discharge = pd.concat(parts, ignore_index=True) + run_label = "; ".join(labels) + + log.info("=== DONE: %s ===", self.aoi_dir.name) + return NextGenResult( + aoi_dir=self.aoi_dir, + hydrofabric=result, + catchments_path=cat_path, + crosswalk_path=xwalk_path, + discharge=discharge, + discharge_csvs=csvs, + run=run_label, + ) + + +def _main() -> None: + import argparse + + from ..logging_utils import configure_cli_logging + + configure_cli_logging() + + p = argparse.ArgumentParser( + description="Resolve an AOI to NextGen hydrofabric catchments and pull " + "ngen/t-route discharge from the CIROH community datastream bucket." + ) + p.add_argument("aoi", help="AOI shapefile / GeoPackage / GeoJSON path") + p.add_argument("--aoi-layer", default=None, help="AOI layer (multi-layer gpkg)") + p.add_argument("--out-dir", default="out", help="Root output dir (default ./out)") + p.add_argument( + "--predicate", + default="intersects", + choices=["intersects", "within"], + help="Catchment selection test (default intersects)", + ) + p.add_argument( + "--nwm-crosswalk", + action="store_true", + help="Also resolve NWM/NHD hf_id (COMID) crosswalk (slower S3 read)", + ) + p.add_argument( + "--model", default=DEFAULT_MODEL, help="ngen model (default cfe_nom)" + ) + p.add_argument( + "--forecast", + default=DEFAULT_FORECAST, + help="Forecast product (short_range, medium_range, analysis_assim_extend)", + ) + p.add_argument("--date", default=None, help="Run date YYYYMMDD (default latest)") + p.add_argument("--cycle", default=None, help="Cycle hour, e.g. 06 (default latest)") + p.add_argument( + "--sortby", + default="maximum", + help="Horizon aggregation: maximum|minimum|mean|none (default maximum)", + ) + p.add_argument("--at-time", default=None, help="Single timestamp to slice") + p.add_argument( + "--no-discharge", + action="store_true", + help="Resolve catchments/ids only; skip discharge download", + ) + args = p.parse_args() + + sortby = None if str(args.sortby).lower() == "none" else args.sortby + res = getNextGenAOI( + args.aoi, + args.out_dir, + aoi_layer=args.aoi_layer, + predicate=args.predicate, + nwm_crosswalk=args.nwm_crosswalk, + model=args.model, + forecast=args.forecast, + date=args.date, + cycle=args.cycle, + sortby=sortby, + at_time=args.at_time, + fetch_discharge=not args.no_discharge, + ) + log.info( + "AOI %s: %d catchments, %d feature_ids, VPU(s) %s", + res.aoi_dir.name, + len(res.catchments), + len(res.feature_ids), + ", ".join(res.hydrofabric.vpus), + ) + log.info("Discharge run(s): %s", res.run) + log.info("Catchments: %s", res.catchments_path) + log.info("Crosswalk: %s", res.crosswalk_path) + for c in res.discharge_csvs: + log.info("Discharge: %s", c) + + +if __name__ == "__main__": + _main() From 12751e44468694762d7caaa649c4cad23f634934 Mon Sep 17 00:00:00 2001 From: Manjila Singh Date: Tue, 28 Jul 2026 08:45:56 -0500 Subject: [PATCH 2/4] expose fimbox.nextgen at top level and ship vpu_bbox.json as package data --- pyproject.toml | 5 +++++ src/fimbox/__init__.py | 26 ++++++++++++++++++++++++++ 2 files changed, 31 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 91362bc..35ac40d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -81,6 +81,11 @@ include-package-data = true where = ["src"] exclude = ["tests*"] +[tool.setuptools.package-data] +# Ship the NextGen per-VPU bounding-box index so AOI->VPU resolution needs no +# network round-trip. (Other runtime data is fetched on demand via pooch/S3.) +"fimbox.nextgen" = ["data/*.json"] + [tool.ruff] line-length = 88 target-version = "py310" diff --git a/src/fimbox/__init__.py b/src/fimbox/__init__.py index 1d011ba..724700e 100644 --- a/src/fimbox/__init__.py +++ b/src/fimbox/__init__.py @@ -284,3 +284,29 @@ # On-the-fly access to the bundled reference/calibration lookup tables, fetched # anonymously from the public SDML S3 bucket via pooch (imported at module top). __all__ += ["fetch_data"] + +# NextGen-in-a-Box bridge: resolve an AOI to the NextGen v2.2 hydrofabric +# catchments and pull ngen/t-route discharge from the public CIROH community +# datastream bucket. Heavy/optional deps (s3fs, xarray, geopandas) are imported +# lazily inside the functions, so guard the import to keep a partial install +# from breaking `import fimbox`. +try: + from .nextgen import ( + AOIHydrofabric, + NextGenAOI, + NextGenDatastream, + NextGenHydrofabric, + NextGenResult, + getNextGenAOI, + ) + + __all__ += [ + "getNextGenAOI", + "NextGenAOI", + "NextGenResult", + "NextGenHydrofabric", + "AOIHydrofabric", + "NextGenDatastream", + ] +except Exception: + pass From 30376c797f9ae2394660ae0d1ce2ab48a7810371 Mon Sep 17 00:00:00 2001 From: Manjila Singh Date: Tue, 28 Jul 2026 08:46:14 -0500 Subject: [PATCH 3/4] deleted test files --- tests/README.md | 36 -- tests/conftest.py | 5 - tests/test_branchprocessing.py | 1008 ------------------------------ tests/test_calibrate_pipeline.py | 289 --------- tests/test_downloaddata.py | 203 ------ tests/test_fimevaluation.py | 123 ---- tests/test_fimgeneration.py | 102 --- tests/test_generate_dem_diff.py | 71 --- tests/test_getallinputdata.py | 94 --- tests/test_nwmstreamflow.py | 87 --- tests/test_preprocessDEM.py | 28 - tests/test_preprocessing_hucs.py | 19 - 12 files changed, 2065 deletions(-) delete mode 100644 tests/README.md delete mode 100644 tests/conftest.py delete mode 100644 tests/test_branchprocessing.py delete mode 100644 tests/test_calibrate_pipeline.py delete mode 100644 tests/test_downloaddata.py delete mode 100644 tests/test_fimevaluation.py delete mode 100644 tests/test_fimgeneration.py delete mode 100644 tests/test_generate_dem_diff.py delete mode 100644 tests/test_getallinputdata.py delete mode 100644 tests/test_nwmstreamflow.py delete mode 100644 tests/test_preprocessDEM.py delete mode 100644 tests/test_preprocessing_hucs.py diff --git a/tests/README.md b/tests/README.md deleted file mode 100644 index e9c5433..0000000 --- a/tests/README.md +++ /dev/null @@ -1,36 +0,0 @@ -### Tests -
- -The test suite doubles as the usage reference for `fimbox`: every stage of the workflow has a test file whose top-level constants (AOI paths, input files, worker counts) are meant to be edited to point at your own data. Tests skip cleanly when the referenced AOI or input file is absent. Most files also keep commented-out variants showing alternative call patterns and optional parameters. - -**Folder contents** - -| File | What it covers | -|---|---| -| `conftest.py` | Suite-wide logging setup (`configure_cli_logging`). | -| `test_preprocessing_hucs.py` | HUC validation with `HUCChecker` (single HUC, lists, .lst/.csv files, strict mode). | -| `test_downloaddata.py` | Individual dataset downloaders: DEM, NHDPlus/NWM hydrography, FEMA NFHL, NLD levees, OSM roads/bridges, USGS gages. | -| `test_getallinputdata.py` | The combined `getAllInputData` pipeline from a boundary shapefile, including bring-your-own flowlines/catchments/DEM. | -| `test_preprocessDEM.py` | `DEMProcessor` fetching and conditioning (resolutions, local DEM, CRS handling). | -| `test_generate_dem_diff.py` | Bridge LiDAR rasters (`generateBridgeRaster`, with `status()` check) and `BridgeDEMDiff` mosaicking. | -| `test_branchprocessing.py` | `BranchDerivation`, `AOIProcessingConfig`, `calculate_allbranches`, plus step-level tests for the BranchZero and CreateHAND substeps. | -| `test_calibrate_pipeline.py` | The full SRC calibration pipeline via one `run_calibration()` call with every `CalibrationConfig` parameter spelled out, plus one test per calibration stage. | -| `test_nwmstreamflow.py` | Streamflow retrieval (`getNWMretrospective`, `getNWMforecast`, `USGSData`), plotting, and KGE/NSE/PBias statistics. | -| `test_fimgeneration.py` | FIM generation from `discharge-inputs/` CSVs with date/range selection and depth output options. | -| `test_fimevaluation.py` | Benchmark FIM query/download via `queryBenchmarkFIM` (FIMbench) and candidate-vs-benchmark evaluation via `evaluateFIM` (FIMeval). | - -### Running -
- -```bash -# from the repo root, with the environment activated -pytest tests/ -v # whole suite -pytest tests/test_calibrate_pipeline.py -v # one stage -pytest tests/test_branchprocessing.py -v -k hand # one test by keyword -``` - -Before running, edit the constants at the top of each test file (for example `AOI_DIR`, `BANKFULL_FLOWS_FILE`, `N_WORKERS`) to match your machine. The calibration lookup tables referenced by the tests ship in the repo [`data/`](../data/) folder. - -The expected order when building an AOI from scratch mirrors the workflow: `test_getallinputdata` (stage inputs), `test_generate_dem_diff` (optional bridge healing), `test_branchprocessing` (HAND + SRC), `test_calibrate_pipeline` (calibration), `test_nwmstreamflow` (discharge), `test_fimgeneration` (flood maps), `test_fimevaluation` (benchmark evaluation). - -**For more usage notes refer to the [docs](../docs/) for the `fimbox` python package.** diff --git a/tests/conftest.py b/tests/conftest.py deleted file mode 100644 index 6d0e03c..0000000 --- a/tests/conftest.py +++ /dev/null @@ -1,5 +0,0 @@ -"""Test-suite-wide logging setup.""" - -from fimbox.logging_utils import configure_cli_logging - -configure_cli_logging() diff --git a/tests/test_branchprocessing.py b/tests/test_branchprocessing.py deleted file mode 100644 index 981aa40..0000000 --- a/tests/test_branchprocessing.py +++ /dev/null @@ -1,1008 +0,0 @@ -""" -Author: Supath Dhital -Date Updated: May 2026 - - -Branch processing tests. - -Run order: - 1. test_branch_derivation — level paths, branch polygons, branch list - 2. test_branch_zero_full — DEM clip, AGREE, pit-fill, D8 flowdir - 3. test_create_hand — full HAND generation (flow accum → split reaches) -""" - -import logging -from pathlib import Path - -# single steps IMPORTS -from fimbox import ( - BranchDerivation, -) - -log = logging.getLogger(__name__) - -# imports used only by the B-series CreateHAND step tests below. -# BranchZero substeps (StreamBooleanRasterizer, HydroenforceDEM, FlowdirDEM, -# HeadwaterRasterizer, LevelPathBooleanRasterizer, rasterize_3d_levee_lines, -# burn_levee_elevations) are exercised indirectly via test_step_Z1, so they -# are not imported here. - -# AOI parameters — point this at any user-supplied AOI working directory. -OUT_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" / "watershed-data" - -# Source-data filename prefix. -IDENTIFIER = "nwmmr" - -# Tunable CreateHAND parameters- All have sensible defaults in CreateHAND itself. -PARAMS_CREATE_HAND = dict( - cost_distance_tolerance=50.0, # m, lateral cost distance - lateral_elevation_threshold=10, # m, lateral thalweg drop cap - max_split_distance_m=1500.0, # m, split-reach max length - slope_min=0.0001, # rise/run floor - lakes_buffer_dist_m=100.0, # m, lake-boundary buffer - # SRC / crosswalk - mannings_n=0.06, # channel roughness - stage_min_m=0.0, # SRC stage ladder start - stage_interval_m=0.3048, # SRC stage step (1 ft) - stage_max_m=25.2984, # SRC stage ladder end (~83 ft) - min_catchment_area=0.25, # km^2, short-reach replace threshold - min_stream_length=0.5, # km, short-reach replace threshold - crosswalk_max_distance_m=100.0, # m, midpoint-to-NWM-flowline cap - # SRC slope source feeding Manning's equation: - # "iris_sword" (default) - IRIS-SWORD slope on order>=4 streams, else DEM - # "dem" - DEM rise/run slope only - # "hfab" - hydrofabric native slope, else DEM fallback - src_slope_source="iris_sword", - # IRIS-SWORD slope table (feature_id, slope_iris_sword). None -> the table - # shipped in fimbox/data is used when src_slope_source == "iris_sword". - iris_slope_csv=None, - # Hydrofabric slope column when it isn't the usual 'Slope'/'So'. - hfab_slope_column=None, -) - -DEM = OUT_DIR / "dem.tif" -STREAMS = OUT_DIR / f"{IDENTIFIER}_subset_streams.gpkg" -BOUNDARY_BUF = OUT_DIR / "wbd_buffered.gpkg" -CATCHMENTS = OUT_DIR / f"{IDENTIFIER}_catchments_proj_subset.gpkg" -HEADWATERS = ( - OUT_DIR / f"{IDENTIFIER}_headwater_points_subset.gpkg" - if (OUT_DIR / f"{IDENTIFIER}_headwater_points_subset.gpkg").is_file() - else OUT_DIR / f"{IDENTIFIER}_headwaters.gpkg" -) -LEVELPATH_EXT = OUT_DIR / f"{IDENTIFIER}_subset_streams_levelPaths_extended.gpkg" -BRIDGE_DIFF = OUT_DIR / "bridge_elev_diff.tif" -NLD_LEVEES = OUT_DIR / "3d_nld_subset_levees_burned.gpkg" - -# optional files -WBD8_CLP = OUT_DIR / "wbd8_clp.gpkg" -LAKES = OUT_DIR / f"{IDENTIFIER}_lakes_proj_subset.gpkg" -LEVEE_AREAS = OUT_DIR / "LeveeProtectedAreas_subset.gpkg" -LEVEE_LP_CSV = OUT_DIR / "levee_levelpaths.csv" - -# branch-zero derived paths -BRANCH_DIR = OUT_DIR / "branches" / "0" -BRANCH_ID = "0" -DEM_BRANCH = BRANCH_DIR / f"dem_{BRANCH_ID}.tif" -FLOWDIR = BRANCH_DIR / f"flowdir_d8_burned_filled_{BRANCH_ID}.tif" -HW_RASTER = BRANCH_DIR / f"headwaters_{BRANCH_ID}.tif" -STREAM_BOOL = BRANCH_DIR / f"flows_grid_boolean_{BRANCH_ID}.tif" -DEM_BURNED = BRANCH_DIR / f"dem_burned_{BRANCH_ID}.tif" -DEM_FILLED = BRANCH_DIR / f"dem_burned_filled_{BRANCH_ID}.tif" - -# REM + filtered catchment paths -REM = BRANCH_DIR / f"rem_{BRANCH_ID}.tif" -REM_ZEROED = BRANCH_DIR / f"rem_zeroed_masked_{BRANCH_ID}.tif" -CATCH_POLY = BRANCH_DIR / f"gw_catchments_reaches_{BRANCH_ID}.gpkg" -FILT_CATCH = ( - BRANCH_DIR / f"gw_catchments_reaches_filtered_addedAttributes_{BRANCH_ID}.gpkg" -) -FILT_FLOWS = BRANCH_DIR / f"demDerived_reaches_split_filtered_{BRANCH_ID}.gpkg" -FILT_TIF = ( - BRANCH_DIR / f"gw_catchments_reaches_filtered_addedAttributes_{BRANCH_ID}.tif" -) - -# SRC / crosswalk / hydroTable outputs (steps 16-21) -SLOPES_MASKED = BRANCH_DIR / f"slopes_d8_dem_meters_masked_{BRANCH_ID}.tif" -STAGE_TXT = BRANCH_DIR / f"stage_{BRANCH_ID}.txt" -CATCHLIST_TXT = BRANCH_DIR / f"catch_list_{BRANCH_ID}.txt" -SRC_BASE_CSV = BRANCH_DIR / f"src_base_{BRANCH_ID}.csv" -XWALK_CATCH = ( - BRANCH_DIR - / f"gw_catchments_reaches_filtered_addedAttributes_crosswalked_{BRANCH_ID}.gpkg" -) -XWALK_FLOWS = ( - BRANCH_DIR - / f"demDerived_reaches_split_filtered_addedAttributes_crosswalked_{BRANCH_ID}.gpkg" -) -SRC_FULL_CSV = BRANCH_DIR / f"src_full_crosswalked_{BRANCH_ID}.csv" -SRC_JSON = BRANCH_DIR / f"src_{BRANCH_ID}.json" -XWALK_CSV = BRANCH_DIR / f"crosswalk_table_{BRANCH_ID}.csv" -HYDRO_TABLE = BRANCH_DIR / f"hydroTable_{BRANCH_ID}.csv" -ROADS_CSV = BRANCH_DIR / f"osm_roads_fimpact_{BRANCH_ID}.csv" -BRIDGES_GPKG = BRANCH_DIR / f"osm_bridge_centroids_{BRANCH_ID}.gpkg" - -# HAND generation derived paths -FLOWACCUM = BRANCH_DIR / f"flowaccum_d8_burned_filled_{BRANCH_ID}.tif" -STREAM_PIX = BRANCH_DIR / f"demDerived_streamPixels_{BRANCH_ID}.tif" -THALWEG_ADJ = BRANCH_DIR / f"dem_lateral_thalweg_adj_{BRANCH_ID}.tif" -FLOWDIR_STR = BRANCH_DIR / f"flowdir_d8_burned_filled_flows_{BRANCH_ID}.tif" -THALWEG_COND = BRANCH_DIR / f"dem_thalwegCond_{BRANCH_ID}.tif" -SLOPES_D8 = BRANCH_DIR / f"slopes_d8_dem_{BRANCH_ID}.tif" -STREAM_ORDER = BRANCH_DIR / f"streamOrder_{BRANCH_ID}.tif" -SN_CATCH = BRANCH_DIR / f"sn_catchments_reaches_{BRANCH_ID}.tif" -DEM_REACHES = BRANCH_DIR / f"demDerived_reaches_{BRANCH_ID}.gpkg" -SPLIT_REACHES = BRANCH_DIR / f"demDerived_reaches_split_{BRANCH_ID}.gpkg" -SPLIT_PTS = BRANCH_DIR / f"demDerived_reaches_split_points_{BRANCH_ID}.gpkg" -GW_REACHES = BRANCH_DIR / f"gw_catchments_reaches_{BRANCH_ID}.tif" -PIXEL_PTS = BRANCH_DIR / f"flows_points_pixels_{BRANCH_ID}.gpkg" -GW_PIXELS = BRANCH_DIR / f"gw_catchments_pixels_{BRANCH_ID}.tif" - - -# ========================== -# COMBINED — the whole branch pipeline in one go, matching the step-by-step -# sequence exactly: -# Step Z0 BranchDerivation — level paths, branch polygons, branch_ids.lst -# Step Z1 BranchZero — whole-AOI DEM clip + AGREE + pit-fill + D8 -# (serial, in the main process, branch_id="0") -# Step B non-zero branches — BranchZero + CreateHAND per branch, in parallel -# -# Every parameter is spelled out so this test doubles as the parameter reference. -# Optional inputs (bridge_diff, levees, headwaters, levelpaths_extended) are -# resolved from OUT_DIR and passed only when the file exists on disk, matching -# the step-by-step behaviour. -# ============================ -def test_branchprocessing_combined(): - """Run the full branch pipeline in one call. - - Order inside calculate_allbranches: - 1. BranchZero for branch 0 — serial in main process, always first. - Branch 0 gets DEM clip / AGREE / pit-fill / D8 flowdir only. - CreateHAND does NOT run on branch 0; its flowdir is the shared - input every non-zero branch reads. After deny-list cleanup - branches/0/ keeps only what is NOT listed in deny_branch_zero.lst - (no hydroTable — that is expected and correct by design). - 2. All non-zero branches in parallel via Dask — BranchZero then the - full 22-step CreateHAND. Non-zero branches produce hydroTable. - 3. Deny-list cleanup removes intermediates from every branch dir. - """ - from fimbox import AOIProcessingConfig, calculate_allbranches - from fimbox._dask import _resolve_n_workers - - BranchDerivation( - out_dir=OUT_DIR, - branch_id_attribute="levpa_id", - reach_id_attribute="ID", - branch_buffer_distance_meters=7000.0, - ).run() - - bridge_diff = BRIDGE_DIFF if BRIDGE_DIFF.exists() else None - levee_gpkg = NLD_LEVEES if NLD_LEVEES.exists() else None - headwaters = HEADWATERS if HEADWATERS.exists() else None - levelpaths_extended = LEVELPATH_EXT if LEVELPATH_EXT.exists() else None - - n_workers = _resolve_n_workers() - cfg = AOIProcessingConfig( - aoi_dir=OUT_DIR, - branch_list_path=OUT_DIR / "branch_ids.lst", - # BranchZero inputs (whole-AOI, branch_id="0") - dem_path=DEM, - streams_gpkg=STREAMS, - boundary_gpkg=BOUNDARY_BUF, - bridge_elev_diff_path=bridge_diff, - levee_gpkg_path=levee_gpkg, - headwaters_gpkg=headwaters, - levelpaths_extended_gpkg=levelpaths_extended, - # AGREE DEM conditioning - agree_buffer_m=15.0, - agree_smooth_drop=10.0, - agree_sharp_drop=1000.0, - # CreateHAND geometry - cost_distance_tolerance=50.0, - lateral_elevation_threshold=10, - max_split_distance_m=1500.0, - slope_min=0.0001, - lakes_buffer_dist_m=100.0, - # SRC / crosswalk - mannings_n=0.06, - stage_min_m=0.0, - stage_interval_m=0.3048, - stage_max_m=25.2984, - min_catchment_area=0.25, - min_stream_length=0.5, - crosswalk_max_distance_m=100.0, - # SRC slope source: "iris_sword" | "dem" | "hfab" - src_slope_source="iris_sword", - iris_slope_csv=None, - hfab_slope_column=None, - # execution - n_workers=n_workers, - keep_failed_branches=True, - delete_deny_list=True, - ) - - result = calculate_allbranches( - cfg, - run_branch_zero=True, - delete_deny_list=True, - deny_unit_list=Path(__file__).resolve().parent.parent - / "config" - / "deny_unit.lst", - branch_ids_csv=OUT_DIR / "branch_ids.csv", - ) - - # Branch 0 now runs BranchZero + full CreateHAND (same as every non-zero branch). - b0 = OUT_DIR / "branches" / "0" - assert (b0 / "branch_zero_complete.txt").exists(), ( - "branch_zero_complete.txt missing" - ) - assert (b0 / "dem_0.tif").exists(), "dem_0.tif missing from branch 0" - assert (b0 / "flowdir_d8_burned_filled_0.tif").exists(), ( - "flowdir missing from branch 0" - ) - assert (b0 / "hydroTable_0.csv").exists(), "hydroTable_0.csv missing from branch 0" - assert result.n_branch_zero_recorded == 1, "branch zero not in branch_ids.csv" - assert result.branch_ids_csv.exists(), "branch_ids.csv not created" - - # branch_results now includes branch 0 at index 0. - b0_res = next((r for r in result.branch_results if r.branch_id == "0"), None) - assert b0_res is not None and b0_res.status == "ok", f"branch 0 status: {b0_res}" - - non_zero = [r for r in result.branch_results if r.branch_id != "0"] - ok = sum(1 for r in non_zero if r.status == "ok") - log.info(f"combined: branch_zero=ok non-zero ok={ok}/{len(non_zero)}") - assert result.n_non_zero_recorded == ok - - # Spot-check one non-zero branch hydroTable. - ok_branches = [r.branch_id for r in non_zero if r.status == "ok"] - if ok_branches: - sample_ht = ( - OUT_DIR / "branches" / ok_branches[0] / f"hydroTable_{ok_branches[0]}.csv" - ) - assert sample_ht.exists(), f"hydroTable missing from branch {ok_branches[0]}" - - -# ============================================================================= -# INDIVIDUAL STEP-BY-STEP TESTS -# Running these in file order rebuilds the full per-branch pipeline -# Layers: -# Z0 BranchDerivation — level paths + branch_list.csv -# Z1 BranchZero — DEM clip + AGREE + pit-fill + D8 -# (wraps stream raster, headwater -# raster, optional levelpath raster, -# optional levee burn, AGREE, -# pit-fill, flowdir) -# B02..B21 CreateHAND steps 2-21 — one test per CreateHAND step -# ============================================================================= - -# # Stage Z — bootstrap. Together they produce every input the B-series tests need. -# def test_step_Z0_branch_derivation(): -# """Derive level paths, branch polygons, and branch list from staged NWM data.""" -# result = BranchDerivation( -# out_dir=OUT_DIR, -# branch_id_attribute="levpa_id", -# reach_id_attribute="ID", -# branch_buffer_distance_meters=7000.0, -# ).run() -# assert result.dissolved_levelpaths.exists(), "dissolved levelpaths not written" -# assert result.branch_polygons.exists(), "branch polygons not written" -# assert result.branch_list.exists(), "branch list file not written" -# assert len(result.branch_dataframe) > 0, "branch dataframe is empty" -# log.info(f"branch count: {len(result.branch_dataframe)}") - - -# def test_step_Z1_branch_zero_full(): -# """Run BranchZero: DEM clip, stream rasterize, optional headwater/levelpath/levee -# rasters, AGREE conditioning, pit-fill, and D8 flowdir for branch 0. - -# This single call wraps the substeps BranchZero already folds together -# (StreamBooleanRasterizer, HeadwaterRasterizer, optional -# LevelPathBooleanRasterizer, optional rasterize_3d_levee_lines + -# burn_levee_elevations, HydroenforceDEM, WhiteboxTools pit-fill, -# FlowdirDEM). Calling the substeps individually would duplicate work the -# class already orchestrates correctly. -# """ -# outputs = BranchZero( -# dem_path=DEM, -# streams_gpkg=STREAMS, -# boundary_gpkg=BOUNDARY_BUF, -# out_dir=OUT_DIR, -# bridge_elev_diff_path=BRIDGE_DIFF if BRIDGE_DIFF.exists() else None, -# levee_gpkg_path=NLD_LEVEES if NLD_LEVEES.exists() else None, -# headwaters_gpkg=HEADWATERS if HEADWATERS.exists() else None, -# levelpaths_extended_gpkg=LEVELPATH_EXT if LEVELPATH_EXT.exists() else None, -# agree_buffer_m=15.0, -# agree_smooth_drop=10.0, -# agree_sharp_drop=1000.0, -# branch_zero_id=BRANCH_ID, -# ).run() -# for key, p in outputs.items(): -# log.info(f" {key:35s} --> {p.name}") -# assert DEM_BRANCH.exists(), "dem_0.tif missing" -# assert STREAM_BOOL.exists(), "flows_grid_boolean_0.tif missing" -# assert DEM_BURNED.exists(), "dem_burned_0.tif missing" -# assert DEM_FILLED.exists(), "dem_burned_filled_0.tif missing" -# assert FLOWDIR.exists(), "flowdir_d8_burned_filled_0.tif missing" - - -# # Stage B — CreateHAND steps 2..21, one isolated test each. -# def test_step_B02_flow_accumulation(): -# """CreateHAND step 2: D8 flow accumulation + stream-pixel mask.""" -# assert FLOWDIR.exists(), "FLOWDIR missing — run step_A6 first" -# if not HW_RASTER.exists(): -# log.warning("skipping flow accumulation — no headwater raster") -# return -# fa_out, sp_out = FlowAccDEM( -# flowdir=FLOWDIR, -# headwaters=HW_RASTER, -# out_flowaccum=FLOWACCUM, -# out_stream_pixels=STREAM_PIX, -# threshold=1.0, -# ).run() -# import rasterio - -# with rasterio.open(str(sp_out)) as src: -# stream_count = int((src.read(1) == 1).sum()) -# log.info(f"stream pixels: {stream_count}") -# assert fa_out.exists() and sp_out.exists() and stream_count > 0 - - -# def test_step_B03_thalweg_adjustment(): -# """CreateHAND step 3: lateral thalweg minimum + flow-conditioned DEM.""" -# for p in (DEM_BRANCH, STREAM_PIX, FLOWDIR): -# assert p.exists(), f"missing: {p}" -# result = ThalwegAdjustment( -# dem=DEM_BRANCH, -# stream_pixels=STREAM_PIX, -# flowdir=FLOWDIR, -# out_thalweg_adj=THALWEG_ADJ, -# out_flowdir_streams=FLOWDIR_STR, -# out_thalweg_cond=THALWEG_COND, -# cost_distance_tolerance=50.0, -# lateral_elevation_threshold=10, -# ).run() -# assert result["thalweg_adj"].exists() and result["thalweg_cond"].exists() - - -# def test_step_B04_d8_slopes(): -# """CreateHAND step 4: D8 slope raster (rise/run from thalweg-adjusted DEM).""" -# assert THALWEG_ADJ.exists() and FLOWDIR.exists() -# import numpy as np, rasterio - -# out = D8SlopeDEM( -# dem=THALWEG_ADJ, flowdir=FLOWDIR, out_path=SLOPES_D8, slope_min=0.0001 -# ).run() -# with rasterio.open(str(out)) as src: -# d = src.read(1) -# nd = src.nodata -# valid = d[(d != nd) & np.isfinite(d)] if nd is not None else d[np.isfinite(d)] -# log.info(f"slope range: [{valid.min():.6f}, {valid.max():.6f}]") -# # slope_min is clamped at 1e-4 in float32; allow a single-precision epsilon -# # of tolerance (~1e-7) so the test doesn't fail on the float32 representation -# # of 1e-4 (which is 9.9999997e-05). -# assert float(valid.min()) >= 0.0001 - 1e-7 - - -# def test_step_B05_streamnet_reaches(): -# """CreateHAND step 5: vectorise stream network into reach polylines.""" -# for p in (FLOWDIR, THALWEG_COND, FLOWACCUM, STREAM_PIX): -# assert p.exists(), f"missing: {p}" -# result = StreamNetReaches( -# flowdir=FLOWDIR, -# dem_thalweg_cond=THALWEG_COND, -# flowaccum=FLOWACCUM, -# stream_pixels=STREAM_PIX, -# out_dir=BRANCH_DIR, -# branch_id=BRANCH_ID, -# ).run() -# import geopandas as gpd - -# reaches = gpd.read_file(str(result["demDerived_reaches"])) -# log.info(f"reaches: {len(reaches)}") -# assert len(reaches) > 0 - - -# def test_step_B06_split_reaches(): -# """CreateHAND step 6: split reaches at length limit + lake boundaries.""" -# for p in (DEM_REACHES, THALWEG_COND, STREAMS): -# assert p.exists(), f"missing: {p}" -# split_gpkg, pts_gpkg = split_derived_reaches( -# reaches_gpkg=DEM_REACHES, -# dem_thalweg_cond=THALWEG_COND, -# nwm_streams_gpkg=STREAMS, -# out_split_gpkg=SPLIT_REACHES, -# out_points_gpkg=SPLIT_PTS, -# wbd8_clp_gpkg=WBD8_CLP if WBD8_CLP.exists() else None, -# lakes_gpkg=LAKES if LAKES.exists() else None, -# # This could be interesting point where based on slope or anyother logic, you can segment the reach--> ultimately gives the corresponsing -# # catchment, meaning shorter the reach length- denser the catchment -# max_length=1500.0, -# slope_min=0.0001, -# lakes_buffer_dist=100.0, -# ) -# import geopandas as gpd - -# split = gpd.read_file(str(split_gpkg)) -# log.info(f"split reaches: {len(split)} columns={list(split.columns)}") -# assert ( -# len(split) > 0 and "HydroID" in split.columns and "NextDownID" in split.columns -# ) - - -# def test_step_B07_gage_watershed_reaches(): -# """CreateHAND step 7: reverse-D8 walk labelling each pixel by its HydroID.""" -# from fimbox import GageCatchments - -# for p in (FLOWDIR, SPLIT_PTS): -# assert p.exists(), f"missing: {p}" -# # declutter=True mirrors CreateHAND step 7: solidify the reach raster -# # (fill pits, de-checkerboard, one piece per HydroID) so it polygonizes clean. -# GageCatchments( -# flowdir=FLOWDIR, -# outlet_points=SPLIT_PTS, -# out_path=GW_REACHES, -# declutter=True, -# ).run() -# assert GW_REACHES.exists() - - -# def test_step_B08_stream_pixel_points(): -# """CreateHAND step 8: vectorise stream-pixel centroids (one point per stream pixel).""" -# from fimbox import stream_pixel_points - -# assert STREAM_PIX.exists() -# stream_pixel_points(stream_pixels=STREAM_PIX, out_gpkg=PIXEL_PTS) -# assert PIXEL_PTS.exists() - - -# def test_step_B09_gage_watershed_pixels(): -# """CreateHAND step 9: reverse-D8 walk labelling each pixel by NWM feature_id.""" -# from fimbox import GageCatchments - -# for p in (FLOWDIR, PIXEL_PTS): -# assert p.exists(), f"missing: {p}" -# GageCatchments( -# flowdir=FLOWDIR, -# outlet_points=PIXEL_PTS, -# out_path=GW_PIXELS, -# ).run() -# assert GW_PIXELS.exists() - - -# def test_step_B10_outlet_backpool_mitigation(): -# """CreateHAND step 10: trim oversized outlet catchments (no-op for branch 0).""" -# from fimbox import OutletBackpoolMitigate - -# for p in (SPLIT_REACHES, GW_PIXELS, GW_REACHES, SPLIT_PTS, STREAMS, THALWEG_COND): -# assert p.exists(), f"missing: {p}" -# OutletBackpoolMitigate( -# branch_dir=BRANCH_DIR, -# catchment_pixels_path=GW_PIXELS, -# catchment_reaches_path=GW_REACHES, -# split_flows_gpkg=SPLIT_REACHES, -# split_points_gpkg=SPLIT_PTS, -# nwm_streams_gpkg=STREAMS, -# dem_path=THALWEG_COND, -# slope_min=0.0001, -# ).run() -# # No new file is asserted — backpool mitigation modifies the existing -# # gw_catchments_pixels/reaches rasters in place for non-zero branches only. -# assert GW_PIXELS.exists() and GW_REACHES.exists() - - -# def test_step_B11_make_rem(): -# """CreateHAND step 11: HAND = pixel_elev - nearest_stream_pixel_elev. - -# Note: the raw REM **can** be negative (pixels lower than the nearest -# downstream stream pixel — happens near floodplain edges and where the -# D8 walk crosses meander cutoffs). Negative values get clipped to zero -# in step 12 (``rem_zeroed_masked``). This test only asserts the raster -# was produced and contains finite values — it does NOT enforce -# non-negativity, which is a step-12 invariant. -# """ -# from fimbox import MakeREM - -# for p in (THALWEG_COND, GW_PIXELS, STREAM_PIX): -# assert p.exists(), f"missing: {p}" -# out = MakeREM( -# dem_thalweg_cond=THALWEG_COND, -# gw_catchments_pixels=GW_PIXELS, -# stream_pixels=STREAM_PIX, -# out_rem=REM, -# ).run() -# import rasterio, numpy as np - -# with rasterio.open(str(out)) as src: -# data = src.read(1) -# nd = src.nodata -# valid = data[data != nd] if nd is not None else data.ravel() -# log.info( -# f"REM range: [{float(valid.min()):.2f}, {float(valid.max()):.2f}] " -# f"({(valid < 0).sum()} negative pixels — clipped by step 12)" -# ) -# assert out.exists() and valid.size > 0 and np.isfinite(valid).all() - - -# def test_step_B11b_rem_nonnegative_after_zero_mask(): -# """Cross-check: after step 12 (rem_zeroed_masked), the REM raster must be -# non-negative and contain no NaN pixels. The reference formula -# ``(A * (A>=0) * (B>0))`` with an explicit NoDataValue treats NaN inputs as -# zero; the fimbox port now matches that behaviour by rewriting NaN to the -# nodata sentinel before the multiply. - -# Lives next to B11 so a failure here points at the zero-mask logic, not at -# MakeREM itself. Skipped silently if step 12 hasn't run yet (run B12 first). -# """ -# import numpy as np -# import rasterio - -# if not REM_ZEROED.exists(): -# log.warning("skipping non-negativity check — run step_B12 first") -# return -# with rasterio.open(str(REM_ZEROED)) as src: -# data = src.read(1) -# nd = src.nodata -# # Strip both the nodata sentinel and any NaN before the min() so the -# # test catches the actual data range, not an IEEE NaN propagating. -# if nd is not None: -# valid_mask = (data != nd) & ~np.isnan(data) -# else: -# valid_mask = ~np.isnan(data) -# valid = data[valid_mask] -# nan_count = int(np.isnan(data).sum()) -# log.info(f"REM zero-mask: {valid.size} valid pixels, {nan_count} NaN pixels") -# assert valid.size > 0 -# assert ( -# nan_count == 0 -# ), f"step 12 leaked {nan_count} NaN pixels into the masked REM raster" -# assert ( -# float(valid.min()) >= 0.0 -# ), f"step 12 left negatives in REM: min={valid.min()}" - - -# def test_step_B12_rem_zeroed_masked(): -# """CreateHAND step 12: clip negative HAND to 0 + mask outside catchments.""" -# from fimbox import rem_zeroed_masked - -# for p in (REM, GW_REACHES): -# assert p.exists(), f"missing: {p}" -# rem_zeroed_masked(REM, GW_REACHES, REM_ZEROED) -# assert REM_ZEROED.exists() - - -# def test_step_B13_polygonize_catchments(): -# """CreateHAND step 13: rasterised catchments --> per-HydroID polygon gpkg.""" -# # Helper lives inside create_hand.py as a private function; import it explicitly. -# from fimbox.preprocessing.calculate_branch.create_hand import ( -# _polygonize_catchments, -# ) - -# assert GW_REACHES.exists() -# _polygonize_catchments(GW_REACHES, CATCH_POLY) -# import geopandas as gpd - -# gdf = gpd.read_file(str(CATCH_POLY)) -# log.info(f"polygonised: {len(gdf)} catchments") -# assert CATCH_POLY.exists() and "HydroID" in gdf.columns and len(gdf) > 0 - - -# def test_step_B14_filter_catchments(): -# """CreateHAND step 14: drop slivers + attach flow attributes per HydroID.""" -# from fimbox import FilterCatchments - -# for p in (CATCH_POLY, SPLIT_REACHES): -# assert p.exists(), f"missing: {p}" -# out_catch, out_flows = FilterCatchments( -# catchments_gpkg=CATCH_POLY, -# flows_gpkg=SPLIT_REACHES, -# out_catchments=FILT_CATCH, -# out_flows=FILT_FLOWS, -# aoi_code=OUT_DIR.parent.name, -# boundary_gpkg=WBD8_CLP if WBD8_CLP.exists() else None, -# ).run() -# import geopandas as gpd - -# catches = gpd.read_file(str(out_catch)) -# flows = gpd.read_file(str(out_flows)) -# log.info(f"filtered catchments: {len(catches)} flows: {len(flows)}") -# assert len(catches) > 0 and "areasqkm" in catches.columns -# assert len(flows) > 0 and "HydroID" in flows.columns - - -# def test_step_B15_rasterize_filtered_catchments(): -# """CreateHAND step 15: burn HydroID back onto the reference raster grid.""" -# from fimbox.preprocessing.calculate_branch.create_hand import ( -# _rasterize_catchments, -# ) - -# for p in (FILT_CATCH, GW_REACHES): -# assert p.exists(), f"missing: {p}" -# _rasterize_catchments(FILT_CATCH, GW_REACHES, FILT_TIF) -# assert FILT_TIF.exists() - - -# def test_step_B16_mask_slopes_to_catchments(): -# """CreateHAND step 16: clip D8 slopes to the filtered catchment mask.""" -# from fimbox import mask_slopes_to_catchments - -# for p in (SLOPES_D8, FILT_TIF): -# assert p.exists(), f"missing: {p}" -# mask_slopes_to_catchments(SLOPES_D8, FILT_TIF, SLOPES_MASKED) -# assert SLOPES_MASKED.exists() - - -# def test_step_B17_stages_and_catchlist(): -# """CreateHAND step 17: write the stage ladder + per-HydroID metadata text files.""" -# from fimbox import make_stages_and_catchlist - -# for p in (FILT_FLOWS, FILT_CATCH): -# assert p.exists(), f"missing: {p}" -# make_stages_and_catchlist( -# flows_gpkg=FILT_FLOWS, -# catchments_gpkg=FILT_CATCH, -# out_stages=STAGE_TXT, -# out_catchlist=CATCHLIST_TXT, -# stages_min=0.0, -# stages_interval=0.3048, -# stages_max=25.2984, -# ) -# assert STAGE_TXT.exists() and CATCHLIST_TXT.exists() - - -# def test_step_B18_build_src_base(): -# """CreateHAND step 18: synthetic rating curve base table.""" -# from fimbox import build_src_base - -# for p in (REM_ZEROED, FILT_TIF, SLOPES_MASKED, CATCHLIST_TXT, STAGE_TXT): -# assert p.exists(), f"missing: {p}" -# build_src_base( -# hand_raster=REM_ZEROED, -# catch_raster=FILT_TIF, -# slope_raster=SLOPES_MASKED, -# catchlist_txt=CATCHLIST_TXT, -# stages_txt=STAGE_TXT, -# out_csv=SRC_BASE_CSV, -# ) -# import pandas as pd - -# df = pd.read_csv(SRC_BASE_CSV) -# log.info(f"src_base: {len(df)} rows HydroIDs={df['CatchId'].nunique()}") -# assert SRC_BASE_CSV.exists() and len(df) > 0 - - -# def test_step_B19_add_crosswalk(): -# """CreateHAND step 19: NWM crosswalk + Manning's hydraulics + hydroTable.""" -# from fimbox import add_crosswalk - -# for p in (FILT_CATCH, FILT_FLOWS, SRC_BASE_CSV, STREAMS): -# assert p.exists(), f"missing: {p}" -# add_crosswalk( -# catchments_gpkg=FILT_CATCH, -# flows_gpkg=FILT_FLOWS, -# src_base_csv=SRC_BASE_CSV, -# nwm_streams_gpkg=STREAMS, -# out_catchments_gpkg=XWALK_CATCH, -# out_flows_gpkg=XWALK_FLOWS, -# out_src_csv=SRC_FULL_CSV, -# out_src_json=SRC_JSON, -# out_crosswalk_csv=XWALK_CSV, -# out_hydro_csv=HYDRO_TABLE, -# boundary_gpkg=WBD8_CLP if WBD8_CLP.exists() else None, -# mannings_n=0.06, -# min_catchment_area=0.25, -# min_stream_length=0.5, -# max_distance_m=100.0, -# small_segments_csv=BRANCH_DIR / f"small_segments_{BRANCH_ID}.csv", -# # SRC slope source (optional): "iris_sword" (default) | "dem" | "hfab". -# src_slope_source="iris_sword", -# iris_slope_csv=None, # None -> packaged fimbox/data table -# hfab_slope_column=None, # name the hydrofabric slope col if not Slope/So -# ) -# import pandas as pd - -# ht = pd.read_csv(HYDRO_TABLE, dtype={"HydroID": str}) -# log.info(f"hydroTable: {len(ht)} rows HydroIDs={ht['HydroID'].nunique()}") -# assert HYDRO_TABLE.exists() and (ht["discharge_cms"] >= 0).all() -# # The three slope variants are carried so the chosen source is transparent. -# for col in ("SLOPE", "SLOPE_RISE_RUN", "SLOPE_IRIS_SWORD"): -# assert col in ht.columns, f"hydroTable missing {col}" - - -# def test_step_B20_heal_bridges_osm(): -# """CreateHAND step 20: raise HAND at OSM bridge decks (in-place REM update).""" -# from fimbox import heal_bridges_osm - -# bridges_gpkg = OUT_DIR / "osm_bridges_subset.gpkg" -# if not bridges_gpkg.exists(): -# log.warning("skipping bridge heal — no OSM bridges gpkg") -# return -# for p in (REM_ZEROED, XWALK_CATCH): -# assert p.exists(), f"missing: {p}" -# bridge_diff = OUT_DIR / "bridge_elev_diff.tif" -# heal_bridges_osm( -# hand_raster=REM_ZEROED, -# bridges_gpkg=bridges_gpkg, -# catchments_gpkg=XWALK_CATCH, -# out_centroids_gpkg=BRIDGES_GPKG, -# bridge_diff_raster=bridge_diff if bridge_diff.exists() else None, -# ) -# assert BRIDGES_GPKG.exists() - - -# def test_step_B21_process_roads_fimpact(): -# """CreateHAND step 21: sample HAND along OSM roads to derive flood thresholds.""" -# from fimbox import process_roads_fimpact - -# roads_gpkg = OUT_DIR / "osm_roads_subset.gpkg" -# if not roads_gpkg.exists(): -# log.warning("skipping road FIMpact — no OSM roads gpkg") -# return -# for p in (REM_ZEROED, XWALK_CATCH): -# assert p.exists(), f"missing: {p}" -# process_roads_fimpact( -# hand_raster=REM_ZEROED, -# roads_gpkg=roads_gpkg, -# catchments_gpkg=XWALK_CATCH, -# out_csv=ROADS_CSV, -# ) -# assert ROADS_CSV.exists() - - -# # Stage C — branch-zero post-CreateHAND steps -# # (download USGS gauges --> AOI-level assignment --> branch-zero crosswalk --> cleanup) - -# # AOI-level path to the staged USGS gages gpkg -# USGS_GAGES = OUT_DIR / "usgs_gages.gpkg" -# USGS_SUBSET = OUT_DIR / "usgs_subset_gages.gpkg" -# USGS_SUBSET_BZERO = OUT_DIR / f"usgs_subset_gages_{BRANCH_ID}.gpkg" -# NWM_LEVELPATHS = OUT_DIR / f"{IDENTIFIER}_subset_streams_levelPaths.gpkg" - - -# def test_step_C20_download_usgs_gages(): -# """Download USGS gauge points inside the AOI from the ArcGIS Online -# FeatureServer. Writes ``usgs_gages.gpkg`` at the AOI root, with the columns -# ``assign_gages_to_branches`` expects: ``location_id``, ``feature_id``, -# ``aoi_id``, ``source``, geometry. -# """ -# from fimbox import DownloadUSGSGages - -# # Use the buffered boundary so gauges just outside the WBD are still -# # captured (they may snap to streams that drain into the AOI). -# boundary = BOUNDARY_BUF if BOUNDARY_BUF.exists() else WBD8_CLP -# assert boundary.exists(), f"missing boundary: {boundary}" - -# gdf = DownloadUSGSGages().download( -# boundary=boundary, -# aoi_id=OUT_DIR.parent.name, -# out_dir=OUT_DIR, -# out_name="usgs_gages.gpkg", -# out_layer="usgs_gages", -# ) -# log.info(f"USGS gauges downloaded: {len(gdf)} features --> {USGS_GAGES.name}") -# # Empty AOI (no gauges in CONUS layer) is acceptable; only assert the -# # file exists when at least one feature came back. -# if len(gdf) > 0: -# assert USGS_GAGES.exists() -# assert {"location_id", "feature_id", "aoi_id", "source"}.issubset(gdf.columns) - - -# def test_step_C21_assign_gages_to_branches(): -# """Stage 1 of the gage crosswalk: tag every gage with a ``feature_id`` + -# ``levpa_id`` (= branch id) and write the AOI-wide + branch-zero gpkgs. - -# Skips if either ``usgs_gages.gpkg`` (from C20) or -# ``nwm_subset_streams_levelPaths.gpkg`` (from BranchDerivation in Z0) is -# missing — both prerequisites get logged so a failure points at the -# right upstream step. -# """ -# from fimbox import assign_gages_to_branches - -# if not USGS_GAGES.exists(): -# log.warning( -# "skipping gage assignment — usgs_gages.gpkg missing (run step_C20 first)" -# ) -# return -# if not NWM_LEVELPATHS.exists(): -# log.warning( -# "skipping gage assignment — nwm_subset_streams_levelPaths.gpkg missing " -# "(run step_Z0_branch_derivation first)" -# ) -# return - -# assign_gages_to_branches( -# usgs_gages_gpkg=USGS_GAGES, -# nwm_streams_levelpaths_gpkg=NWM_LEVELPATHS, -# aoi_id=OUT_DIR.parent.name, -# out_dir=OUT_DIR, -# # DownloadUSGSGages writes "aoi_id"; the default filter column ("HUC8") -# # would not find anything in that gpkg. -# aoi_filter_column="aoi_id", -# branch_zero_id=BRANCH_ID, -# ) -# # When the AOI actually contains gauges both files exist; on empty AOIs -# # neither is written and the function returns None (logged a warning). -# if USGS_SUBSET.exists(): -# log.info( -# f"AOI-wide gages --> {USGS_SUBSET.name} | " -# f"branch-zero --> {USGS_SUBSET_BZERO.name}" -# ) -# assert USGS_SUBSET_BZERO.exists() - - -# def test_step_C22_usgs_crosswalk_branch_zero(): -# """Stage 2 of the gage crosswalk for branch zero. - -# Snaps every branch-zero gage to its DEM-derived thalweg and samples the -# DEM + thalweg-conditioned DEM to populate ``dem_elevation`` and -# ``dem_adj_elevation`` on the gage table. Output: -# ``branches/0/usgs_elev_table.csv``. - -# Prerequisites: ``usgs_subset_gages_0.gpkg`` (from C21) and the per-branch -# CreateHAND outputs (from Z1 + the B-series). -# """ -# from fimbox import run_branch_crosswalk - -# if not USGS_SUBSET_BZERO.exists(): -# log.warning( -# "skipping USGS crosswalk — usgs_subset_gages_0.gpkg missing " -# "(run step_C20 + step_C21 first to produce it)" -# ) -# return -# bzero_gages = USGS_SUBSET_BZERO - -# # dem_meters_{B}.tif is the inundation-mapping name; fimbox writes dem_{B}.tif -# # via BranchZero. Use whichever exists. -# dem_b = BRANCH_DIR / f"dem_meters_{BRANCH_ID}.tif" -# if not dem_b.exists(): -# dem_b = DEM_BRANCH - -# for p in (XWALK_CATCH, FILT_FLOWS, dem_b, THALWEG_COND): -# assert p.exists(), f"missing: {p}" - -# out = run_branch_crosswalk( -# aoi_gages_gpkg=bzero_gages, -# branch_catchments_gpkg=XWALK_CATCH, -# branch_flows_gpkg=FILT_FLOWS, -# dem_path=dem_b, -# dem_thalweg_path=THALWEG_COND, -# branch_id=BRANCH_ID, -# out_dir=BRANCH_DIR, -# ) -# usgs_table = BRANCH_DIR / "usgs_elev_table.csv" -# log.info(f"USGS crosswalk wrote: {[p for p in out.values() if p]}") -# # usgs_elev_table.csv only exists when the AOI has gages — log either way. -# if usgs_table.exists(): -# import pandas as pd - -# df = pd.read_csv(usgs_table) -# log.info(f"usgs_elev_table.csv rows: {len(df)}") - - -# def test_step_C23_outputs_cleanup_branch_zero(): -# """Apply the deny-list cleanup to ``branches/0/``. - -# Default behaviour deletes every intermediate raster + vector listed in -# --> fimbox/config/deny_branch_zero.lst. -# """ -# import os - -# from fimbox import remove_deny_list_files - -# deny_path = ( -# Path(__file__).resolve().parent.parent / "config" / "deny_branch_zero.lst" -# ) -# assert deny_path.is_file(), f"deny list missing: {deny_path}" - -# # API sanity checks that always run (never touch real files). -# assert remove_deny_list_files(BRANCH_DIR, "NONE", BRANCH_ID) == 0 -# assert remove_deny_list_files(BRANCH_DIR, "none", BRANCH_ID) == 0 - -# if os.environ.get("FIMBOX_KEEP_BRANCH_ZERO"): -# n_patterns = sum( -# 1 -# for L in deny_path.read_text().splitlines() -# if L.strip() and not L.lstrip().startswith("#") -# ) -# log.info( -# f"step_C23: skipping cleanup (FIMBOX_KEEP_BRANCH_ZERO set). " -# f"{deny_path.name} has {n_patterns} active patterns; " -# "unset the env var to enable cleanup." -# ) -# return - -# # The branch-0 directory may be empty when only later steps have been -# # populated, or when an earlier C23 run already cleaned it. Skip cleanly -# # if there's nothing to do. -# if not BRANCH_DIR.exists(): -# log.warning(f"skipping cleanup — branch dir {BRANCH_DIR} missing") -# return - -# n = remove_deny_list_files( -# src_dir=BRANCH_DIR, -# deny_list=deny_path, -# branch_id=BRANCH_ID, -# verbose=True, -# ) -# log.info(f"step_C23: removed {n} files from {BRANCH_DIR}") - - -# def test_step_C24_calculate_allbranches(tmp_path): -# """Fast wrapper check without launching real branch workers.""" -# from fimbox import AOIProcessingConfig, calculate_allbranches - -# aoi_dir = tmp_path / "aoi" -# aoi_dir.mkdir() -# # Match BranchDerivation's actual output: branch_ids.lst (one id per line). -# # Empty file = branch-zero-only run, which is what this wrapper test exercises. -# branch_list_path = aoi_dir / "branch_ids.lst" -# branch_list_path.write_text("") - -# deny_unit_list = tmp_path / "deny_unit.lst" -# deny_unit_list.write_text("temporary_{}.tif\n") -# # aoi_id defaults to the AOI folder name ("aoi") when not passed. -# removable = aoi_dir / f"temporary_{aoi_dir.name}.tif" -# removable.write_bytes(b"x") - -# result = calculate_allbranches( -# AOIProcessingConfig( -# aoi_dir=aoi_dir, -# branch_list_path=branch_list_path, -# n_workers=1, -# ), -# delete_deny_list=True, -# deny_unit_list=deny_unit_list, -# branch_ids_csv=aoi_dir / "branch_ids.csv", -# ) - -# assert result.n_branch_zero_recorded == 1 -# assert result.n_non_zero_recorded == 0 -# assert result.n_unit_files_removed == 1 -# assert not removable.exists() - - -# def test_step_C25_calculate_allbranches_live_run(): -# """Live run for the real non-zero branch loop. - -# Set FIMBOX_KEEP_UNIT=1 to skip AOI-level cleanup. -# Set FIMBOX_SKIP_ALLBRANCHES=1 to skip this test (e.g. during quick CI -# smoke runs); by default it always runs. -# """ -# from fimbox import AOIProcessingConfig, calculate_allbranches - -# if os.environ.get("FIMBOX_SKIP_ALLBRANCHES"): -# pytest.skip("FIMBOX_SKIP_ALLBRANCHES set — skipping live branch loop") - -# # BranchDerivation writes branch_ids.lst -# branch_list_path = OUT_DIR / "branch_ids.lst" - -# deny_unit_list = Path(__file__).resolve().parent.parent / "config" / "deny_unit.lst" -# assert deny_unit_list.is_file(), f"deny_unit.lst missing: {deny_unit_list}" - -# # Reuse branch-zero tuning so both paths stay in sync. -# # Auto-size workers to the machine; set FIMBOX_DASK_WORKERS=1 for serial. -# from fimbox._dask import _resolve_n_workers - -# n_workers = _resolve_n_workers() -# log.info(f"Branch processing with n_workers={n_workers}") - -# cfg = AOIProcessingConfig( -# aoi_dir=OUT_DIR, -# branch_list_path=branch_list_path, -# n_workers=n_workers, # auto-sized; FIMBOX_DASK_WORKERS=1 forces serial -# keep_failed_branches=True, # keep a failed branch dir for inspection -# delete_deny_list=True, -# **PARAMS_CREATE_HAND, -# ) - -# delete_deny_list = True -# result = calculate_allbranches( -# cfg, -# delete_deny_list=delete_deny_list, -# deny_unit_list=deny_unit_list if delete_deny_list else None, -# branch_ids_csv=OUT_DIR / "branch_ids.csv", -# ) - -# assert result.n_branch_zero_recorded == 1 -# assert result.branch_ids_csv.exists(), "branch_ids.csv was not created" -# assert result.n_non_zero_recorded == sum( -# 1 for r in result.branch_results if r.status == "ok" -# ) diff --git a/tests/test_calibrate_pipeline.py b/tests/test_calibrate_pipeline.py deleted file mode 100644 index 7212435..0000000 --- a/tests/test_calibrate_pipeline.py +++ /dev/null @@ -1,289 +0,0 @@ -""" -Author: Supath Dhital -Date Updated: June 2026 - -Tests for the synthetic rating curve (SRC) calibration pipeline. - -Two layers: - - COMBINED ......... test_calibrate_full_pipeline runs the whole thing in a - single run_calibration() call against the live AOI, with EVERY optional - CalibrationConfig parameter spelled out so the full surface is visible - in one place. - - STEP BY STEP ..... one test per stage (thalweg, longitudinal, bathymetry, - bankfull, subdiv, nonmonotonic, usgs, spatial, log scan) so any single - step can be run / debugged alone. - -It will point into the working version of the AOI and skip when it is absent. -""" - -from __future__ import annotations - -from pathlib import Path - -import pandas as pd -import pytest - -from fimbox import CalibrationConfig, run_calibration -from fimbox._dask import _resolve_n_workers -from fimbox.datasets import fetch_data - -# Live AOI + input files. Edit these to point at your data; tests skip when the AOI is absent. -AOI_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" - -# Calibration lookup tables, fetched on demand from the public SDML S3 bucket -# (anonymous, cached locally) via the pooch registry in ``fimbox.datasets``. - -# Bankfull recurrence flows (NWM v3) -BANKFULL_FLOWS_FILE = fetch_data("nwm3_high_water_threshold") - -# Optimized variable-roughness Manning's n table (per feature_id channel/overbank n) -VMANN_INPUT_FILE = fetch_data("mannings_optz") - -# USGS rating-curve calibration. Rating curve + NWM recurrence flows (v3) are -# required; the acceptable-gage quality filter refines which gages qualify. -USGS_RATING_CURVE_CSV = fetch_data("usgs_rating_curves") -NWM_RECUR_FILE = fetch_data("nwm3_recurrence_flows") -USGS_ACCEPTABLE_GAGES = fetch_data("acceptable_gages") - -# Bathymetry: eHydro surveyed channels (.gpkg). -BATHY_EHYDRO_FILE = fetch_data("bathymetry_ehydro_ohrfc") - -# Spatial-observation calibration: per-AOI benchmark points (.parquet). Yet to -# be added — left unset for now. -CALIB_POINTS_FILE = None - -# Manual calibration: per-feature_id coefficient CSV. Yet to be added. -MAN_CALB_FILE = None - -# Worker count for the branch-parallel routines — auto-sized to the device -JOB_BRANCH_LIMIT = _resolve_n_workers() - -_BRANCHES = ( - AOI_DIR / "watershed-data" / "branches" - if (AOI_DIR / "watershed-data" / "branches").is_dir() - else AOI_DIR / "branches" -) -_skip_no_aoi = pytest.mark.skipif( - not _BRANCHES.is_dir(), reason=f"AOI not present: {_BRANCHES}" -) -_skip_no_bankfull = pytest.mark.skipif( - not BANKFULL_FLOWS_FILE.is_file(), reason=f"missing {BANKFULL_FLOWS_FILE}" -) -_skip_no_bathy = pytest.mark.skipif( - BATHY_EHYDRO_FILE is None or not Path(BATHY_EHYDRO_FILE).is_file(), - reason=f"bathy eHydro file not set: {BATHY_EHYDRO_FILE}", -) -_skip_no_usgs = pytest.mark.skipif( - not USGS_RATING_CURVE_CSV.is_file(), reason=f"missing {USGS_RATING_CURVE_CSV}" -) -_skip_no_manual = pytest.mark.skipif( - MAN_CALB_FILE is None or not Path(MAN_CALB_FILE).is_file(), - reason=f"manual calib file not set: {MAN_CALB_FILE}", -) -_skip_no_spatial = pytest.mark.skipif( - CALIB_POINTS_FILE is None or not Path(CALIB_POINTS_FILE).is_file(), - reason=f"spatial calib points not set: {CALIB_POINTS_FILE}", -) - - -# COMBINED — the whole calibration pipeline in one call, matching the step-by-step -# sequence: reset -> aggregate_pre -> thalweg -> longitudinal -> bathymetry -> -# bankfull -> subdiv -> nonmonotonic -> usgs -> spatial -> manual -> aggregate_post -> log_scan. -# File-dependent steps (bathy, usgs, spatial, manual) self-skip when their -# input file is absent, matching the step-by-step skip decorators. -@_skip_no_aoi -@_skip_no_bankfull -def test_calibrate_full_pipeline(): - """One run_calibration() call driving the full default pipeline. - Every CalibrationConfig parameter is spelled out, grouped by step.""" - cfg = CalibrationConfig( - # reset — revert hydroTables to uncalibrated baseline before re-applying. - # Set True when re-calibrating an AOI that was already calibrated. - calibration_rerun=True, - # aggregate_pre — assemble usgs/ras2fim elev tables before adjustments - aggregate_pre=True, - # thalweg — remove thalweg-notch artifact rows, refill stage ladder - thalweg_notches_adjustment=True, - # longitudinal — smooth hydraulic geometry along reach chains - longitudinal_filter=True, - # bathymetry — add missing in-channel area below the DEM (needs bathy_file_ehydro) - bathymetry_adjust=True, - bathy_file_ehydro=BATHY_EHYDRO_FILE, - # bankfull — identify bankfull stage in every branch SRC - src_bankfull_toggle=True, - bankfull_flows_file=BANKFULL_FLOWS_FILE, - include_branch_zero=True, - # subdiv — channel/overbank subdivision (needs vmann + bankfull on) - src_subdiv_toggle=True, - vmann_input_file=VMANN_INPUT_FILE, - default_channel_n=0.06, # used when feature_id missing from vmann table - default_overbank_n=0.12, - # nonmonotonic — force monotonic in-channel rating curves - nonmonotonic_src_adjustment=True, - nonmonotonic_stream_order_min=4, - # usgs — calibrate SRCs against USGS rating curves at NWM recurrence flows - src_adjust_usgs=True, - usgs_rating_curve_csv=USGS_RATING_CURVE_CSV, - usgs_acceptable_gages=USGS_ACCEPTABLE_GAGES, - nwm_recur_file=NWM_RECUR_FILE, - # spatial — calibrate SRCs against benchmark inundation points - src_adjust_spatial=True, - calib_points_file=CALIB_POINTS_FILE, # None -> step self-skips - # manual — apply a per-feature_id coefficient table - manual_calb_toggle=True, - man_calb_file=MAN_CALB_FILE, # None -> step self-skips - # aggregate_post — publish htable + bridge + road to AOI root - aggregate_post=True, - # log scan — collect error/warning lines into per-AOI summary files - scan_logs=True, - # execution - job_branch_limit=JOB_BRANCH_LIMIT, - skip_unimplemented=True, # warn instead of raising on stubs - ) - run_calibration(AOI_DIR, cfg) - - # Subdivision rewrites the per-branch hydroTable with subdiv columns. - sample_ht = next(_BRANCHES.glob("*/hydroTable_*.csv")) - cols = pd.read_csv(sample_ht, nrows=1).columns - assert "subdiv_discharge_cms" in cols - assert "channel_n" in cols - - -# # STEP BY STEP — each stage on its own. -# @_skip_no_aoi -# def test_step_reset(): -# """Reset per-branch hydroTable + src_full_crosswalked to baseline. -# Needed only for reruns; a no-op on a fresh AOI. Runs before aggregation.""" -# HydroTableReset(aoi_dir=AOI_DIR).run() - - -# @_skip_no_aoi -# def test_step_aggregate_pre(): -# """Pre-calibration aggregation: usgs/ras2fim elev tables if available (not integrated yet) -> AOI root.""" -# BranchAggregator(aoi_dir=AOI_DIR, usgs_elev=True, ras_elev=True).run() - - -# @_skip_no_aoi -# def test_step_thalweg_notches(): -# """Remove thalweg-notch artifact rows and refill the stage ladder.""" -# results = ThalwegNotchesAdjustment( -# aoi_dir=AOI_DIR, -# n_workers=JOB_BRANCH_LIMIT, # branch-parallel -# stage_interval_m=0.3048, # SRC stage step -# n_stages=84, # full ladder length -# extrap_rows=3, # trailing rows fit for extrapolation -# ).run() -# assert results - - -# @_skip_no_aoi -# def test_step_longitudinal(): -# """Smooth hydraulic geometry along reach chains, recompute discharge.""" -# results = LongitudinalFlowFilter( -# aoi_dir=AOI_DIR, n_workers=JOB_BRANCH_LIMIT, n_stages=84 -# ).run() -# assert results - - -# @_skip_no_aoi -# @_skip_no_bathy -# def test_step_bathymetry(): -# """Add missing in-channel area below the DEM from eHydro surveys, then -# recompute discharge.""" -# results = BathymetricAdjustment( -# aoi_dir=AOI_DIR, -# bathy_file_ehydro=BATHY_EHYDRO_FILE, -# ).run() -# assert results - - -# @_skip_no_aoi -# @_skip_no_bankfull -# def test_step_bankfull(): -# """Identify bankfull stage in every branch SRC.""" -# results = SrcBankfull( -# aoi_dir=AOI_DIR, -# bankfull_flows_file=BANKFULL_FLOWS_FILE, -# n_workers=JOB_BRANCH_LIMIT, -# include_branch_zero=True, -# ).run() -# assert results # dict of branch_id -> status string - - -# @_skip_no_aoi -# @_skip_no_bankfull -# def test_step_subdiv(): -# """Channel/overbank subdivision. Depends on bankfull having run, so run -# it first within this test to keep the step self-contained.""" -# SrcBankfull( -# aoi_dir=AOI_DIR, bankfull_flows_file=BANKFULL_FLOWS_FILE, n_workers=1 -# ).run() -# results = SrcSubdiv( -# aoi_dir=AOI_DIR, -# vmann_table=VMANN_INPUT_FILE, -# n_workers=JOB_BRANCH_LIMIT, -# default_channel_n=0.06, # used when feature_id missing from vmann table -# default_overbank_n=0.12, -# ).run() -# assert results - - -# @_skip_no_aoi -# def test_step_nonmonotonic(): -# """Force monotonic in-channel rating curves.""" -# results = SrcNonmonotonic( -# aoi_dir=AOI_DIR, stream_order_min=4, include_branch_zero=True -# ).run() -# assert results - - -# @_skip_no_aoi -# @_skip_no_usgs -# def test_step_usgs(): -# """Calibrate SRCs against USGS rating curves at NWM recurrence flows. -# Needs usgs_elev_table.csv at the AOI root; self-skips when inputs are absent.""" -# results = UsgsRatingCalibrator( -# aoi_dir=AOI_DIR, -# usgs_rating_curve_csv=USGS_RATING_CURVE_CSV, -# usgs_acceptable_gages=USGS_ACCEPTABLE_GAGES, -# nwm_recur_file=NWM_RECUR_FILE, -# n_workers=JOB_BRANCH_LIMIT, -# ).run() -# assert results is not None - - -# @_skip_no_aoi -# @_skip_no_spatial -# def test_step_spatial(): -# """Calibrate SRCs against benchmark inundation points. Samples HAND/HydroID -# rasters at each point; self-skips when the points file is absent.""" -# results = SpatialObsCalibrator( -# aoi_dir=AOI_DIR, -# calib_points_file=CALIB_POINTS_FILE, -# n_workers=JOB_BRANCH_LIMIT, -# ).run() -# assert results is not None - - -# @_skip_no_aoi -# @_skip_no_manual -# def test_step_manual(): -# """Apply a per-feature_id coefficient table to each branch hydroTable. -# Needs MAN_CALB_FILE (aoi_id, feature_id, calb_coef_manual); no-op when the -# AOI has no matching entry.""" -# ManualCalibrator(aoi_dir=AOI_DIR, calibration_file=MAN_CALB_FILE).run() - - -# @_skip_no_aoi -# def test_step_aggregate_post(): -# """Post-calibration aggregation: htable + bridge + road -> AOI root.""" -# BranchAggregator(aoi_dir=AOI_DIR, htable=True, bridge=True, road=True).run() - - -# @_skip_no_aoi -# def test_step_log_scan(): -# """Scan logs/ for error / warning lines into per-AOI summary files.""" -# out = LogScanner(aoi_dir=AOI_DIR, calibration_rerun=False).run() -# assert set(out) == {"errors", "warnings"} diff --git a/tests/test_downloaddata.py b/tests/test_downloaddata.py deleted file mode 100644 index b42987e..0000000 --- a/tests/test_downloaddata.py +++ /dev/null @@ -1,203 +0,0 @@ -# Example Usage: -from pathlib import Path - -import fimbox - -PKG_ROOT = Path(__file__).resolve().parents[1] -REPO_ROOT = Path(__file__).resolve().parents[2] - -test_boundary = PKG_ROOT / "docs" / "test_boundary" / "test_smallB.shp" -OUT_DIR = REPO_ROOT / "out" - -# # Testing the entire NHDPlus data extraction process along with National Flood Hazard Layer data extraction -# # This is OLDER VERSION using EPA AWS S3 Bucket which will get for whole HUC6 region--> not very effective -# def test_getNHDdata(): -# nhd_data = fimbox.getNHDPlusData( -# NHDglobalBoundary = NHDboundary, #Contains all NHDPlus VPU/RPU boundaries -# # inputs_dir = None, #Directory to save input data, if None, direct folder directory will be created -# boundary_path = test_boundary, #Path to the boundary shapefile for which NHDPlus data is to be extracted, OR HUC8 ID -# # huc8: Optional[str] = None, -# # epsg: Optional[int] = None, -# # out_dir: Optional[str] = None, #Directory to save output data, if None, direct folder directory will be created -# # auto_run= True -# ) -# # nhd_data.process_flowlines() -# print(f"Process successful!") - -# def test_get_nfhl(): -# fimbox.DownloadFEMANFHL( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# out_name="fema_nfhl.gpkg", -# # log_path=None, -# ) - -# def test_download_nld(): -# fimbox.DownloadNLD( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# lines_name="NLD_Lines.gpkg", # default; override as needed -# polys_name="NLD_Polygons.gpkg", # default; override as needed -# ) - -##This is for the medium range -# def test_get_nhddata(): -# fimbox.NWMFlowlinesDownloader().download( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# out_name="nwm_subset_streams.gpkg", -# out_layer="flowlines", -# ) - -##Medium range -# def test_get_catchments(): -# fimbox.NWMCatchmentsDownloader().download( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# out_name="nwm_subset_catchments.gpkg", -# out_layer="catchments", -# ) - -# def test_get_lakes(): -# fimbox.NWMLakesDownloader().download( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# out_name="nwm_subset_lakes.gpkg", -# out_layer="lakes", -# ) - - -# # Get all NHD Plus Data -# def test_get_nhd_all(): -# fimbox.getNHDPlusData( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# download_flowlines=True, -# download_catchments=True, -# download_lakes=True, -# resolution="medium", # "high" -> NHDPlus HR flowlines/catchments via pynhd; "medium" (default) -> NWM. Lakes always NWM. -# identifier="nwmmr", # filename prefix; default "nwm" -> nwm_subset_streams.gpkg etc. -# ) - - -# High-resolution flowlines + catchments only (NHDPlus HR via pynhd). -# def test_get_nhd_hr(): -# fimbox.getNHDPlusHRData( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# download_flowlines=True, -# download_catchments=True, -# identifier="nwm", # prefix for saved files -# ) - - -# Bring-your-own flowlines/catchments: map your column names to the canonical -# schema (streams: ID, order_, levpa_id, feature_id[=ID]; catchments: ID). -# def test_normalize_byo_flowlines_catchments(): -# fl = fimbox.normalize_flowlines( -# "path/to/my_flowlines.gpkg", -# field_map={"ID": "nhdplusid", "order_": "streamorde", "levpa_id": "levelpathi"}, -# ) -# cat = fimbox.normalize_catchments( -# "path/to/my_catchments.gpkg", field_map={"ID": "nhdplusid"} -# ) -# assert {"ID", "order_", "levpa_id", "feature_id"}.issubset(fl.columns) -# assert "ID" in cat.columns - - -# Download + process a 3DEP DEM. Reads only the AOI window straight from the -# Planetary Computer COGs over HTTP (no national VRT parse), one snapped -# reprojection, then hole-fill + clip. ~8x faster than the old py3dep path. -# Resolutions: 10 / 30 m seamless nationwide; 1 / 3 m from project lidar where -# it covers the AOI; 60 m is Alaska-only. A resolution with no data for the AOI -# logs + raises DEMResolutionUnavailable (default stays 10 m). -def test_get_dem(): - fimbox.DEMProcessor( - boundary=test_boundary, - output_dir=OUT_DIR, - resolution=10, # 1, 3, 10 (default), 30, 60 - out_name="dem.tif", # default is 3dep_dem_m.tif - # epsg=None, # output CRS; None -> auto UTM zone - # layer=None, # layer name if boundary has multiple - # use_dask=True, # dask chunking for reproject/heal - # chunksize=None, # None -> auto from CPU count; or set px - ) - - -# "give me 1 m, else just 10 m": fallback_to_10m downgrades to 10 m (and logs) -# when the requested resolution isn't available for the AOI. -# def test_get_dem_fallback(): -# fimbox.DEMProcessor( -# boundary=test_boundary, -# output_dir=OUT_DIR, -# resolution=1, # 1 m where lidar exists, else fall back -# out_name="dem.tif", -# fallback_to_10m=True, # default False -> raises if unavailable -# ) - -# Bring your own DEM: pass dem_file and it gets the SAME conditioning as a -# downloaded one (reproject -> hole-fill -> clip to boundary). -# def test_process_byo_dem(): -# fimbox.DEMProcessor( -# boundary=test_boundary, -# output_dir=OUT_DIR, -# out_name="dem.tif", -# resolution=10, -# dem_file="path/to/my_dem.tif", -# ) - -# def test_get_osm_roads(): -# fimbox.DownloadOSMRoads().download( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# out_name="osm_roads.gpkg", -# out_layer="osm_roads", -# ) - -# def test_get_osm_bridges(): -# fimbox.DownloadOSMBridges().download( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# out_name="osm_bridges.gpkg", -# out_layer="osm_bridges", -# ) - - -# # Get the HUC8 information -# def test_get_huc8_info(): -# huc8_info = fimbox.getHUC8Info( -# boundary=test_boundary, -# calc_overlap=True, -# save=True, -# out_dir=OUT_DIR, -# ) -# logging.getLogger(__name__).info(f"HUC8 info:\n{huc8_info}") - - -# USGS gauge points — downloads from the ArcGIS Online FeatureServer. -# Uncomment to run live; the smoke test below always runs. - -# def test_download_usgs_gages(): -# """Download USGS gauges inside the test boundary into ../out/usgs_gages.gpkg.""" -# gdf = fimbox.DownloadUSGSGages().download( -# boundary=test_boundary, -# aoi_id="08060202", -# out_dir=OUT_DIR, -# out_name="usgs_gages.gpkg", -# out_layer="usgs_gages", -# ) -# log = logging.getLogger(__name__) -# log.info(f"USGS gauges downloaded: {len(gdf)} features") -# assert {"location_id", "feature_id", "aoi_id", "source"}.issubset(gdf.columns) - - -# def test_usgs_gages_signature(): -# """Smoke test: DownloadUSGSGages is exported and has the documented API.""" -# import inspect - -# assert hasattr(fimbox, "DownloadUSGSGages") -# sig = inspect.signature(fimbox.DownloadUSGSGages.download) -# expected = {"boundary", "aoi_id", "out_dir", "out_name", "out_layer"} -# assert expected.issubset(sig.parameters.keys()), ( -# f"DownloadUSGSGages.download missing kwargs: {expected - set(sig.parameters)}" -# ) diff --git a/tests/test_fimevaluation.py b/tests/test_fimevaluation.py deleted file mode 100644 index 7f899ba..0000000 --- a/tests/test_fimevaluation.py +++ /dev/null @@ -1,123 +0,0 @@ -""" -Author: Supath Dhital -Date Created: July 2026 - -FIM evaluation tests: query benchmark FIMs from the FIMbench database -(queryBenchmarkFIM) and evaluate candidate FIMs against them with FIMeval -(evaluateFIM). - -Point AOI_DIR at a working directory that already has flood maps in -fim-outputs/ (produced by test_fimgeneration). fimbench and fimeval install -together with fimbox. -""" - -from __future__ import annotations - -from pathlib import Path - -AOI_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" - -# Boundary-extraction method: "smallest_extent" | "convex_hull" | "AOI" -METHOD_NAME = "smallest_extent" - -# Optional query filters (edit to match your AOI / event). -EVENT_DATE = None # e.g. "2017-08-30" -START, END = "2016-04-01", "2026-01-01" -TIER = None # e.g. "HWM", "tier1" - - -# COMBINED — the whole evaluation pipeline in one go: query the FIMbench -# catalog with the AOI's newest flood extent, download the matched benchmark -# assets into /benchmark-data/, then run FIMeval (EvaluateFIM + -# contingency maps + metric plots) on the staged case. -def test_fimevaluation_combined(): - from fimbox import evaluateFIM, queryBenchmarkFIM - - query = queryBenchmarkFIM( - AOI_DIR, # footprint = newest extent raster in /fim-outputs/ - # raster_path="my_fim.tif", # explicit candidate raster instead - # boundary_path="my_aoi.gpkg", # or an AOI boundary vector - # huc8="03020201", # narrow by basin - tier=TIER, # narrow by benchmark tier; None -> all tiers - event_date=EVENT_DATE, # exact event date; None -> ignore - start_date=START, # date-range start - end_date=END, # date-range end - area=True, # add overlap % / km^2 per match - download=True, # fetch GeoTIFF + GeoPackage - # out_dir="downloads/", # default /benchmark-data/ - ) - print(query) # pretty match summary - - result = evaluateFIM( - AOI_DIR, - # candidate=None, # default: every extent .tif in fim-outputs/ - # benchmark=None, # default: newest .tif in benchmark-data/ - # case_name=None, # default: first candidate's stem - method_name=METHOD_NAME, - # aoi_boundary="my_aoi.gpkg", # required when method_name="AOI" - # pwb_dir="my_pwb.gpkg", # own permanent-water-bodies vector - # target_crs="EPSG:32633", # outside CONUS (default EPSG:5070) - # target_resolution=10, # m, when resolutions differ - contingency_map=True, - plot_metrics=True, - # building_footprint=True, # building-level agreement (GEE auth) - ) - assert result.case_dir.is_dir() - assert result.output_dir.is_dir() - assert "benchmark" in result.benchmark.name.lower() - assert result.metrics_files, "FIMeval produced no metrics CSV" - print(result.metrics) - - -# # STEP BY STEP — each stage on its own. - -# def test_step_query_catalog_only(): -# """Catalog-only discovery: what benchmarks exist in a date range -# (no AOI, no download). Returns plain dicts straight from the catalog.""" -# from fimbox import queryBenchmarkFIM - -# response = queryBenchmarkFIM(start_date=START, end_date=END) -# print(response) -# for record in response.records: -# print(record) - - -# def test_step_query_by_filename(): -# """Direct download by exact catalog filename.""" -# from fimbox import queryBenchmarkFIM - -# response = queryBenchmarkFIM( -# file_name="HWM_10_0m_20160928_20161009_780051W352232N_BM.tif", -# download=True, -# out_dir=AOI_DIR / "benchmark-data", -# ) -# assert response.downloads - - -# def test_step_evaluate_own_benchmark(): -# """Evaluate against a benchmark raster you already have on disk -# (skips the FIMbench query entirely).""" -# from fimbox import evaluateFIM - -# result = evaluateFIM( -# AOI_DIR, -# benchmark="path/to/my_benchmark.tif", -# case_name="own_benchmark", -# method_name="convex_hull", -# ) -# assert result.metrics_files - - -# def test_step_building_footprint(): -# """Building-level agreement analysis. Uses the Microsoft global building -# footprints via Google Earth Engine by default (pops a GEE auth prompt), -# or pass building_footprint_file= to use your own vector.""" -# from fimbox import evaluateFIM - -# result = evaluateFIM( -# AOI_DIR, -# method_name=METHOD_NAME, -# building_footprint=True, -# # building_footprint_file="path/to/footprints.gpkg", -# ) -# assert result.output_dir.is_dir() diff --git a/tests/test_fimgeneration.py b/tests/test_fimgeneration.py deleted file mode 100644 index c13d1a2..0000000 --- a/tests/test_fimgeneration.py +++ /dev/null @@ -1,102 +0,0 @@ -""" -Author: Supath Dhital -Date Updated: June 2026 - -FIM generation driven by the CSVs in /discharge-inputs/. - -generateFIM(aoi_dir).from_discharge_inputs(...) selects which discharge CSVs to -run and generates an inundation extent raster for each (named after the input -CSV) in /fim-outputs/. Pass depth=True to also write a depth raster. - -Selection modes: - * nothing -> every CSV in discharge-inputs/ - * csv= -> that one CSV - * date="YYYY-MM-DD" -> CSVs whose filename carries that day/instant stamp - * start=.., end=.. -> CSVs whose YYYYMMDD token falls in the range - -Point AOI_DIR at a working directory that already has branches and at least one -discharge CSV (produced by the streamflow pipeline, e.g. getNWMretrospective). -""" - -from __future__ import annotations - -from pathlib import Path - -import pytest - -from fimbox import generateFIM - -AOI_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" -N_WORKERS = 4 - -# Optional selection filters (edit to match the CSVs you have). -EVENT_DATE = "2020-05-20 12:00:00" -START = "2020-05-19" -END = "2020-05-22" - -_BRANCHES_DIR = ( - AOI_DIR / "watershed-data" / "branches" - if (AOI_DIR / "watershed-data" / "branches").is_dir() - else AOI_DIR / "branches" -) -_skip_no_branches = pytest.mark.skipif( - not _BRANCHES_DIR.is_dir(), reason=f"AOI not present: {_BRANCHES_DIR}" -) - - -# @_skip_no_branches -# def test_extract_feature_ids(): -# out_csv = extract_feature_ids(AOI_DIR) -# assert out_csv.is_file() -# print(f"\nfeature_id.csv -> {out_csv}") - - -# default: generate FIM for every discharge CSV in the AOI -@_skip_no_branches -def test_fim_all_discharge_inputs(): - results = generateFIM( - AOI_DIR, n_workers=N_WORKERS, depth=True - ).from_discharge_inputs() - assert results - for r in results: - print(f" extent={r.extent_path}") - assert r.extent_path is not None and Path(r.extent_path).is_file() - - -# # a specific CSV -# @_skip_no_branches -# def test_fim_specific_csv(): -# csvs = sorted((AOI_DIR / "discharge-inputs").glob("*.csv")) -# if not csvs: -# pytest.skip("no discharge CSVs to pick from") -# results = generateFIM(AOI_DIR, n_workers=N_WORKERS).from_discharge_inputs(csv=csvs[0]) -# assert len(results) == 1 - - -# # match by date stamp in the filename -# @_skip_no_branches -# def test_fim_by_date(): -# results = generateFIM(AOI_DIR, n_workers=N_WORKERS).from_discharge_inputs( -# date=EVENT_DATE -# ) -# assert results - - -# # match by date range -# @_skip_no_branches -# def test_fim_by_range(): -# results = generateFIM(AOI_DIR, n_workers=N_WORKERS).from_discharge_inputs( -# start=START, end=END -# ) -# assert results - - -# # also write the depth raster -# @_skip_no_branches -# def test_fim_with_depth(): -# results = generateFIM( -# AOI_DIR, n_workers=N_WORKERS, depth=True -# ).from_discharge_inputs(date=EVENT_DATE) -# assert results -# r = results[0] -# assert r.depth_path is not None and Path(r.depth_path).is_file() diff --git a/tests/test_generate_dem_diff.py b/tests/test_generate_dem_diff.py deleted file mode 100644 index 6f30194..0000000 --- a/tests/test_generate_dem_diff.py +++ /dev/null @@ -1,71 +0,0 @@ -""" -Tests for bridge DEM processing pipeline. -Step 1 (generateBridgeRaster): streams USGS LiDAR and writes per-bridge .tif -Step 2 (BridgeDEMDiff): computes lidar_elev - dem_elev and saves bridge_elev_diff.tif -""" - -import logging -from pathlib import Path - -import fimbox - -log = logging.getLogger(__name__) - -OUT_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" / "watershed-data" -bridge_gpkg = OUT_DIR / "osm_bridges_subset.gpkg" -dem_path = OUT_DIR / "dem.tif" -out_dir = OUT_DIR - - -# check which bridges already have rasters vs still pending (safe to run anytime) -def test_bridge_raster_status(): - info = fimbox.generateBridgeRaster( - bridge_gpkg=bridge_gpkg, - out_dir=out_dir, - ).status() - assert "total" in info - - -# download LiDAR and build per-bridge elevation tifs -def test_generate_bridge_raster(): - tif_dir = fimbox.generateBridgeRaster( - bridge_gpkg=bridge_gpkg, - out_dir=out_dir, - resolution=10.0, - buffer_m=10.0, - n_workers=4, - # id_col="my_id", # only needed if gpkg has no 'osmid' column - ).run() - log.info(f"Per-bridge tifs --> {tif_dir}") - - -# compute difference raster -def test_bridge_dem_diff(): - out_path = fimbox.BridgeDEMDiff( - dem_path=dem_path, - lidar_tif_dir=OUT_DIR / "bridge_dem" / "lidar_osm_rasters", - bridge_gpkg=bridge_gpkg, - out_dir=out_dir, - out_name="bridge_elev_diff.tif", - n_workers=4, - ).run() - log.info(f"Bridge diff raster --> {out_path}") - - -# Run both steps end-to-end -# def test_full_pipeline(): -# tif_dir = fimbox.generateBridgeRaster( -# bridge_gpkg=bridge_gpkg, -# out_dir=out_dir, -# resolution=10.0, -# n_workers=4, -# ).run() -# -# out_path = fimbox.BridgeDEMDiff( -# dem_path=dem_path, -# lidar_tif_dir=tif_dir, -# bridge_gpkg=bridge_gpkg, -# out_dir=out_dir, -# n_workers=4, -# ).run() -# print(f"Done: {out_path}") diff --git a/tests/test_getallinputdata.py b/tests/test_getallinputdata.py deleted file mode 100644 index 4989dba..0000000 --- a/tests/test_getallinputdata.py +++ /dev/null @@ -1,94 +0,0 @@ -# Example Usage: -from pathlib import Path - -import fimbox - -PKG_ROOT = Path(__file__).resolve().parents[1] -REPO_ROOT = Path(__file__).resolve().parents[2] - -test_boundary = PKG_ROOT / "docs" / "test_boundary" / "test_smallB.shp" -OUT_DIR = REPO_ROOT / "out" -test_huc8 = "08060202" # Yazoo River basin, MS - - -# Combined preprocessing pipeline tests -# Run full pipeline from a boundary shapefile -def test_preprocess_all_from_boundary(): - pp = fimbox.getAllInputData( - boundary=test_boundary, - out_dir=OUT_DIR, - buffer_m=5000, # metres to buffer boundary for data downloads - headwater_buffer_cells=8, # pixels to shrink buffer for headwater clip - get_flowlines=True, # set False to use your own flowlines and corresponding catchments - get_catchments=True, # set False to skip NWM catchments--> use - resolution="medium", # "high" -> NHDPlus HR flowlines/catchments via pynhd; "medium" (default) -> NWM. Lakes always NWM. - identifier="nwmmr", # filename prefix for ALL source files; flows download->processing. Default "nwm". - ) - pp.run() - - -# Bring your own flowlines / catchments / DEM (any column names, any source). -# Pass the file paths + field maps; flowlines/catchments are normalised to the -# pipeline schema (streams: ID, order_, levpa_id, feature_id[=ID]; catchments: ID), -# the DEM is reprojected/clipped/hole-filled, and all files are saved under the -# chosen identifier prefix so the whole pipeline picks them up automatically. -# def test_preprocess_byo_inputs(): -# pp = fimbox.getAllInputData( -# boundary=test_boundary, -# out_dir=OUT_DIR, -# flowlines="path/to/my_flowlines.gpkg", -# catchments="path/to/my_catchments.gpkg", -# stream_fields={"ID": "nhdplusid", "order_": "streamorde", "levpa_id": "levelpathi"}, -# catchment_fields={"ID": "nhdplusid"}, # must match the flowline reach id -# dem="path/to/my_dem.tif", # reprojected, clipped, and hole-filled like a downloaded DEM -# identifier="3dhp", # files saved as 3dhp_subset_streams.gpkg etc.; whole pipeline follows it -# ) -# pp.run() - - -# # Run full pipeline from a HUC8 ID -# # get_flowlines / get_catchments default to True (downloads everything, -# # including OSM bridges). Set either to False to skip that dataset and use -# # your own instead. -# def test_preprocess_all_from_huc8(): -# pp = fimbox.getAllInputData( -# huc8=test_huc8, -# out_dir=OUT_DIR, -# buffer_m=2000, -# headwater_buffer_cells=8, -# get_flowlines=True, # set False to use your own flowlines and corresponding catchments -# get_catchments=True, # set False to skip NWM catchments--> use your own in later steps -# ) -# pp.run() - - -# Same pipeline, but bring your own flowlines/catchments -# (skips the NWM flowline + catchment downloads; everything else still runs) -# def test_preprocess_all_byo_flowlines_catchments(): -# pp = fimbox.getAllInputData( -# huc8=test_huc8, -# out_dir=OUT_DIR, -# buffer_m=2000, -# headwater_buffer_cells=8, -# get_flowlines=False, -# get_catchments=False, -# ) -# pp.run() - - -# Run individual steps -# def test_preprocess_dem_only(): -# pp = fimbox.getAllInputData(boundary=test_boundary, out_dir=OUT_DIR) -# pp.run_dem() - -# def test_preprocess_nhd_only(): -# pp = fimbox.getAllInputData(boundary=test_boundary, out_dir=OUT_DIR) -# pp.run_nhd() - -# def test_preprocess_nld_only(): -# pp = fimbox.getAllInputData(boundary=test_boundary, out_dir=OUT_DIR) -# pp.run_nld() - -# def test_preprocess_osm_only(): -# pp = fimbox.getAllInputData(boundary=test_boundary, out_dir=OUT_DIR) -# pp.run_osm() diff --git a/tests/test_nwmstreamflow.py b/tests/test_nwmstreamflow.py deleted file mode 100644 index 4169e18..0000000 --- a/tests/test_nwmstreamflow.py +++ /dev/null @@ -1,87 +0,0 @@ -""" -Author: Supath Dhital -Date Created: June 2026 - -Streamflow retrieval / plot / statistics — minimal, call-the-function tests. - -Point AOI_DIR at a working directory whose feature_id.csv exists (or pass a -feature_ids list / CSV per call). Edit the dates and USGS site to your basin, -then run the functions you want. -""" - -from __future__ import annotations - -from pathlib import Path - -from fimbox import ( - getNWMretrospective, -) - -AOI_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" - -START = "2016-10-05" -END = "2016-10-20" -EVENT = "2020-10-10 21:00:00" - - -# retrospective — different extraction combinations -def test_retrospective_event_date(): - getNWMretrospective(AOI_DIR, date=EVENT) - - -# def test_retrospective_range_continuous(): -# # start + end, nothing else -> one CSV per hour -# getNWMretrospective(AOI_DIR, start=START, end=END) - - -# def test_retrospective_range_sortby(): -# # start + end + sortby -> one aggregated CSV -# getNWMretrospective(AOI_DIR, start=START, end=END, sortby="maximum") - - -# def test_retrospective_feature_ids_list(): -# # pass feature_ids directly instead of relying on the AOI's feature_id.csv -# getNWMretrospective(AOI_DIR, feature_ids=[FEATURE_ID], date=EVENT) - - -# # forecast — different combinations -# def test_forecast_shortrange(): -# getNWMforecast(AOI_DIR, "shortrange") - - -# def test_forecast_mediumrange_maxsort(): -# getNWMforecast(AOI_DIR, "mediumrange", sort_by="maximum") - - -# def test_forecast_specific_cycle(): -# getNWMforecast(AOI_DIR, "shortrange", forecast_date="2024-06-01", hour=12) - - -# # USGS observations -# def test_usgs_fetch(): -# USGSData(AOI_DIR).fetch([USGS_SITE], START, END) - - -# def test_usgs_feature_id_pairs(): -# # which USGS gage falls on which reach (feature_id) within the AOI -# pairs = get_usgs_fid_pairs(AOI_DIR) -# print(pairs) - - -# # plots -# def test_plot_feature_id(): -# plot_nwm(AOI_DIR, [FEATURE_ID], START, END) - - -# def test_plot_usgs(): -# plot_usgs(AOI_DIR, [USGS_SITE], START, END) - - -# def test_plot_usgs_and_feature_id(): -# # time series overlay of USGS and the NWM feature_id together -# plot_comparison(AOI_DIR, FEATURE_ID, USGS_SITE, START, END) - - -# # statistics -# def test_statistics_usgs_vs_nwm(): -# calculate_statistics(AOI_DIR, FEATURE_ID, USGS_SITE, START, END) diff --git a/tests/test_preprocessDEM.py b/tests/test_preprocessDEM.py deleted file mode 100644 index 3e7c636..0000000 --- a/tests/test_preprocessDEM.py +++ /dev/null @@ -1,28 +0,0 @@ -# Example Usage: -import logging -from pathlib import Path - -import fimbox - -log = logging.getLogger(__name__) - -PKG_ROOT = Path(__file__).resolve().parents[1] -REPO_ROOT = Path(__file__).resolve().parents[2] - -boundary = PKG_ROOT / "docs" / "test_boundary" / "test_smallB.shp" -OUT_DIR = REPO_ROOT / "out" - - -def test_process_dem(): - output_path = fimbox.DEMProcessor( - boundary=boundary, - resolution=10, # 3DEP resolution in m: 1, 3, 10 (default), 30, 60 - output_dir=OUT_DIR / "dem_test", - # layer=None, # if boundary is a geopackage with multiple layers - # dem_file=None, # local DEM to condition instead of fetching - # epsg=None, # output CRS EPSG; None auto-detects the UTM zone - # fallback_to_10m=False, # if resolution unavailable, use 10m not raise - # use_dask=True, # dask chunking for the reproject/heal stage - # chunksize=None, # dask chunk edge in px; None -> auto from CPU count - ).result_path - log.info(f"3DEP DEM --> {output_path}") diff --git a/tests/test_preprocessing_hucs.py b/tests/test_preprocessing_hucs.py deleted file mode 100644 index 31f9aea..0000000 --- a/tests/test_preprocessing_hucs.py +++ /dev/null @@ -1,19 +0,0 @@ -# importing the fimbox preprocessing module to test HUCChecker -import logging - -import fimbox - -log = logging.getLogger(__name__) -checker = fimbox.HUCChecker() - - -def test_huc_checker(): - # Single HUC Query - r = checker.check_any("03020202", strict=False) - log.info(f"total={r.n_total} found={r.n_found} missing={r.n_missing}") - log.info(f"missing: {r.missing_hucs}") - - # List of HUCs Query - r = checker.check_any(["01010001", "99999999"], strict=False) - log.info(f"total={r.n_total} found={r.n_found} missing={r.n_missing}") - log.info(f"missing: {r.missing_hucs}") From 7d2e645e593981d59ab1783d4d31f39292448e7a Mon Sep 17 00:00:00 2001 From: Manjila Singh Date: Tue, 28 Jul 2026 08:52:09 -0500 Subject: [PATCH 4/4] restore tests folder --- tests/README.md | 36 ++ tests/conftest.py | 5 + tests/test_branchprocessing.py | 1008 ++++++++++++++++++++++++++++++ tests/test_calibrate_pipeline.py | 289 +++++++++ tests/test_downloaddata.py | 203 ++++++ tests/test_fimevaluation.py | 123 ++++ tests/test_fimgeneration.py | 102 +++ tests/test_generate_dem_diff.py | 71 +++ tests/test_getallinputdata.py | 94 +++ tests/test_nwmstreamflow.py | 87 +++ tests/test_preprocessDEM.py | 28 + tests/test_preprocessing_hucs.py | 19 + 12 files changed, 2065 insertions(+) create mode 100644 tests/README.md create mode 100644 tests/conftest.py create mode 100644 tests/test_branchprocessing.py create mode 100644 tests/test_calibrate_pipeline.py create mode 100644 tests/test_downloaddata.py create mode 100644 tests/test_fimevaluation.py create mode 100644 tests/test_fimgeneration.py create mode 100644 tests/test_generate_dem_diff.py create mode 100644 tests/test_getallinputdata.py create mode 100644 tests/test_nwmstreamflow.py create mode 100644 tests/test_preprocessDEM.py create mode 100644 tests/test_preprocessing_hucs.py diff --git a/tests/README.md b/tests/README.md new file mode 100644 index 0000000..e9c5433 --- /dev/null +++ b/tests/README.md @@ -0,0 +1,36 @@ +### Tests +
+ +The test suite doubles as the usage reference for `fimbox`: every stage of the workflow has a test file whose top-level constants (AOI paths, input files, worker counts) are meant to be edited to point at your own data. Tests skip cleanly when the referenced AOI or input file is absent. Most files also keep commented-out variants showing alternative call patterns and optional parameters. + +**Folder contents** + +| File | What it covers | +|---|---| +| `conftest.py` | Suite-wide logging setup (`configure_cli_logging`). | +| `test_preprocessing_hucs.py` | HUC validation with `HUCChecker` (single HUC, lists, .lst/.csv files, strict mode). | +| `test_downloaddata.py` | Individual dataset downloaders: DEM, NHDPlus/NWM hydrography, FEMA NFHL, NLD levees, OSM roads/bridges, USGS gages. | +| `test_getallinputdata.py` | The combined `getAllInputData` pipeline from a boundary shapefile, including bring-your-own flowlines/catchments/DEM. | +| `test_preprocessDEM.py` | `DEMProcessor` fetching and conditioning (resolutions, local DEM, CRS handling). | +| `test_generate_dem_diff.py` | Bridge LiDAR rasters (`generateBridgeRaster`, with `status()` check) and `BridgeDEMDiff` mosaicking. | +| `test_branchprocessing.py` | `BranchDerivation`, `AOIProcessingConfig`, `calculate_allbranches`, plus step-level tests for the BranchZero and CreateHAND substeps. | +| `test_calibrate_pipeline.py` | The full SRC calibration pipeline via one `run_calibration()` call with every `CalibrationConfig` parameter spelled out, plus one test per calibration stage. | +| `test_nwmstreamflow.py` | Streamflow retrieval (`getNWMretrospective`, `getNWMforecast`, `USGSData`), plotting, and KGE/NSE/PBias statistics. | +| `test_fimgeneration.py` | FIM generation from `discharge-inputs/` CSVs with date/range selection and depth output options. | +| `test_fimevaluation.py` | Benchmark FIM query/download via `queryBenchmarkFIM` (FIMbench) and candidate-vs-benchmark evaluation via `evaluateFIM` (FIMeval). | + +### Running +
+ +```bash +# from the repo root, with the environment activated +pytest tests/ -v # whole suite +pytest tests/test_calibrate_pipeline.py -v # one stage +pytest tests/test_branchprocessing.py -v -k hand # one test by keyword +``` + +Before running, edit the constants at the top of each test file (for example `AOI_DIR`, `BANKFULL_FLOWS_FILE`, `N_WORKERS`) to match your machine. The calibration lookup tables referenced by the tests ship in the repo [`data/`](../data/) folder. + +The expected order when building an AOI from scratch mirrors the workflow: `test_getallinputdata` (stage inputs), `test_generate_dem_diff` (optional bridge healing), `test_branchprocessing` (HAND + SRC), `test_calibrate_pipeline` (calibration), `test_nwmstreamflow` (discharge), `test_fimgeneration` (flood maps), `test_fimevaluation` (benchmark evaluation). + +**For more usage notes refer to the [docs](../docs/) for the `fimbox` python package.** diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..6d0e03c --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,5 @@ +"""Test-suite-wide logging setup.""" + +from fimbox.logging_utils import configure_cli_logging + +configure_cli_logging() diff --git a/tests/test_branchprocessing.py b/tests/test_branchprocessing.py new file mode 100644 index 0000000..981aa40 --- /dev/null +++ b/tests/test_branchprocessing.py @@ -0,0 +1,1008 @@ +""" +Author: Supath Dhital +Date Updated: May 2026 + + +Branch processing tests. + +Run order: + 1. test_branch_derivation — level paths, branch polygons, branch list + 2. test_branch_zero_full — DEM clip, AGREE, pit-fill, D8 flowdir + 3. test_create_hand — full HAND generation (flow accum → split reaches) +""" + +import logging +from pathlib import Path + +# single steps IMPORTS +from fimbox import ( + BranchDerivation, +) + +log = logging.getLogger(__name__) + +# imports used only by the B-series CreateHAND step tests below. +# BranchZero substeps (StreamBooleanRasterizer, HydroenforceDEM, FlowdirDEM, +# HeadwaterRasterizer, LevelPathBooleanRasterizer, rasterize_3d_levee_lines, +# burn_levee_elevations) are exercised indirectly via test_step_Z1, so they +# are not imported here. + +# AOI parameters — point this at any user-supplied AOI working directory. +OUT_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" / "watershed-data" + +# Source-data filename prefix. +IDENTIFIER = "nwmmr" + +# Tunable CreateHAND parameters- All have sensible defaults in CreateHAND itself. +PARAMS_CREATE_HAND = dict( + cost_distance_tolerance=50.0, # m, lateral cost distance + lateral_elevation_threshold=10, # m, lateral thalweg drop cap + max_split_distance_m=1500.0, # m, split-reach max length + slope_min=0.0001, # rise/run floor + lakes_buffer_dist_m=100.0, # m, lake-boundary buffer + # SRC / crosswalk + mannings_n=0.06, # channel roughness + stage_min_m=0.0, # SRC stage ladder start + stage_interval_m=0.3048, # SRC stage step (1 ft) + stage_max_m=25.2984, # SRC stage ladder end (~83 ft) + min_catchment_area=0.25, # km^2, short-reach replace threshold + min_stream_length=0.5, # km, short-reach replace threshold + crosswalk_max_distance_m=100.0, # m, midpoint-to-NWM-flowline cap + # SRC slope source feeding Manning's equation: + # "iris_sword" (default) - IRIS-SWORD slope on order>=4 streams, else DEM + # "dem" - DEM rise/run slope only + # "hfab" - hydrofabric native slope, else DEM fallback + src_slope_source="iris_sword", + # IRIS-SWORD slope table (feature_id, slope_iris_sword). None -> the table + # shipped in fimbox/data is used when src_slope_source == "iris_sword". + iris_slope_csv=None, + # Hydrofabric slope column when it isn't the usual 'Slope'/'So'. + hfab_slope_column=None, +) + +DEM = OUT_DIR / "dem.tif" +STREAMS = OUT_DIR / f"{IDENTIFIER}_subset_streams.gpkg" +BOUNDARY_BUF = OUT_DIR / "wbd_buffered.gpkg" +CATCHMENTS = OUT_DIR / f"{IDENTIFIER}_catchments_proj_subset.gpkg" +HEADWATERS = ( + OUT_DIR / f"{IDENTIFIER}_headwater_points_subset.gpkg" + if (OUT_DIR / f"{IDENTIFIER}_headwater_points_subset.gpkg").is_file() + else OUT_DIR / f"{IDENTIFIER}_headwaters.gpkg" +) +LEVELPATH_EXT = OUT_DIR / f"{IDENTIFIER}_subset_streams_levelPaths_extended.gpkg" +BRIDGE_DIFF = OUT_DIR / "bridge_elev_diff.tif" +NLD_LEVEES = OUT_DIR / "3d_nld_subset_levees_burned.gpkg" + +# optional files +WBD8_CLP = OUT_DIR / "wbd8_clp.gpkg" +LAKES = OUT_DIR / f"{IDENTIFIER}_lakes_proj_subset.gpkg" +LEVEE_AREAS = OUT_DIR / "LeveeProtectedAreas_subset.gpkg" +LEVEE_LP_CSV = OUT_DIR / "levee_levelpaths.csv" + +# branch-zero derived paths +BRANCH_DIR = OUT_DIR / "branches" / "0" +BRANCH_ID = "0" +DEM_BRANCH = BRANCH_DIR / f"dem_{BRANCH_ID}.tif" +FLOWDIR = BRANCH_DIR / f"flowdir_d8_burned_filled_{BRANCH_ID}.tif" +HW_RASTER = BRANCH_DIR / f"headwaters_{BRANCH_ID}.tif" +STREAM_BOOL = BRANCH_DIR / f"flows_grid_boolean_{BRANCH_ID}.tif" +DEM_BURNED = BRANCH_DIR / f"dem_burned_{BRANCH_ID}.tif" +DEM_FILLED = BRANCH_DIR / f"dem_burned_filled_{BRANCH_ID}.tif" + +# REM + filtered catchment paths +REM = BRANCH_DIR / f"rem_{BRANCH_ID}.tif" +REM_ZEROED = BRANCH_DIR / f"rem_zeroed_masked_{BRANCH_ID}.tif" +CATCH_POLY = BRANCH_DIR / f"gw_catchments_reaches_{BRANCH_ID}.gpkg" +FILT_CATCH = ( + BRANCH_DIR / f"gw_catchments_reaches_filtered_addedAttributes_{BRANCH_ID}.gpkg" +) +FILT_FLOWS = BRANCH_DIR / f"demDerived_reaches_split_filtered_{BRANCH_ID}.gpkg" +FILT_TIF = ( + BRANCH_DIR / f"gw_catchments_reaches_filtered_addedAttributes_{BRANCH_ID}.tif" +) + +# SRC / crosswalk / hydroTable outputs (steps 16-21) +SLOPES_MASKED = BRANCH_DIR / f"slopes_d8_dem_meters_masked_{BRANCH_ID}.tif" +STAGE_TXT = BRANCH_DIR / f"stage_{BRANCH_ID}.txt" +CATCHLIST_TXT = BRANCH_DIR / f"catch_list_{BRANCH_ID}.txt" +SRC_BASE_CSV = BRANCH_DIR / f"src_base_{BRANCH_ID}.csv" +XWALK_CATCH = ( + BRANCH_DIR + / f"gw_catchments_reaches_filtered_addedAttributes_crosswalked_{BRANCH_ID}.gpkg" +) +XWALK_FLOWS = ( + BRANCH_DIR + / f"demDerived_reaches_split_filtered_addedAttributes_crosswalked_{BRANCH_ID}.gpkg" +) +SRC_FULL_CSV = BRANCH_DIR / f"src_full_crosswalked_{BRANCH_ID}.csv" +SRC_JSON = BRANCH_DIR / f"src_{BRANCH_ID}.json" +XWALK_CSV = BRANCH_DIR / f"crosswalk_table_{BRANCH_ID}.csv" +HYDRO_TABLE = BRANCH_DIR / f"hydroTable_{BRANCH_ID}.csv" +ROADS_CSV = BRANCH_DIR / f"osm_roads_fimpact_{BRANCH_ID}.csv" +BRIDGES_GPKG = BRANCH_DIR / f"osm_bridge_centroids_{BRANCH_ID}.gpkg" + +# HAND generation derived paths +FLOWACCUM = BRANCH_DIR / f"flowaccum_d8_burned_filled_{BRANCH_ID}.tif" +STREAM_PIX = BRANCH_DIR / f"demDerived_streamPixels_{BRANCH_ID}.tif" +THALWEG_ADJ = BRANCH_DIR / f"dem_lateral_thalweg_adj_{BRANCH_ID}.tif" +FLOWDIR_STR = BRANCH_DIR / f"flowdir_d8_burned_filled_flows_{BRANCH_ID}.tif" +THALWEG_COND = BRANCH_DIR / f"dem_thalwegCond_{BRANCH_ID}.tif" +SLOPES_D8 = BRANCH_DIR / f"slopes_d8_dem_{BRANCH_ID}.tif" +STREAM_ORDER = BRANCH_DIR / f"streamOrder_{BRANCH_ID}.tif" +SN_CATCH = BRANCH_DIR / f"sn_catchments_reaches_{BRANCH_ID}.tif" +DEM_REACHES = BRANCH_DIR / f"demDerived_reaches_{BRANCH_ID}.gpkg" +SPLIT_REACHES = BRANCH_DIR / f"demDerived_reaches_split_{BRANCH_ID}.gpkg" +SPLIT_PTS = BRANCH_DIR / f"demDerived_reaches_split_points_{BRANCH_ID}.gpkg" +GW_REACHES = BRANCH_DIR / f"gw_catchments_reaches_{BRANCH_ID}.tif" +PIXEL_PTS = BRANCH_DIR / f"flows_points_pixels_{BRANCH_ID}.gpkg" +GW_PIXELS = BRANCH_DIR / f"gw_catchments_pixels_{BRANCH_ID}.tif" + + +# ========================== +# COMBINED — the whole branch pipeline in one go, matching the step-by-step +# sequence exactly: +# Step Z0 BranchDerivation — level paths, branch polygons, branch_ids.lst +# Step Z1 BranchZero — whole-AOI DEM clip + AGREE + pit-fill + D8 +# (serial, in the main process, branch_id="0") +# Step B non-zero branches — BranchZero + CreateHAND per branch, in parallel +# +# Every parameter is spelled out so this test doubles as the parameter reference. +# Optional inputs (bridge_diff, levees, headwaters, levelpaths_extended) are +# resolved from OUT_DIR and passed only when the file exists on disk, matching +# the step-by-step behaviour. +# ============================ +def test_branchprocessing_combined(): + """Run the full branch pipeline in one call. + + Order inside calculate_allbranches: + 1. BranchZero for branch 0 — serial in main process, always first. + Branch 0 gets DEM clip / AGREE / pit-fill / D8 flowdir only. + CreateHAND does NOT run on branch 0; its flowdir is the shared + input every non-zero branch reads. After deny-list cleanup + branches/0/ keeps only what is NOT listed in deny_branch_zero.lst + (no hydroTable — that is expected and correct by design). + 2. All non-zero branches in parallel via Dask — BranchZero then the + full 22-step CreateHAND. Non-zero branches produce hydroTable. + 3. Deny-list cleanup removes intermediates from every branch dir. + """ + from fimbox import AOIProcessingConfig, calculate_allbranches + from fimbox._dask import _resolve_n_workers + + BranchDerivation( + out_dir=OUT_DIR, + branch_id_attribute="levpa_id", + reach_id_attribute="ID", + branch_buffer_distance_meters=7000.0, + ).run() + + bridge_diff = BRIDGE_DIFF if BRIDGE_DIFF.exists() else None + levee_gpkg = NLD_LEVEES if NLD_LEVEES.exists() else None + headwaters = HEADWATERS if HEADWATERS.exists() else None + levelpaths_extended = LEVELPATH_EXT if LEVELPATH_EXT.exists() else None + + n_workers = _resolve_n_workers() + cfg = AOIProcessingConfig( + aoi_dir=OUT_DIR, + branch_list_path=OUT_DIR / "branch_ids.lst", + # BranchZero inputs (whole-AOI, branch_id="0") + dem_path=DEM, + streams_gpkg=STREAMS, + boundary_gpkg=BOUNDARY_BUF, + bridge_elev_diff_path=bridge_diff, + levee_gpkg_path=levee_gpkg, + headwaters_gpkg=headwaters, + levelpaths_extended_gpkg=levelpaths_extended, + # AGREE DEM conditioning + agree_buffer_m=15.0, + agree_smooth_drop=10.0, + agree_sharp_drop=1000.0, + # CreateHAND geometry + cost_distance_tolerance=50.0, + lateral_elevation_threshold=10, + max_split_distance_m=1500.0, + slope_min=0.0001, + lakes_buffer_dist_m=100.0, + # SRC / crosswalk + mannings_n=0.06, + stage_min_m=0.0, + stage_interval_m=0.3048, + stage_max_m=25.2984, + min_catchment_area=0.25, + min_stream_length=0.5, + crosswalk_max_distance_m=100.0, + # SRC slope source: "iris_sword" | "dem" | "hfab" + src_slope_source="iris_sword", + iris_slope_csv=None, + hfab_slope_column=None, + # execution + n_workers=n_workers, + keep_failed_branches=True, + delete_deny_list=True, + ) + + result = calculate_allbranches( + cfg, + run_branch_zero=True, + delete_deny_list=True, + deny_unit_list=Path(__file__).resolve().parent.parent + / "config" + / "deny_unit.lst", + branch_ids_csv=OUT_DIR / "branch_ids.csv", + ) + + # Branch 0 now runs BranchZero + full CreateHAND (same as every non-zero branch). + b0 = OUT_DIR / "branches" / "0" + assert (b0 / "branch_zero_complete.txt").exists(), ( + "branch_zero_complete.txt missing" + ) + assert (b0 / "dem_0.tif").exists(), "dem_0.tif missing from branch 0" + assert (b0 / "flowdir_d8_burned_filled_0.tif").exists(), ( + "flowdir missing from branch 0" + ) + assert (b0 / "hydroTable_0.csv").exists(), "hydroTable_0.csv missing from branch 0" + assert result.n_branch_zero_recorded == 1, "branch zero not in branch_ids.csv" + assert result.branch_ids_csv.exists(), "branch_ids.csv not created" + + # branch_results now includes branch 0 at index 0. + b0_res = next((r for r in result.branch_results if r.branch_id == "0"), None) + assert b0_res is not None and b0_res.status == "ok", f"branch 0 status: {b0_res}" + + non_zero = [r for r in result.branch_results if r.branch_id != "0"] + ok = sum(1 for r in non_zero if r.status == "ok") + log.info(f"combined: branch_zero=ok non-zero ok={ok}/{len(non_zero)}") + assert result.n_non_zero_recorded == ok + + # Spot-check one non-zero branch hydroTable. + ok_branches = [r.branch_id for r in non_zero if r.status == "ok"] + if ok_branches: + sample_ht = ( + OUT_DIR / "branches" / ok_branches[0] / f"hydroTable_{ok_branches[0]}.csv" + ) + assert sample_ht.exists(), f"hydroTable missing from branch {ok_branches[0]}" + + +# ============================================================================= +# INDIVIDUAL STEP-BY-STEP TESTS +# Running these in file order rebuilds the full per-branch pipeline +# Layers: +# Z0 BranchDerivation — level paths + branch_list.csv +# Z1 BranchZero — DEM clip + AGREE + pit-fill + D8 +# (wraps stream raster, headwater +# raster, optional levelpath raster, +# optional levee burn, AGREE, +# pit-fill, flowdir) +# B02..B21 CreateHAND steps 2-21 — one test per CreateHAND step +# ============================================================================= + +# # Stage Z — bootstrap. Together they produce every input the B-series tests need. +# def test_step_Z0_branch_derivation(): +# """Derive level paths, branch polygons, and branch list from staged NWM data.""" +# result = BranchDerivation( +# out_dir=OUT_DIR, +# branch_id_attribute="levpa_id", +# reach_id_attribute="ID", +# branch_buffer_distance_meters=7000.0, +# ).run() +# assert result.dissolved_levelpaths.exists(), "dissolved levelpaths not written" +# assert result.branch_polygons.exists(), "branch polygons not written" +# assert result.branch_list.exists(), "branch list file not written" +# assert len(result.branch_dataframe) > 0, "branch dataframe is empty" +# log.info(f"branch count: {len(result.branch_dataframe)}") + + +# def test_step_Z1_branch_zero_full(): +# """Run BranchZero: DEM clip, stream rasterize, optional headwater/levelpath/levee +# rasters, AGREE conditioning, pit-fill, and D8 flowdir for branch 0. + +# This single call wraps the substeps BranchZero already folds together +# (StreamBooleanRasterizer, HeadwaterRasterizer, optional +# LevelPathBooleanRasterizer, optional rasterize_3d_levee_lines + +# burn_levee_elevations, HydroenforceDEM, WhiteboxTools pit-fill, +# FlowdirDEM). Calling the substeps individually would duplicate work the +# class already orchestrates correctly. +# """ +# outputs = BranchZero( +# dem_path=DEM, +# streams_gpkg=STREAMS, +# boundary_gpkg=BOUNDARY_BUF, +# out_dir=OUT_DIR, +# bridge_elev_diff_path=BRIDGE_DIFF if BRIDGE_DIFF.exists() else None, +# levee_gpkg_path=NLD_LEVEES if NLD_LEVEES.exists() else None, +# headwaters_gpkg=HEADWATERS if HEADWATERS.exists() else None, +# levelpaths_extended_gpkg=LEVELPATH_EXT if LEVELPATH_EXT.exists() else None, +# agree_buffer_m=15.0, +# agree_smooth_drop=10.0, +# agree_sharp_drop=1000.0, +# branch_zero_id=BRANCH_ID, +# ).run() +# for key, p in outputs.items(): +# log.info(f" {key:35s} --> {p.name}") +# assert DEM_BRANCH.exists(), "dem_0.tif missing" +# assert STREAM_BOOL.exists(), "flows_grid_boolean_0.tif missing" +# assert DEM_BURNED.exists(), "dem_burned_0.tif missing" +# assert DEM_FILLED.exists(), "dem_burned_filled_0.tif missing" +# assert FLOWDIR.exists(), "flowdir_d8_burned_filled_0.tif missing" + + +# # Stage B — CreateHAND steps 2..21, one isolated test each. +# def test_step_B02_flow_accumulation(): +# """CreateHAND step 2: D8 flow accumulation + stream-pixel mask.""" +# assert FLOWDIR.exists(), "FLOWDIR missing — run step_A6 first" +# if not HW_RASTER.exists(): +# log.warning("skipping flow accumulation — no headwater raster") +# return +# fa_out, sp_out = FlowAccDEM( +# flowdir=FLOWDIR, +# headwaters=HW_RASTER, +# out_flowaccum=FLOWACCUM, +# out_stream_pixels=STREAM_PIX, +# threshold=1.0, +# ).run() +# import rasterio + +# with rasterio.open(str(sp_out)) as src: +# stream_count = int((src.read(1) == 1).sum()) +# log.info(f"stream pixels: {stream_count}") +# assert fa_out.exists() and sp_out.exists() and stream_count > 0 + + +# def test_step_B03_thalweg_adjustment(): +# """CreateHAND step 3: lateral thalweg minimum + flow-conditioned DEM.""" +# for p in (DEM_BRANCH, STREAM_PIX, FLOWDIR): +# assert p.exists(), f"missing: {p}" +# result = ThalwegAdjustment( +# dem=DEM_BRANCH, +# stream_pixels=STREAM_PIX, +# flowdir=FLOWDIR, +# out_thalweg_adj=THALWEG_ADJ, +# out_flowdir_streams=FLOWDIR_STR, +# out_thalweg_cond=THALWEG_COND, +# cost_distance_tolerance=50.0, +# lateral_elevation_threshold=10, +# ).run() +# assert result["thalweg_adj"].exists() and result["thalweg_cond"].exists() + + +# def test_step_B04_d8_slopes(): +# """CreateHAND step 4: D8 slope raster (rise/run from thalweg-adjusted DEM).""" +# assert THALWEG_ADJ.exists() and FLOWDIR.exists() +# import numpy as np, rasterio + +# out = D8SlopeDEM( +# dem=THALWEG_ADJ, flowdir=FLOWDIR, out_path=SLOPES_D8, slope_min=0.0001 +# ).run() +# with rasterio.open(str(out)) as src: +# d = src.read(1) +# nd = src.nodata +# valid = d[(d != nd) & np.isfinite(d)] if nd is not None else d[np.isfinite(d)] +# log.info(f"slope range: [{valid.min():.6f}, {valid.max():.6f}]") +# # slope_min is clamped at 1e-4 in float32; allow a single-precision epsilon +# # of tolerance (~1e-7) so the test doesn't fail on the float32 representation +# # of 1e-4 (which is 9.9999997e-05). +# assert float(valid.min()) >= 0.0001 - 1e-7 + + +# def test_step_B05_streamnet_reaches(): +# """CreateHAND step 5: vectorise stream network into reach polylines.""" +# for p in (FLOWDIR, THALWEG_COND, FLOWACCUM, STREAM_PIX): +# assert p.exists(), f"missing: {p}" +# result = StreamNetReaches( +# flowdir=FLOWDIR, +# dem_thalweg_cond=THALWEG_COND, +# flowaccum=FLOWACCUM, +# stream_pixels=STREAM_PIX, +# out_dir=BRANCH_DIR, +# branch_id=BRANCH_ID, +# ).run() +# import geopandas as gpd + +# reaches = gpd.read_file(str(result["demDerived_reaches"])) +# log.info(f"reaches: {len(reaches)}") +# assert len(reaches) > 0 + + +# def test_step_B06_split_reaches(): +# """CreateHAND step 6: split reaches at length limit + lake boundaries.""" +# for p in (DEM_REACHES, THALWEG_COND, STREAMS): +# assert p.exists(), f"missing: {p}" +# split_gpkg, pts_gpkg = split_derived_reaches( +# reaches_gpkg=DEM_REACHES, +# dem_thalweg_cond=THALWEG_COND, +# nwm_streams_gpkg=STREAMS, +# out_split_gpkg=SPLIT_REACHES, +# out_points_gpkg=SPLIT_PTS, +# wbd8_clp_gpkg=WBD8_CLP if WBD8_CLP.exists() else None, +# lakes_gpkg=LAKES if LAKES.exists() else None, +# # This could be interesting point where based on slope or anyother logic, you can segment the reach--> ultimately gives the corresponsing +# # catchment, meaning shorter the reach length- denser the catchment +# max_length=1500.0, +# slope_min=0.0001, +# lakes_buffer_dist=100.0, +# ) +# import geopandas as gpd + +# split = gpd.read_file(str(split_gpkg)) +# log.info(f"split reaches: {len(split)} columns={list(split.columns)}") +# assert ( +# len(split) > 0 and "HydroID" in split.columns and "NextDownID" in split.columns +# ) + + +# def test_step_B07_gage_watershed_reaches(): +# """CreateHAND step 7: reverse-D8 walk labelling each pixel by its HydroID.""" +# from fimbox import GageCatchments + +# for p in (FLOWDIR, SPLIT_PTS): +# assert p.exists(), f"missing: {p}" +# # declutter=True mirrors CreateHAND step 7: solidify the reach raster +# # (fill pits, de-checkerboard, one piece per HydroID) so it polygonizes clean. +# GageCatchments( +# flowdir=FLOWDIR, +# outlet_points=SPLIT_PTS, +# out_path=GW_REACHES, +# declutter=True, +# ).run() +# assert GW_REACHES.exists() + + +# def test_step_B08_stream_pixel_points(): +# """CreateHAND step 8: vectorise stream-pixel centroids (one point per stream pixel).""" +# from fimbox import stream_pixel_points + +# assert STREAM_PIX.exists() +# stream_pixel_points(stream_pixels=STREAM_PIX, out_gpkg=PIXEL_PTS) +# assert PIXEL_PTS.exists() + + +# def test_step_B09_gage_watershed_pixels(): +# """CreateHAND step 9: reverse-D8 walk labelling each pixel by NWM feature_id.""" +# from fimbox import GageCatchments + +# for p in (FLOWDIR, PIXEL_PTS): +# assert p.exists(), f"missing: {p}" +# GageCatchments( +# flowdir=FLOWDIR, +# outlet_points=PIXEL_PTS, +# out_path=GW_PIXELS, +# ).run() +# assert GW_PIXELS.exists() + + +# def test_step_B10_outlet_backpool_mitigation(): +# """CreateHAND step 10: trim oversized outlet catchments (no-op for branch 0).""" +# from fimbox import OutletBackpoolMitigate + +# for p in (SPLIT_REACHES, GW_PIXELS, GW_REACHES, SPLIT_PTS, STREAMS, THALWEG_COND): +# assert p.exists(), f"missing: {p}" +# OutletBackpoolMitigate( +# branch_dir=BRANCH_DIR, +# catchment_pixels_path=GW_PIXELS, +# catchment_reaches_path=GW_REACHES, +# split_flows_gpkg=SPLIT_REACHES, +# split_points_gpkg=SPLIT_PTS, +# nwm_streams_gpkg=STREAMS, +# dem_path=THALWEG_COND, +# slope_min=0.0001, +# ).run() +# # No new file is asserted — backpool mitigation modifies the existing +# # gw_catchments_pixels/reaches rasters in place for non-zero branches only. +# assert GW_PIXELS.exists() and GW_REACHES.exists() + + +# def test_step_B11_make_rem(): +# """CreateHAND step 11: HAND = pixel_elev - nearest_stream_pixel_elev. + +# Note: the raw REM **can** be negative (pixels lower than the nearest +# downstream stream pixel — happens near floodplain edges and where the +# D8 walk crosses meander cutoffs). Negative values get clipped to zero +# in step 12 (``rem_zeroed_masked``). This test only asserts the raster +# was produced and contains finite values — it does NOT enforce +# non-negativity, which is a step-12 invariant. +# """ +# from fimbox import MakeREM + +# for p in (THALWEG_COND, GW_PIXELS, STREAM_PIX): +# assert p.exists(), f"missing: {p}" +# out = MakeREM( +# dem_thalweg_cond=THALWEG_COND, +# gw_catchments_pixels=GW_PIXELS, +# stream_pixels=STREAM_PIX, +# out_rem=REM, +# ).run() +# import rasterio, numpy as np + +# with rasterio.open(str(out)) as src: +# data = src.read(1) +# nd = src.nodata +# valid = data[data != nd] if nd is not None else data.ravel() +# log.info( +# f"REM range: [{float(valid.min()):.2f}, {float(valid.max()):.2f}] " +# f"({(valid < 0).sum()} negative pixels — clipped by step 12)" +# ) +# assert out.exists() and valid.size > 0 and np.isfinite(valid).all() + + +# def test_step_B11b_rem_nonnegative_after_zero_mask(): +# """Cross-check: after step 12 (rem_zeroed_masked), the REM raster must be +# non-negative and contain no NaN pixels. The reference formula +# ``(A * (A>=0) * (B>0))`` with an explicit NoDataValue treats NaN inputs as +# zero; the fimbox port now matches that behaviour by rewriting NaN to the +# nodata sentinel before the multiply. + +# Lives next to B11 so a failure here points at the zero-mask logic, not at +# MakeREM itself. Skipped silently if step 12 hasn't run yet (run B12 first). +# """ +# import numpy as np +# import rasterio + +# if not REM_ZEROED.exists(): +# log.warning("skipping non-negativity check — run step_B12 first") +# return +# with rasterio.open(str(REM_ZEROED)) as src: +# data = src.read(1) +# nd = src.nodata +# # Strip both the nodata sentinel and any NaN before the min() so the +# # test catches the actual data range, not an IEEE NaN propagating. +# if nd is not None: +# valid_mask = (data != nd) & ~np.isnan(data) +# else: +# valid_mask = ~np.isnan(data) +# valid = data[valid_mask] +# nan_count = int(np.isnan(data).sum()) +# log.info(f"REM zero-mask: {valid.size} valid pixels, {nan_count} NaN pixels") +# assert valid.size > 0 +# assert ( +# nan_count == 0 +# ), f"step 12 leaked {nan_count} NaN pixels into the masked REM raster" +# assert ( +# float(valid.min()) >= 0.0 +# ), f"step 12 left negatives in REM: min={valid.min()}" + + +# def test_step_B12_rem_zeroed_masked(): +# """CreateHAND step 12: clip negative HAND to 0 + mask outside catchments.""" +# from fimbox import rem_zeroed_masked + +# for p in (REM, GW_REACHES): +# assert p.exists(), f"missing: {p}" +# rem_zeroed_masked(REM, GW_REACHES, REM_ZEROED) +# assert REM_ZEROED.exists() + + +# def test_step_B13_polygonize_catchments(): +# """CreateHAND step 13: rasterised catchments --> per-HydroID polygon gpkg.""" +# # Helper lives inside create_hand.py as a private function; import it explicitly. +# from fimbox.preprocessing.calculate_branch.create_hand import ( +# _polygonize_catchments, +# ) + +# assert GW_REACHES.exists() +# _polygonize_catchments(GW_REACHES, CATCH_POLY) +# import geopandas as gpd + +# gdf = gpd.read_file(str(CATCH_POLY)) +# log.info(f"polygonised: {len(gdf)} catchments") +# assert CATCH_POLY.exists() and "HydroID" in gdf.columns and len(gdf) > 0 + + +# def test_step_B14_filter_catchments(): +# """CreateHAND step 14: drop slivers + attach flow attributes per HydroID.""" +# from fimbox import FilterCatchments + +# for p in (CATCH_POLY, SPLIT_REACHES): +# assert p.exists(), f"missing: {p}" +# out_catch, out_flows = FilterCatchments( +# catchments_gpkg=CATCH_POLY, +# flows_gpkg=SPLIT_REACHES, +# out_catchments=FILT_CATCH, +# out_flows=FILT_FLOWS, +# aoi_code=OUT_DIR.parent.name, +# boundary_gpkg=WBD8_CLP if WBD8_CLP.exists() else None, +# ).run() +# import geopandas as gpd + +# catches = gpd.read_file(str(out_catch)) +# flows = gpd.read_file(str(out_flows)) +# log.info(f"filtered catchments: {len(catches)} flows: {len(flows)}") +# assert len(catches) > 0 and "areasqkm" in catches.columns +# assert len(flows) > 0 and "HydroID" in flows.columns + + +# def test_step_B15_rasterize_filtered_catchments(): +# """CreateHAND step 15: burn HydroID back onto the reference raster grid.""" +# from fimbox.preprocessing.calculate_branch.create_hand import ( +# _rasterize_catchments, +# ) + +# for p in (FILT_CATCH, GW_REACHES): +# assert p.exists(), f"missing: {p}" +# _rasterize_catchments(FILT_CATCH, GW_REACHES, FILT_TIF) +# assert FILT_TIF.exists() + + +# def test_step_B16_mask_slopes_to_catchments(): +# """CreateHAND step 16: clip D8 slopes to the filtered catchment mask.""" +# from fimbox import mask_slopes_to_catchments + +# for p in (SLOPES_D8, FILT_TIF): +# assert p.exists(), f"missing: {p}" +# mask_slopes_to_catchments(SLOPES_D8, FILT_TIF, SLOPES_MASKED) +# assert SLOPES_MASKED.exists() + + +# def test_step_B17_stages_and_catchlist(): +# """CreateHAND step 17: write the stage ladder + per-HydroID metadata text files.""" +# from fimbox import make_stages_and_catchlist + +# for p in (FILT_FLOWS, FILT_CATCH): +# assert p.exists(), f"missing: {p}" +# make_stages_and_catchlist( +# flows_gpkg=FILT_FLOWS, +# catchments_gpkg=FILT_CATCH, +# out_stages=STAGE_TXT, +# out_catchlist=CATCHLIST_TXT, +# stages_min=0.0, +# stages_interval=0.3048, +# stages_max=25.2984, +# ) +# assert STAGE_TXT.exists() and CATCHLIST_TXT.exists() + + +# def test_step_B18_build_src_base(): +# """CreateHAND step 18: synthetic rating curve base table.""" +# from fimbox import build_src_base + +# for p in (REM_ZEROED, FILT_TIF, SLOPES_MASKED, CATCHLIST_TXT, STAGE_TXT): +# assert p.exists(), f"missing: {p}" +# build_src_base( +# hand_raster=REM_ZEROED, +# catch_raster=FILT_TIF, +# slope_raster=SLOPES_MASKED, +# catchlist_txt=CATCHLIST_TXT, +# stages_txt=STAGE_TXT, +# out_csv=SRC_BASE_CSV, +# ) +# import pandas as pd + +# df = pd.read_csv(SRC_BASE_CSV) +# log.info(f"src_base: {len(df)} rows HydroIDs={df['CatchId'].nunique()}") +# assert SRC_BASE_CSV.exists() and len(df) > 0 + + +# def test_step_B19_add_crosswalk(): +# """CreateHAND step 19: NWM crosswalk + Manning's hydraulics + hydroTable.""" +# from fimbox import add_crosswalk + +# for p in (FILT_CATCH, FILT_FLOWS, SRC_BASE_CSV, STREAMS): +# assert p.exists(), f"missing: {p}" +# add_crosswalk( +# catchments_gpkg=FILT_CATCH, +# flows_gpkg=FILT_FLOWS, +# src_base_csv=SRC_BASE_CSV, +# nwm_streams_gpkg=STREAMS, +# out_catchments_gpkg=XWALK_CATCH, +# out_flows_gpkg=XWALK_FLOWS, +# out_src_csv=SRC_FULL_CSV, +# out_src_json=SRC_JSON, +# out_crosswalk_csv=XWALK_CSV, +# out_hydro_csv=HYDRO_TABLE, +# boundary_gpkg=WBD8_CLP if WBD8_CLP.exists() else None, +# mannings_n=0.06, +# min_catchment_area=0.25, +# min_stream_length=0.5, +# max_distance_m=100.0, +# small_segments_csv=BRANCH_DIR / f"small_segments_{BRANCH_ID}.csv", +# # SRC slope source (optional): "iris_sword" (default) | "dem" | "hfab". +# src_slope_source="iris_sword", +# iris_slope_csv=None, # None -> packaged fimbox/data table +# hfab_slope_column=None, # name the hydrofabric slope col if not Slope/So +# ) +# import pandas as pd + +# ht = pd.read_csv(HYDRO_TABLE, dtype={"HydroID": str}) +# log.info(f"hydroTable: {len(ht)} rows HydroIDs={ht['HydroID'].nunique()}") +# assert HYDRO_TABLE.exists() and (ht["discharge_cms"] >= 0).all() +# # The three slope variants are carried so the chosen source is transparent. +# for col in ("SLOPE", "SLOPE_RISE_RUN", "SLOPE_IRIS_SWORD"): +# assert col in ht.columns, f"hydroTable missing {col}" + + +# def test_step_B20_heal_bridges_osm(): +# """CreateHAND step 20: raise HAND at OSM bridge decks (in-place REM update).""" +# from fimbox import heal_bridges_osm + +# bridges_gpkg = OUT_DIR / "osm_bridges_subset.gpkg" +# if not bridges_gpkg.exists(): +# log.warning("skipping bridge heal — no OSM bridges gpkg") +# return +# for p in (REM_ZEROED, XWALK_CATCH): +# assert p.exists(), f"missing: {p}" +# bridge_diff = OUT_DIR / "bridge_elev_diff.tif" +# heal_bridges_osm( +# hand_raster=REM_ZEROED, +# bridges_gpkg=bridges_gpkg, +# catchments_gpkg=XWALK_CATCH, +# out_centroids_gpkg=BRIDGES_GPKG, +# bridge_diff_raster=bridge_diff if bridge_diff.exists() else None, +# ) +# assert BRIDGES_GPKG.exists() + + +# def test_step_B21_process_roads_fimpact(): +# """CreateHAND step 21: sample HAND along OSM roads to derive flood thresholds.""" +# from fimbox import process_roads_fimpact + +# roads_gpkg = OUT_DIR / "osm_roads_subset.gpkg" +# if not roads_gpkg.exists(): +# log.warning("skipping road FIMpact — no OSM roads gpkg") +# return +# for p in (REM_ZEROED, XWALK_CATCH): +# assert p.exists(), f"missing: {p}" +# process_roads_fimpact( +# hand_raster=REM_ZEROED, +# roads_gpkg=roads_gpkg, +# catchments_gpkg=XWALK_CATCH, +# out_csv=ROADS_CSV, +# ) +# assert ROADS_CSV.exists() + + +# # Stage C — branch-zero post-CreateHAND steps +# # (download USGS gauges --> AOI-level assignment --> branch-zero crosswalk --> cleanup) + +# # AOI-level path to the staged USGS gages gpkg +# USGS_GAGES = OUT_DIR / "usgs_gages.gpkg" +# USGS_SUBSET = OUT_DIR / "usgs_subset_gages.gpkg" +# USGS_SUBSET_BZERO = OUT_DIR / f"usgs_subset_gages_{BRANCH_ID}.gpkg" +# NWM_LEVELPATHS = OUT_DIR / f"{IDENTIFIER}_subset_streams_levelPaths.gpkg" + + +# def test_step_C20_download_usgs_gages(): +# """Download USGS gauge points inside the AOI from the ArcGIS Online +# FeatureServer. Writes ``usgs_gages.gpkg`` at the AOI root, with the columns +# ``assign_gages_to_branches`` expects: ``location_id``, ``feature_id``, +# ``aoi_id``, ``source``, geometry. +# """ +# from fimbox import DownloadUSGSGages + +# # Use the buffered boundary so gauges just outside the WBD are still +# # captured (they may snap to streams that drain into the AOI). +# boundary = BOUNDARY_BUF if BOUNDARY_BUF.exists() else WBD8_CLP +# assert boundary.exists(), f"missing boundary: {boundary}" + +# gdf = DownloadUSGSGages().download( +# boundary=boundary, +# aoi_id=OUT_DIR.parent.name, +# out_dir=OUT_DIR, +# out_name="usgs_gages.gpkg", +# out_layer="usgs_gages", +# ) +# log.info(f"USGS gauges downloaded: {len(gdf)} features --> {USGS_GAGES.name}") +# # Empty AOI (no gauges in CONUS layer) is acceptable; only assert the +# # file exists when at least one feature came back. +# if len(gdf) > 0: +# assert USGS_GAGES.exists() +# assert {"location_id", "feature_id", "aoi_id", "source"}.issubset(gdf.columns) + + +# def test_step_C21_assign_gages_to_branches(): +# """Stage 1 of the gage crosswalk: tag every gage with a ``feature_id`` + +# ``levpa_id`` (= branch id) and write the AOI-wide + branch-zero gpkgs. + +# Skips if either ``usgs_gages.gpkg`` (from C20) or +# ``nwm_subset_streams_levelPaths.gpkg`` (from BranchDerivation in Z0) is +# missing — both prerequisites get logged so a failure points at the +# right upstream step. +# """ +# from fimbox import assign_gages_to_branches + +# if not USGS_GAGES.exists(): +# log.warning( +# "skipping gage assignment — usgs_gages.gpkg missing (run step_C20 first)" +# ) +# return +# if not NWM_LEVELPATHS.exists(): +# log.warning( +# "skipping gage assignment — nwm_subset_streams_levelPaths.gpkg missing " +# "(run step_Z0_branch_derivation first)" +# ) +# return + +# assign_gages_to_branches( +# usgs_gages_gpkg=USGS_GAGES, +# nwm_streams_levelpaths_gpkg=NWM_LEVELPATHS, +# aoi_id=OUT_DIR.parent.name, +# out_dir=OUT_DIR, +# # DownloadUSGSGages writes "aoi_id"; the default filter column ("HUC8") +# # would not find anything in that gpkg. +# aoi_filter_column="aoi_id", +# branch_zero_id=BRANCH_ID, +# ) +# # When the AOI actually contains gauges both files exist; on empty AOIs +# # neither is written and the function returns None (logged a warning). +# if USGS_SUBSET.exists(): +# log.info( +# f"AOI-wide gages --> {USGS_SUBSET.name} | " +# f"branch-zero --> {USGS_SUBSET_BZERO.name}" +# ) +# assert USGS_SUBSET_BZERO.exists() + + +# def test_step_C22_usgs_crosswalk_branch_zero(): +# """Stage 2 of the gage crosswalk for branch zero. + +# Snaps every branch-zero gage to its DEM-derived thalweg and samples the +# DEM + thalweg-conditioned DEM to populate ``dem_elevation`` and +# ``dem_adj_elevation`` on the gage table. Output: +# ``branches/0/usgs_elev_table.csv``. + +# Prerequisites: ``usgs_subset_gages_0.gpkg`` (from C21) and the per-branch +# CreateHAND outputs (from Z1 + the B-series). +# """ +# from fimbox import run_branch_crosswalk + +# if not USGS_SUBSET_BZERO.exists(): +# log.warning( +# "skipping USGS crosswalk — usgs_subset_gages_0.gpkg missing " +# "(run step_C20 + step_C21 first to produce it)" +# ) +# return +# bzero_gages = USGS_SUBSET_BZERO + +# # dem_meters_{B}.tif is the inundation-mapping name; fimbox writes dem_{B}.tif +# # via BranchZero. Use whichever exists. +# dem_b = BRANCH_DIR / f"dem_meters_{BRANCH_ID}.tif" +# if not dem_b.exists(): +# dem_b = DEM_BRANCH + +# for p in (XWALK_CATCH, FILT_FLOWS, dem_b, THALWEG_COND): +# assert p.exists(), f"missing: {p}" + +# out = run_branch_crosswalk( +# aoi_gages_gpkg=bzero_gages, +# branch_catchments_gpkg=XWALK_CATCH, +# branch_flows_gpkg=FILT_FLOWS, +# dem_path=dem_b, +# dem_thalweg_path=THALWEG_COND, +# branch_id=BRANCH_ID, +# out_dir=BRANCH_DIR, +# ) +# usgs_table = BRANCH_DIR / "usgs_elev_table.csv" +# log.info(f"USGS crosswalk wrote: {[p for p in out.values() if p]}") +# # usgs_elev_table.csv only exists when the AOI has gages — log either way. +# if usgs_table.exists(): +# import pandas as pd + +# df = pd.read_csv(usgs_table) +# log.info(f"usgs_elev_table.csv rows: {len(df)}") + + +# def test_step_C23_outputs_cleanup_branch_zero(): +# """Apply the deny-list cleanup to ``branches/0/``. + +# Default behaviour deletes every intermediate raster + vector listed in +# --> fimbox/config/deny_branch_zero.lst. +# """ +# import os + +# from fimbox import remove_deny_list_files + +# deny_path = ( +# Path(__file__).resolve().parent.parent / "config" / "deny_branch_zero.lst" +# ) +# assert deny_path.is_file(), f"deny list missing: {deny_path}" + +# # API sanity checks that always run (never touch real files). +# assert remove_deny_list_files(BRANCH_DIR, "NONE", BRANCH_ID) == 0 +# assert remove_deny_list_files(BRANCH_DIR, "none", BRANCH_ID) == 0 + +# if os.environ.get("FIMBOX_KEEP_BRANCH_ZERO"): +# n_patterns = sum( +# 1 +# for L in deny_path.read_text().splitlines() +# if L.strip() and not L.lstrip().startswith("#") +# ) +# log.info( +# f"step_C23: skipping cleanup (FIMBOX_KEEP_BRANCH_ZERO set). " +# f"{deny_path.name} has {n_patterns} active patterns; " +# "unset the env var to enable cleanup." +# ) +# return + +# # The branch-0 directory may be empty when only later steps have been +# # populated, or when an earlier C23 run already cleaned it. Skip cleanly +# # if there's nothing to do. +# if not BRANCH_DIR.exists(): +# log.warning(f"skipping cleanup — branch dir {BRANCH_DIR} missing") +# return + +# n = remove_deny_list_files( +# src_dir=BRANCH_DIR, +# deny_list=deny_path, +# branch_id=BRANCH_ID, +# verbose=True, +# ) +# log.info(f"step_C23: removed {n} files from {BRANCH_DIR}") + + +# def test_step_C24_calculate_allbranches(tmp_path): +# """Fast wrapper check without launching real branch workers.""" +# from fimbox import AOIProcessingConfig, calculate_allbranches + +# aoi_dir = tmp_path / "aoi" +# aoi_dir.mkdir() +# # Match BranchDerivation's actual output: branch_ids.lst (one id per line). +# # Empty file = branch-zero-only run, which is what this wrapper test exercises. +# branch_list_path = aoi_dir / "branch_ids.lst" +# branch_list_path.write_text("") + +# deny_unit_list = tmp_path / "deny_unit.lst" +# deny_unit_list.write_text("temporary_{}.tif\n") +# # aoi_id defaults to the AOI folder name ("aoi") when not passed. +# removable = aoi_dir / f"temporary_{aoi_dir.name}.tif" +# removable.write_bytes(b"x") + +# result = calculate_allbranches( +# AOIProcessingConfig( +# aoi_dir=aoi_dir, +# branch_list_path=branch_list_path, +# n_workers=1, +# ), +# delete_deny_list=True, +# deny_unit_list=deny_unit_list, +# branch_ids_csv=aoi_dir / "branch_ids.csv", +# ) + +# assert result.n_branch_zero_recorded == 1 +# assert result.n_non_zero_recorded == 0 +# assert result.n_unit_files_removed == 1 +# assert not removable.exists() + + +# def test_step_C25_calculate_allbranches_live_run(): +# """Live run for the real non-zero branch loop. + +# Set FIMBOX_KEEP_UNIT=1 to skip AOI-level cleanup. +# Set FIMBOX_SKIP_ALLBRANCHES=1 to skip this test (e.g. during quick CI +# smoke runs); by default it always runs. +# """ +# from fimbox import AOIProcessingConfig, calculate_allbranches + +# if os.environ.get("FIMBOX_SKIP_ALLBRANCHES"): +# pytest.skip("FIMBOX_SKIP_ALLBRANCHES set — skipping live branch loop") + +# # BranchDerivation writes branch_ids.lst +# branch_list_path = OUT_DIR / "branch_ids.lst" + +# deny_unit_list = Path(__file__).resolve().parent.parent / "config" / "deny_unit.lst" +# assert deny_unit_list.is_file(), f"deny_unit.lst missing: {deny_unit_list}" + +# # Reuse branch-zero tuning so both paths stay in sync. +# # Auto-size workers to the machine; set FIMBOX_DASK_WORKERS=1 for serial. +# from fimbox._dask import _resolve_n_workers + +# n_workers = _resolve_n_workers() +# log.info(f"Branch processing with n_workers={n_workers}") + +# cfg = AOIProcessingConfig( +# aoi_dir=OUT_DIR, +# branch_list_path=branch_list_path, +# n_workers=n_workers, # auto-sized; FIMBOX_DASK_WORKERS=1 forces serial +# keep_failed_branches=True, # keep a failed branch dir for inspection +# delete_deny_list=True, +# **PARAMS_CREATE_HAND, +# ) + +# delete_deny_list = True +# result = calculate_allbranches( +# cfg, +# delete_deny_list=delete_deny_list, +# deny_unit_list=deny_unit_list if delete_deny_list else None, +# branch_ids_csv=OUT_DIR / "branch_ids.csv", +# ) + +# assert result.n_branch_zero_recorded == 1 +# assert result.branch_ids_csv.exists(), "branch_ids.csv was not created" +# assert result.n_non_zero_recorded == sum( +# 1 for r in result.branch_results if r.status == "ok" +# ) diff --git a/tests/test_calibrate_pipeline.py b/tests/test_calibrate_pipeline.py new file mode 100644 index 0000000..7212435 --- /dev/null +++ b/tests/test_calibrate_pipeline.py @@ -0,0 +1,289 @@ +""" +Author: Supath Dhital +Date Updated: June 2026 + +Tests for the synthetic rating curve (SRC) calibration pipeline. + +Two layers: + + COMBINED ......... test_calibrate_full_pipeline runs the whole thing in a + single run_calibration() call against the live AOI, with EVERY optional + CalibrationConfig parameter spelled out so the full surface is visible + in one place. + + STEP BY STEP ..... one test per stage (thalweg, longitudinal, bathymetry, + bankfull, subdiv, nonmonotonic, usgs, spatial, log scan) so any single + step can be run / debugged alone. + +It will point into the working version of the AOI and skip when it is absent. +""" + +from __future__ import annotations + +from pathlib import Path + +import pandas as pd +import pytest + +from fimbox import CalibrationConfig, run_calibration +from fimbox._dask import _resolve_n_workers +from fimbox.datasets import fetch_data + +# Live AOI + input files. Edit these to point at your data; tests skip when the AOI is absent. +AOI_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" + +# Calibration lookup tables, fetched on demand from the public SDML S3 bucket +# (anonymous, cached locally) via the pooch registry in ``fimbox.datasets``. + +# Bankfull recurrence flows (NWM v3) +BANKFULL_FLOWS_FILE = fetch_data("nwm3_high_water_threshold") + +# Optimized variable-roughness Manning's n table (per feature_id channel/overbank n) +VMANN_INPUT_FILE = fetch_data("mannings_optz") + +# USGS rating-curve calibration. Rating curve + NWM recurrence flows (v3) are +# required; the acceptable-gage quality filter refines which gages qualify. +USGS_RATING_CURVE_CSV = fetch_data("usgs_rating_curves") +NWM_RECUR_FILE = fetch_data("nwm3_recurrence_flows") +USGS_ACCEPTABLE_GAGES = fetch_data("acceptable_gages") + +# Bathymetry: eHydro surveyed channels (.gpkg). +BATHY_EHYDRO_FILE = fetch_data("bathymetry_ehydro_ohrfc") + +# Spatial-observation calibration: per-AOI benchmark points (.parquet). Yet to +# be added — left unset for now. +CALIB_POINTS_FILE = None + +# Manual calibration: per-feature_id coefficient CSV. Yet to be added. +MAN_CALB_FILE = None + +# Worker count for the branch-parallel routines — auto-sized to the device +JOB_BRANCH_LIMIT = _resolve_n_workers() + +_BRANCHES = ( + AOI_DIR / "watershed-data" / "branches" + if (AOI_DIR / "watershed-data" / "branches").is_dir() + else AOI_DIR / "branches" +) +_skip_no_aoi = pytest.mark.skipif( + not _BRANCHES.is_dir(), reason=f"AOI not present: {_BRANCHES}" +) +_skip_no_bankfull = pytest.mark.skipif( + not BANKFULL_FLOWS_FILE.is_file(), reason=f"missing {BANKFULL_FLOWS_FILE}" +) +_skip_no_bathy = pytest.mark.skipif( + BATHY_EHYDRO_FILE is None or not Path(BATHY_EHYDRO_FILE).is_file(), + reason=f"bathy eHydro file not set: {BATHY_EHYDRO_FILE}", +) +_skip_no_usgs = pytest.mark.skipif( + not USGS_RATING_CURVE_CSV.is_file(), reason=f"missing {USGS_RATING_CURVE_CSV}" +) +_skip_no_manual = pytest.mark.skipif( + MAN_CALB_FILE is None or not Path(MAN_CALB_FILE).is_file(), + reason=f"manual calib file not set: {MAN_CALB_FILE}", +) +_skip_no_spatial = pytest.mark.skipif( + CALIB_POINTS_FILE is None or not Path(CALIB_POINTS_FILE).is_file(), + reason=f"spatial calib points not set: {CALIB_POINTS_FILE}", +) + + +# COMBINED — the whole calibration pipeline in one call, matching the step-by-step +# sequence: reset -> aggregate_pre -> thalweg -> longitudinal -> bathymetry -> +# bankfull -> subdiv -> nonmonotonic -> usgs -> spatial -> manual -> aggregate_post -> log_scan. +# File-dependent steps (bathy, usgs, spatial, manual) self-skip when their +# input file is absent, matching the step-by-step skip decorators. +@_skip_no_aoi +@_skip_no_bankfull +def test_calibrate_full_pipeline(): + """One run_calibration() call driving the full default pipeline. + Every CalibrationConfig parameter is spelled out, grouped by step.""" + cfg = CalibrationConfig( + # reset — revert hydroTables to uncalibrated baseline before re-applying. + # Set True when re-calibrating an AOI that was already calibrated. + calibration_rerun=True, + # aggregate_pre — assemble usgs/ras2fim elev tables before adjustments + aggregate_pre=True, + # thalweg — remove thalweg-notch artifact rows, refill stage ladder + thalweg_notches_adjustment=True, + # longitudinal — smooth hydraulic geometry along reach chains + longitudinal_filter=True, + # bathymetry — add missing in-channel area below the DEM (needs bathy_file_ehydro) + bathymetry_adjust=True, + bathy_file_ehydro=BATHY_EHYDRO_FILE, + # bankfull — identify bankfull stage in every branch SRC + src_bankfull_toggle=True, + bankfull_flows_file=BANKFULL_FLOWS_FILE, + include_branch_zero=True, + # subdiv — channel/overbank subdivision (needs vmann + bankfull on) + src_subdiv_toggle=True, + vmann_input_file=VMANN_INPUT_FILE, + default_channel_n=0.06, # used when feature_id missing from vmann table + default_overbank_n=0.12, + # nonmonotonic — force monotonic in-channel rating curves + nonmonotonic_src_adjustment=True, + nonmonotonic_stream_order_min=4, + # usgs — calibrate SRCs against USGS rating curves at NWM recurrence flows + src_adjust_usgs=True, + usgs_rating_curve_csv=USGS_RATING_CURVE_CSV, + usgs_acceptable_gages=USGS_ACCEPTABLE_GAGES, + nwm_recur_file=NWM_RECUR_FILE, + # spatial — calibrate SRCs against benchmark inundation points + src_adjust_spatial=True, + calib_points_file=CALIB_POINTS_FILE, # None -> step self-skips + # manual — apply a per-feature_id coefficient table + manual_calb_toggle=True, + man_calb_file=MAN_CALB_FILE, # None -> step self-skips + # aggregate_post — publish htable + bridge + road to AOI root + aggregate_post=True, + # log scan — collect error/warning lines into per-AOI summary files + scan_logs=True, + # execution + job_branch_limit=JOB_BRANCH_LIMIT, + skip_unimplemented=True, # warn instead of raising on stubs + ) + run_calibration(AOI_DIR, cfg) + + # Subdivision rewrites the per-branch hydroTable with subdiv columns. + sample_ht = next(_BRANCHES.glob("*/hydroTable_*.csv")) + cols = pd.read_csv(sample_ht, nrows=1).columns + assert "subdiv_discharge_cms" in cols + assert "channel_n" in cols + + +# # STEP BY STEP — each stage on its own. +# @_skip_no_aoi +# def test_step_reset(): +# """Reset per-branch hydroTable + src_full_crosswalked to baseline. +# Needed only for reruns; a no-op on a fresh AOI. Runs before aggregation.""" +# HydroTableReset(aoi_dir=AOI_DIR).run() + + +# @_skip_no_aoi +# def test_step_aggregate_pre(): +# """Pre-calibration aggregation: usgs/ras2fim elev tables if available (not integrated yet) -> AOI root.""" +# BranchAggregator(aoi_dir=AOI_DIR, usgs_elev=True, ras_elev=True).run() + + +# @_skip_no_aoi +# def test_step_thalweg_notches(): +# """Remove thalweg-notch artifact rows and refill the stage ladder.""" +# results = ThalwegNotchesAdjustment( +# aoi_dir=AOI_DIR, +# n_workers=JOB_BRANCH_LIMIT, # branch-parallel +# stage_interval_m=0.3048, # SRC stage step +# n_stages=84, # full ladder length +# extrap_rows=3, # trailing rows fit for extrapolation +# ).run() +# assert results + + +# @_skip_no_aoi +# def test_step_longitudinal(): +# """Smooth hydraulic geometry along reach chains, recompute discharge.""" +# results = LongitudinalFlowFilter( +# aoi_dir=AOI_DIR, n_workers=JOB_BRANCH_LIMIT, n_stages=84 +# ).run() +# assert results + + +# @_skip_no_aoi +# @_skip_no_bathy +# def test_step_bathymetry(): +# """Add missing in-channel area below the DEM from eHydro surveys, then +# recompute discharge.""" +# results = BathymetricAdjustment( +# aoi_dir=AOI_DIR, +# bathy_file_ehydro=BATHY_EHYDRO_FILE, +# ).run() +# assert results + + +# @_skip_no_aoi +# @_skip_no_bankfull +# def test_step_bankfull(): +# """Identify bankfull stage in every branch SRC.""" +# results = SrcBankfull( +# aoi_dir=AOI_DIR, +# bankfull_flows_file=BANKFULL_FLOWS_FILE, +# n_workers=JOB_BRANCH_LIMIT, +# include_branch_zero=True, +# ).run() +# assert results # dict of branch_id -> status string + + +# @_skip_no_aoi +# @_skip_no_bankfull +# def test_step_subdiv(): +# """Channel/overbank subdivision. Depends on bankfull having run, so run +# it first within this test to keep the step self-contained.""" +# SrcBankfull( +# aoi_dir=AOI_DIR, bankfull_flows_file=BANKFULL_FLOWS_FILE, n_workers=1 +# ).run() +# results = SrcSubdiv( +# aoi_dir=AOI_DIR, +# vmann_table=VMANN_INPUT_FILE, +# n_workers=JOB_BRANCH_LIMIT, +# default_channel_n=0.06, # used when feature_id missing from vmann table +# default_overbank_n=0.12, +# ).run() +# assert results + + +# @_skip_no_aoi +# def test_step_nonmonotonic(): +# """Force monotonic in-channel rating curves.""" +# results = SrcNonmonotonic( +# aoi_dir=AOI_DIR, stream_order_min=4, include_branch_zero=True +# ).run() +# assert results + + +# @_skip_no_aoi +# @_skip_no_usgs +# def test_step_usgs(): +# """Calibrate SRCs against USGS rating curves at NWM recurrence flows. +# Needs usgs_elev_table.csv at the AOI root; self-skips when inputs are absent.""" +# results = UsgsRatingCalibrator( +# aoi_dir=AOI_DIR, +# usgs_rating_curve_csv=USGS_RATING_CURVE_CSV, +# usgs_acceptable_gages=USGS_ACCEPTABLE_GAGES, +# nwm_recur_file=NWM_RECUR_FILE, +# n_workers=JOB_BRANCH_LIMIT, +# ).run() +# assert results is not None + + +# @_skip_no_aoi +# @_skip_no_spatial +# def test_step_spatial(): +# """Calibrate SRCs against benchmark inundation points. Samples HAND/HydroID +# rasters at each point; self-skips when the points file is absent.""" +# results = SpatialObsCalibrator( +# aoi_dir=AOI_DIR, +# calib_points_file=CALIB_POINTS_FILE, +# n_workers=JOB_BRANCH_LIMIT, +# ).run() +# assert results is not None + + +# @_skip_no_aoi +# @_skip_no_manual +# def test_step_manual(): +# """Apply a per-feature_id coefficient table to each branch hydroTable. +# Needs MAN_CALB_FILE (aoi_id, feature_id, calb_coef_manual); no-op when the +# AOI has no matching entry.""" +# ManualCalibrator(aoi_dir=AOI_DIR, calibration_file=MAN_CALB_FILE).run() + + +# @_skip_no_aoi +# def test_step_aggregate_post(): +# """Post-calibration aggregation: htable + bridge + road -> AOI root.""" +# BranchAggregator(aoi_dir=AOI_DIR, htable=True, bridge=True, road=True).run() + + +# @_skip_no_aoi +# def test_step_log_scan(): +# """Scan logs/ for error / warning lines into per-AOI summary files.""" +# out = LogScanner(aoi_dir=AOI_DIR, calibration_rerun=False).run() +# assert set(out) == {"errors", "warnings"} diff --git a/tests/test_downloaddata.py b/tests/test_downloaddata.py new file mode 100644 index 0000000..b42987e --- /dev/null +++ b/tests/test_downloaddata.py @@ -0,0 +1,203 @@ +# Example Usage: +from pathlib import Path + +import fimbox + +PKG_ROOT = Path(__file__).resolve().parents[1] +REPO_ROOT = Path(__file__).resolve().parents[2] + +test_boundary = PKG_ROOT / "docs" / "test_boundary" / "test_smallB.shp" +OUT_DIR = REPO_ROOT / "out" + +# # Testing the entire NHDPlus data extraction process along with National Flood Hazard Layer data extraction +# # This is OLDER VERSION using EPA AWS S3 Bucket which will get for whole HUC6 region--> not very effective +# def test_getNHDdata(): +# nhd_data = fimbox.getNHDPlusData( +# NHDglobalBoundary = NHDboundary, #Contains all NHDPlus VPU/RPU boundaries +# # inputs_dir = None, #Directory to save input data, if None, direct folder directory will be created +# boundary_path = test_boundary, #Path to the boundary shapefile for which NHDPlus data is to be extracted, OR HUC8 ID +# # huc8: Optional[str] = None, +# # epsg: Optional[int] = None, +# # out_dir: Optional[str] = None, #Directory to save output data, if None, direct folder directory will be created +# # auto_run= True +# ) +# # nhd_data.process_flowlines() +# print(f"Process successful!") + +# def test_get_nfhl(): +# fimbox.DownloadFEMANFHL( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# out_name="fema_nfhl.gpkg", +# # log_path=None, +# ) + +# def test_download_nld(): +# fimbox.DownloadNLD( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# lines_name="NLD_Lines.gpkg", # default; override as needed +# polys_name="NLD_Polygons.gpkg", # default; override as needed +# ) + +##This is for the medium range +# def test_get_nhddata(): +# fimbox.NWMFlowlinesDownloader().download( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# out_name="nwm_subset_streams.gpkg", +# out_layer="flowlines", +# ) + +##Medium range +# def test_get_catchments(): +# fimbox.NWMCatchmentsDownloader().download( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# out_name="nwm_subset_catchments.gpkg", +# out_layer="catchments", +# ) + +# def test_get_lakes(): +# fimbox.NWMLakesDownloader().download( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# out_name="nwm_subset_lakes.gpkg", +# out_layer="lakes", +# ) + + +# # Get all NHD Plus Data +# def test_get_nhd_all(): +# fimbox.getNHDPlusData( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# download_flowlines=True, +# download_catchments=True, +# download_lakes=True, +# resolution="medium", # "high" -> NHDPlus HR flowlines/catchments via pynhd; "medium" (default) -> NWM. Lakes always NWM. +# identifier="nwmmr", # filename prefix; default "nwm" -> nwm_subset_streams.gpkg etc. +# ) + + +# High-resolution flowlines + catchments only (NHDPlus HR via pynhd). +# def test_get_nhd_hr(): +# fimbox.getNHDPlusHRData( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# download_flowlines=True, +# download_catchments=True, +# identifier="nwm", # prefix for saved files +# ) + + +# Bring-your-own flowlines/catchments: map your column names to the canonical +# schema (streams: ID, order_, levpa_id, feature_id[=ID]; catchments: ID). +# def test_normalize_byo_flowlines_catchments(): +# fl = fimbox.normalize_flowlines( +# "path/to/my_flowlines.gpkg", +# field_map={"ID": "nhdplusid", "order_": "streamorde", "levpa_id": "levelpathi"}, +# ) +# cat = fimbox.normalize_catchments( +# "path/to/my_catchments.gpkg", field_map={"ID": "nhdplusid"} +# ) +# assert {"ID", "order_", "levpa_id", "feature_id"}.issubset(fl.columns) +# assert "ID" in cat.columns + + +# Download + process a 3DEP DEM. Reads only the AOI window straight from the +# Planetary Computer COGs over HTTP (no national VRT parse), one snapped +# reprojection, then hole-fill + clip. ~8x faster than the old py3dep path. +# Resolutions: 10 / 30 m seamless nationwide; 1 / 3 m from project lidar where +# it covers the AOI; 60 m is Alaska-only. A resolution with no data for the AOI +# logs + raises DEMResolutionUnavailable (default stays 10 m). +def test_get_dem(): + fimbox.DEMProcessor( + boundary=test_boundary, + output_dir=OUT_DIR, + resolution=10, # 1, 3, 10 (default), 30, 60 + out_name="dem.tif", # default is 3dep_dem_m.tif + # epsg=None, # output CRS; None -> auto UTM zone + # layer=None, # layer name if boundary has multiple + # use_dask=True, # dask chunking for reproject/heal + # chunksize=None, # None -> auto from CPU count; or set px + ) + + +# "give me 1 m, else just 10 m": fallback_to_10m downgrades to 10 m (and logs) +# when the requested resolution isn't available for the AOI. +# def test_get_dem_fallback(): +# fimbox.DEMProcessor( +# boundary=test_boundary, +# output_dir=OUT_DIR, +# resolution=1, # 1 m where lidar exists, else fall back +# out_name="dem.tif", +# fallback_to_10m=True, # default False -> raises if unavailable +# ) + +# Bring your own DEM: pass dem_file and it gets the SAME conditioning as a +# downloaded one (reproject -> hole-fill -> clip to boundary). +# def test_process_byo_dem(): +# fimbox.DEMProcessor( +# boundary=test_boundary, +# output_dir=OUT_DIR, +# out_name="dem.tif", +# resolution=10, +# dem_file="path/to/my_dem.tif", +# ) + +# def test_get_osm_roads(): +# fimbox.DownloadOSMRoads().download( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# out_name="osm_roads.gpkg", +# out_layer="osm_roads", +# ) + +# def test_get_osm_bridges(): +# fimbox.DownloadOSMBridges().download( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# out_name="osm_bridges.gpkg", +# out_layer="osm_bridges", +# ) + + +# # Get the HUC8 information +# def test_get_huc8_info(): +# huc8_info = fimbox.getHUC8Info( +# boundary=test_boundary, +# calc_overlap=True, +# save=True, +# out_dir=OUT_DIR, +# ) +# logging.getLogger(__name__).info(f"HUC8 info:\n{huc8_info}") + + +# USGS gauge points — downloads from the ArcGIS Online FeatureServer. +# Uncomment to run live; the smoke test below always runs. + +# def test_download_usgs_gages(): +# """Download USGS gauges inside the test boundary into ../out/usgs_gages.gpkg.""" +# gdf = fimbox.DownloadUSGSGages().download( +# boundary=test_boundary, +# aoi_id="08060202", +# out_dir=OUT_DIR, +# out_name="usgs_gages.gpkg", +# out_layer="usgs_gages", +# ) +# log = logging.getLogger(__name__) +# log.info(f"USGS gauges downloaded: {len(gdf)} features") +# assert {"location_id", "feature_id", "aoi_id", "source"}.issubset(gdf.columns) + + +# def test_usgs_gages_signature(): +# """Smoke test: DownloadUSGSGages is exported and has the documented API.""" +# import inspect + +# assert hasattr(fimbox, "DownloadUSGSGages") +# sig = inspect.signature(fimbox.DownloadUSGSGages.download) +# expected = {"boundary", "aoi_id", "out_dir", "out_name", "out_layer"} +# assert expected.issubset(sig.parameters.keys()), ( +# f"DownloadUSGSGages.download missing kwargs: {expected - set(sig.parameters)}" +# ) diff --git a/tests/test_fimevaluation.py b/tests/test_fimevaluation.py new file mode 100644 index 0000000..7f899ba --- /dev/null +++ b/tests/test_fimevaluation.py @@ -0,0 +1,123 @@ +""" +Author: Supath Dhital +Date Created: July 2026 + +FIM evaluation tests: query benchmark FIMs from the FIMbench database +(queryBenchmarkFIM) and evaluate candidate FIMs against them with FIMeval +(evaluateFIM). + +Point AOI_DIR at a working directory that already has flood maps in +fim-outputs/ (produced by test_fimgeneration). fimbench and fimeval install +together with fimbox. +""" + +from __future__ import annotations + +from pathlib import Path + +AOI_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" + +# Boundary-extraction method: "smallest_extent" | "convex_hull" | "AOI" +METHOD_NAME = "smallest_extent" + +# Optional query filters (edit to match your AOI / event). +EVENT_DATE = None # e.g. "2017-08-30" +START, END = "2016-04-01", "2026-01-01" +TIER = None # e.g. "HWM", "tier1" + + +# COMBINED — the whole evaluation pipeline in one go: query the FIMbench +# catalog with the AOI's newest flood extent, download the matched benchmark +# assets into /benchmark-data/, then run FIMeval (EvaluateFIM + +# contingency maps + metric plots) on the staged case. +def test_fimevaluation_combined(): + from fimbox import evaluateFIM, queryBenchmarkFIM + + query = queryBenchmarkFIM( + AOI_DIR, # footprint = newest extent raster in /fim-outputs/ + # raster_path="my_fim.tif", # explicit candidate raster instead + # boundary_path="my_aoi.gpkg", # or an AOI boundary vector + # huc8="03020201", # narrow by basin + tier=TIER, # narrow by benchmark tier; None -> all tiers + event_date=EVENT_DATE, # exact event date; None -> ignore + start_date=START, # date-range start + end_date=END, # date-range end + area=True, # add overlap % / km^2 per match + download=True, # fetch GeoTIFF + GeoPackage + # out_dir="downloads/", # default /benchmark-data/ + ) + print(query) # pretty match summary + + result = evaluateFIM( + AOI_DIR, + # candidate=None, # default: every extent .tif in fim-outputs/ + # benchmark=None, # default: newest .tif in benchmark-data/ + # case_name=None, # default: first candidate's stem + method_name=METHOD_NAME, + # aoi_boundary="my_aoi.gpkg", # required when method_name="AOI" + # pwb_dir="my_pwb.gpkg", # own permanent-water-bodies vector + # target_crs="EPSG:32633", # outside CONUS (default EPSG:5070) + # target_resolution=10, # m, when resolutions differ + contingency_map=True, + plot_metrics=True, + # building_footprint=True, # building-level agreement (GEE auth) + ) + assert result.case_dir.is_dir() + assert result.output_dir.is_dir() + assert "benchmark" in result.benchmark.name.lower() + assert result.metrics_files, "FIMeval produced no metrics CSV" + print(result.metrics) + + +# # STEP BY STEP — each stage on its own. + +# def test_step_query_catalog_only(): +# """Catalog-only discovery: what benchmarks exist in a date range +# (no AOI, no download). Returns plain dicts straight from the catalog.""" +# from fimbox import queryBenchmarkFIM + +# response = queryBenchmarkFIM(start_date=START, end_date=END) +# print(response) +# for record in response.records: +# print(record) + + +# def test_step_query_by_filename(): +# """Direct download by exact catalog filename.""" +# from fimbox import queryBenchmarkFIM + +# response = queryBenchmarkFIM( +# file_name="HWM_10_0m_20160928_20161009_780051W352232N_BM.tif", +# download=True, +# out_dir=AOI_DIR / "benchmark-data", +# ) +# assert response.downloads + + +# def test_step_evaluate_own_benchmark(): +# """Evaluate against a benchmark raster you already have on disk +# (skips the FIMbench query entirely).""" +# from fimbox import evaluateFIM + +# result = evaluateFIM( +# AOI_DIR, +# benchmark="path/to/my_benchmark.tif", +# case_name="own_benchmark", +# method_name="convex_hull", +# ) +# assert result.metrics_files + + +# def test_step_building_footprint(): +# """Building-level agreement analysis. Uses the Microsoft global building +# footprints via Google Earth Engine by default (pops a GEE auth prompt), +# or pass building_footprint_file= to use your own vector.""" +# from fimbox import evaluateFIM + +# result = evaluateFIM( +# AOI_DIR, +# method_name=METHOD_NAME, +# building_footprint=True, +# # building_footprint_file="path/to/footprints.gpkg", +# ) +# assert result.output_dir.is_dir() diff --git a/tests/test_fimgeneration.py b/tests/test_fimgeneration.py new file mode 100644 index 0000000..c13d1a2 --- /dev/null +++ b/tests/test_fimgeneration.py @@ -0,0 +1,102 @@ +""" +Author: Supath Dhital +Date Updated: June 2026 + +FIM generation driven by the CSVs in /discharge-inputs/. + +generateFIM(aoi_dir).from_discharge_inputs(...) selects which discharge CSVs to +run and generates an inundation extent raster for each (named after the input +CSV) in /fim-outputs/. Pass depth=True to also write a depth raster. + +Selection modes: + * nothing -> every CSV in discharge-inputs/ + * csv= -> that one CSV + * date="YYYY-MM-DD" -> CSVs whose filename carries that day/instant stamp + * start=.., end=.. -> CSVs whose YYYYMMDD token falls in the range + +Point AOI_DIR at a working directory that already has branches and at least one +discharge CSV (produced by the streamflow pipeline, e.g. getNWMretrospective). +""" + +from __future__ import annotations + +from pathlib import Path + +import pytest + +from fimbox import generateFIM + +AOI_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" +N_WORKERS = 4 + +# Optional selection filters (edit to match the CSVs you have). +EVENT_DATE = "2020-05-20 12:00:00" +START = "2020-05-19" +END = "2020-05-22" + +_BRANCHES_DIR = ( + AOI_DIR / "watershed-data" / "branches" + if (AOI_DIR / "watershed-data" / "branches").is_dir() + else AOI_DIR / "branches" +) +_skip_no_branches = pytest.mark.skipif( + not _BRANCHES_DIR.is_dir(), reason=f"AOI not present: {_BRANCHES_DIR}" +) + + +# @_skip_no_branches +# def test_extract_feature_ids(): +# out_csv = extract_feature_ids(AOI_DIR) +# assert out_csv.is_file() +# print(f"\nfeature_id.csv -> {out_csv}") + + +# default: generate FIM for every discharge CSV in the AOI +@_skip_no_branches +def test_fim_all_discharge_inputs(): + results = generateFIM( + AOI_DIR, n_workers=N_WORKERS, depth=True + ).from_discharge_inputs() + assert results + for r in results: + print(f" extent={r.extent_path}") + assert r.extent_path is not None and Path(r.extent_path).is_file() + + +# # a specific CSV +# @_skip_no_branches +# def test_fim_specific_csv(): +# csvs = sorted((AOI_DIR / "discharge-inputs").glob("*.csv")) +# if not csvs: +# pytest.skip("no discharge CSVs to pick from") +# results = generateFIM(AOI_DIR, n_workers=N_WORKERS).from_discharge_inputs(csv=csvs[0]) +# assert len(results) == 1 + + +# # match by date stamp in the filename +# @_skip_no_branches +# def test_fim_by_date(): +# results = generateFIM(AOI_DIR, n_workers=N_WORKERS).from_discharge_inputs( +# date=EVENT_DATE +# ) +# assert results + + +# # match by date range +# @_skip_no_branches +# def test_fim_by_range(): +# results = generateFIM(AOI_DIR, n_workers=N_WORKERS).from_discharge_inputs( +# start=START, end=END +# ) +# assert results + + +# # also write the depth raster +# @_skip_no_branches +# def test_fim_with_depth(): +# results = generateFIM( +# AOI_DIR, n_workers=N_WORKERS, depth=True +# ).from_discharge_inputs(date=EVENT_DATE) +# assert results +# r = results[0] +# assert r.depth_path is not None and Path(r.depth_path).is_file() diff --git a/tests/test_generate_dem_diff.py b/tests/test_generate_dem_diff.py new file mode 100644 index 0000000..6f30194 --- /dev/null +++ b/tests/test_generate_dem_diff.py @@ -0,0 +1,71 @@ +""" +Tests for bridge DEM processing pipeline. +Step 1 (generateBridgeRaster): streams USGS LiDAR and writes per-bridge .tif +Step 2 (BridgeDEMDiff): computes lidar_elev - dem_elev and saves bridge_elev_diff.tif +""" + +import logging +from pathlib import Path + +import fimbox + +log = logging.getLogger(__name__) + +OUT_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" / "watershed-data" +bridge_gpkg = OUT_DIR / "osm_bridges_subset.gpkg" +dem_path = OUT_DIR / "dem.tif" +out_dir = OUT_DIR + + +# check which bridges already have rasters vs still pending (safe to run anytime) +def test_bridge_raster_status(): + info = fimbox.generateBridgeRaster( + bridge_gpkg=bridge_gpkg, + out_dir=out_dir, + ).status() + assert "total" in info + + +# download LiDAR and build per-bridge elevation tifs +def test_generate_bridge_raster(): + tif_dir = fimbox.generateBridgeRaster( + bridge_gpkg=bridge_gpkg, + out_dir=out_dir, + resolution=10.0, + buffer_m=10.0, + n_workers=4, + # id_col="my_id", # only needed if gpkg has no 'osmid' column + ).run() + log.info(f"Per-bridge tifs --> {tif_dir}") + + +# compute difference raster +def test_bridge_dem_diff(): + out_path = fimbox.BridgeDEMDiff( + dem_path=dem_path, + lidar_tif_dir=OUT_DIR / "bridge_dem" / "lidar_osm_rasters", + bridge_gpkg=bridge_gpkg, + out_dir=out_dir, + out_name="bridge_elev_diff.tif", + n_workers=4, + ).run() + log.info(f"Bridge diff raster --> {out_path}") + + +# Run both steps end-to-end +# def test_full_pipeline(): +# tif_dir = fimbox.generateBridgeRaster( +# bridge_gpkg=bridge_gpkg, +# out_dir=out_dir, +# resolution=10.0, +# n_workers=4, +# ).run() +# +# out_path = fimbox.BridgeDEMDiff( +# dem_path=dem_path, +# lidar_tif_dir=tif_dir, +# bridge_gpkg=bridge_gpkg, +# out_dir=out_dir, +# n_workers=4, +# ).run() +# print(f"Done: {out_path}") diff --git a/tests/test_getallinputdata.py b/tests/test_getallinputdata.py new file mode 100644 index 0000000..4989dba --- /dev/null +++ b/tests/test_getallinputdata.py @@ -0,0 +1,94 @@ +# Example Usage: +from pathlib import Path + +import fimbox + +PKG_ROOT = Path(__file__).resolve().parents[1] +REPO_ROOT = Path(__file__).resolve().parents[2] + +test_boundary = PKG_ROOT / "docs" / "test_boundary" / "test_smallB.shp" +OUT_DIR = REPO_ROOT / "out" +test_huc8 = "08060202" # Yazoo River basin, MS + + +# Combined preprocessing pipeline tests +# Run full pipeline from a boundary shapefile +def test_preprocess_all_from_boundary(): + pp = fimbox.getAllInputData( + boundary=test_boundary, + out_dir=OUT_DIR, + buffer_m=5000, # metres to buffer boundary for data downloads + headwater_buffer_cells=8, # pixels to shrink buffer for headwater clip + get_flowlines=True, # set False to use your own flowlines and corresponding catchments + get_catchments=True, # set False to skip NWM catchments--> use + resolution="medium", # "high" -> NHDPlus HR flowlines/catchments via pynhd; "medium" (default) -> NWM. Lakes always NWM. + identifier="nwmmr", # filename prefix for ALL source files; flows download->processing. Default "nwm". + ) + pp.run() + + +# Bring your own flowlines / catchments / DEM (any column names, any source). +# Pass the file paths + field maps; flowlines/catchments are normalised to the +# pipeline schema (streams: ID, order_, levpa_id, feature_id[=ID]; catchments: ID), +# the DEM is reprojected/clipped/hole-filled, and all files are saved under the +# chosen identifier prefix so the whole pipeline picks them up automatically. +# def test_preprocess_byo_inputs(): +# pp = fimbox.getAllInputData( +# boundary=test_boundary, +# out_dir=OUT_DIR, +# flowlines="path/to/my_flowlines.gpkg", +# catchments="path/to/my_catchments.gpkg", +# stream_fields={"ID": "nhdplusid", "order_": "streamorde", "levpa_id": "levelpathi"}, +# catchment_fields={"ID": "nhdplusid"}, # must match the flowline reach id +# dem="path/to/my_dem.tif", # reprojected, clipped, and hole-filled like a downloaded DEM +# identifier="3dhp", # files saved as 3dhp_subset_streams.gpkg etc.; whole pipeline follows it +# ) +# pp.run() + + +# # Run full pipeline from a HUC8 ID +# # get_flowlines / get_catchments default to True (downloads everything, +# # including OSM bridges). Set either to False to skip that dataset and use +# # your own instead. +# def test_preprocess_all_from_huc8(): +# pp = fimbox.getAllInputData( +# huc8=test_huc8, +# out_dir=OUT_DIR, +# buffer_m=2000, +# headwater_buffer_cells=8, +# get_flowlines=True, # set False to use your own flowlines and corresponding catchments +# get_catchments=True, # set False to skip NWM catchments--> use your own in later steps +# ) +# pp.run() + + +# Same pipeline, but bring your own flowlines/catchments +# (skips the NWM flowline + catchment downloads; everything else still runs) +# def test_preprocess_all_byo_flowlines_catchments(): +# pp = fimbox.getAllInputData( +# huc8=test_huc8, +# out_dir=OUT_DIR, +# buffer_m=2000, +# headwater_buffer_cells=8, +# get_flowlines=False, +# get_catchments=False, +# ) +# pp.run() + + +# Run individual steps +# def test_preprocess_dem_only(): +# pp = fimbox.getAllInputData(boundary=test_boundary, out_dir=OUT_DIR) +# pp.run_dem() + +# def test_preprocess_nhd_only(): +# pp = fimbox.getAllInputData(boundary=test_boundary, out_dir=OUT_DIR) +# pp.run_nhd() + +# def test_preprocess_nld_only(): +# pp = fimbox.getAllInputData(boundary=test_boundary, out_dir=OUT_DIR) +# pp.run_nld() + +# def test_preprocess_osm_only(): +# pp = fimbox.getAllInputData(boundary=test_boundary, out_dir=OUT_DIR) +# pp.run_osm() diff --git a/tests/test_nwmstreamflow.py b/tests/test_nwmstreamflow.py new file mode 100644 index 0000000..4169e18 --- /dev/null +++ b/tests/test_nwmstreamflow.py @@ -0,0 +1,87 @@ +""" +Author: Supath Dhital +Date Created: June 2026 + +Streamflow retrieval / plot / statistics — minimal, call-the-function tests. + +Point AOI_DIR at a working directory whose feature_id.csv exists (or pass a +feature_ids list / CSV per call). Edit the dates and USGS site to your basin, +then run the functions you want. +""" + +from __future__ import annotations + +from pathlib import Path + +from fimbox import ( + getNWMretrospective, +) + +AOI_DIR = Path(__file__).resolve().parents[2] / "out" / "test_smallB" + +START = "2016-10-05" +END = "2016-10-20" +EVENT = "2020-10-10 21:00:00" + + +# retrospective — different extraction combinations +def test_retrospective_event_date(): + getNWMretrospective(AOI_DIR, date=EVENT) + + +# def test_retrospective_range_continuous(): +# # start + end, nothing else -> one CSV per hour +# getNWMretrospective(AOI_DIR, start=START, end=END) + + +# def test_retrospective_range_sortby(): +# # start + end + sortby -> one aggregated CSV +# getNWMretrospective(AOI_DIR, start=START, end=END, sortby="maximum") + + +# def test_retrospective_feature_ids_list(): +# # pass feature_ids directly instead of relying on the AOI's feature_id.csv +# getNWMretrospective(AOI_DIR, feature_ids=[FEATURE_ID], date=EVENT) + + +# # forecast — different combinations +# def test_forecast_shortrange(): +# getNWMforecast(AOI_DIR, "shortrange") + + +# def test_forecast_mediumrange_maxsort(): +# getNWMforecast(AOI_DIR, "mediumrange", sort_by="maximum") + + +# def test_forecast_specific_cycle(): +# getNWMforecast(AOI_DIR, "shortrange", forecast_date="2024-06-01", hour=12) + + +# # USGS observations +# def test_usgs_fetch(): +# USGSData(AOI_DIR).fetch([USGS_SITE], START, END) + + +# def test_usgs_feature_id_pairs(): +# # which USGS gage falls on which reach (feature_id) within the AOI +# pairs = get_usgs_fid_pairs(AOI_DIR) +# print(pairs) + + +# # plots +# def test_plot_feature_id(): +# plot_nwm(AOI_DIR, [FEATURE_ID], START, END) + + +# def test_plot_usgs(): +# plot_usgs(AOI_DIR, [USGS_SITE], START, END) + + +# def test_plot_usgs_and_feature_id(): +# # time series overlay of USGS and the NWM feature_id together +# plot_comparison(AOI_DIR, FEATURE_ID, USGS_SITE, START, END) + + +# # statistics +# def test_statistics_usgs_vs_nwm(): +# calculate_statistics(AOI_DIR, FEATURE_ID, USGS_SITE, START, END) diff --git a/tests/test_preprocessDEM.py b/tests/test_preprocessDEM.py new file mode 100644 index 0000000..3e7c636 --- /dev/null +++ b/tests/test_preprocessDEM.py @@ -0,0 +1,28 @@ +# Example Usage: +import logging +from pathlib import Path + +import fimbox + +log = logging.getLogger(__name__) + +PKG_ROOT = Path(__file__).resolve().parents[1] +REPO_ROOT = Path(__file__).resolve().parents[2] + +boundary = PKG_ROOT / "docs" / "test_boundary" / "test_smallB.shp" +OUT_DIR = REPO_ROOT / "out" + + +def test_process_dem(): + output_path = fimbox.DEMProcessor( + boundary=boundary, + resolution=10, # 3DEP resolution in m: 1, 3, 10 (default), 30, 60 + output_dir=OUT_DIR / "dem_test", + # layer=None, # if boundary is a geopackage with multiple layers + # dem_file=None, # local DEM to condition instead of fetching + # epsg=None, # output CRS EPSG; None auto-detects the UTM zone + # fallback_to_10m=False, # if resolution unavailable, use 10m not raise + # use_dask=True, # dask chunking for the reproject/heal stage + # chunksize=None, # dask chunk edge in px; None -> auto from CPU count + ).result_path + log.info(f"3DEP DEM --> {output_path}") diff --git a/tests/test_preprocessing_hucs.py b/tests/test_preprocessing_hucs.py new file mode 100644 index 0000000..31f9aea --- /dev/null +++ b/tests/test_preprocessing_hucs.py @@ -0,0 +1,19 @@ +# importing the fimbox preprocessing module to test HUCChecker +import logging + +import fimbox + +log = logging.getLogger(__name__) +checker = fimbox.HUCChecker() + + +def test_huc_checker(): + # Single HUC Query + r = checker.check_any("03020202", strict=False) + log.info(f"total={r.n_total} found={r.n_found} missing={r.n_missing}") + log.info(f"missing: {r.missing_hucs}") + + # List of HUCs Query + r = checker.check_any(["01010001", "99999999"], strict=False) + log.info(f"total={r.n_total} found={r.n_found} missing={r.n_missing}") + log.info(f"missing: {r.missing_hucs}")