From c841d64780f805f25a14133ede096bcc17932ea6 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 20 Jul 2026 13:48:51 +0200 Subject: [PATCH 1/5] Add a pytest suite covering the nhflotools code 09pwnmodel2 uses The boundary conditions and the layer model of the active PWN model were untested: of the thirteen entry points 01_pwnmodel2.py imports, only interface_elevation and the Tata wells had any coverage. This adds 126 tests over that whole closure, running in about four seconds. Tests build their own DISV vertex grid (tests/util.py) instead of calling gridgen, so the unit tests need no executables, no network and no large datasets, and every cell area is an exact float that makes first-principles expected values possible. Live services -- HHNK, REGIS, RWS, background tiles -- are monkeypatched by name; the rest of nlmod runs for real, so the suite doubles as a compatibility canary for nlmod@dev. One mf6-marked test runs MODFLOW on 3x3x2 cells and checks budget closure plus an analytic head field, fetching the executable through nlmod.util.get_exe_path on first use. Every file was mutation-checked while being written: a plausible bug was introduced into the source, the test confirmed to fail, and the source restored. That covers the panden conductance formula reverting to resistance aggregation (a silent ~1e4x error), the polder fallback conductance being zeroed, the northsea guard being dropped, numeric_only being removed from the well median, extrapolate_ds's in-place contract breaking, the nearest-donor metric swapping axes, the budget threshold inverting, and the head field mirroring on the x axis. pyarrow becomes a runtime dependency: well.py reads a feather file via pandas, which needs it in production as well as in the tests. Also adds a test workflow keyed on the current nlmod@dev and nhflodata@main commits so the moving branches cannot serve a stale environment, a README overview separating the modules 09pwnmodel2 uses from the untested rest, and TEST_PLAN.md recording the design, the mutation evidence and the deviations. Ten defects found while writing these tests are filed as #60-#69 rather than pinned by a test, so none of them is entrenched as expected behaviour. The REGIS-completeness guard (#65) is covered by an xfail(strict=True) test that turns into a visible failure once the guard is fixed. Removes tests/test_hhnk.py and tests/test_nhflo_utils.py, which held no tests while suggesting hhnk.py and nhflo_utils.py were covered. --- .github/workflows/test.yml | 53 +++ README.md | 59 ++- TEST_PLAN.md | 674 +++++++++++++++++++++++++++++ pyproject.toml | 11 +- tests/conftest.py | 41 ++ tests/test_hhnk.py | 1 - tests/test_major_surface_waters.py | 226 ++++++++++ tests/test_mf6_smoke.py | 114 +++++ tests/test_nhflo_utils.py | 0 tests/test_nhflodata_contract.py | 113 +++++ tests/test_nhi_chloride.py | 178 ++++++++ tests/test_panden.py | 184 ++++++++ tests/test_polder.py | 136 ++++++ tests/test_postprocessing.py | 238 ++++++++-- tests/test_pwnlayers3_layers.py | 586 +++++++++++++++++++++++++ tests/test_pwnlayers3_plot.py | 182 ++++++++ tests/test_pwnlayers_get_top.py | 152 +++++++ tests/test_well.py | 129 ++++++ tests/util.py | 203 +++++++++ 19 files changed, 3246 insertions(+), 34 deletions(-) create mode 100644 .github/workflows/test.yml create mode 100644 TEST_PLAN.md create mode 100644 tests/conftest.py delete mode 100644 tests/test_hhnk.py create mode 100644 tests/test_major_surface_waters.py create mode 100644 tests/test_mf6_smoke.py delete mode 100644 tests/test_nhflo_utils.py create mode 100644 tests/test_nhflodata_contract.py create mode 100644 tests/test_nhi_chloride.py create mode 100644 tests/test_polder.py create mode 100644 tests/test_pwnlayers3_layers.py create mode 100644 tests/test_pwnlayers3_plot.py create mode 100644 tests/test_pwnlayers_get_top.py create mode 100644 tests/util.py diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml new file mode 100644 index 0000000..99a4de5 --- /dev/null +++ b/.github/workflows/test.yml @@ -0,0 +1,53 @@ +name: Test python files + +on: + push: + branches: + - main + pull_request: + branches: + - main + +jobs: + test: + name: Run pytest + runs-on: ubuntu-latest + concurrency: + group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.sha }} + cancel-in-progress: true + env: + HATCH_VERBOSE: 1 + steps: + - uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version-file: "pyproject.toml" + + - name: Install uv + uses: astral-sh/setup-uv@v3 + + - name: Install hatch + run: uv tool install hatch + + # nlmod and nhflodata track moving branches, so key the cache on their current + # head commits: a new upstream commit rebuilds the env, otherwise it is reused. + - name: Resolve moving dependencies + id: deps + run: | + echo "data_sha=$(git ls-remote https://github.com/NHFLO/data.git main | cut -f1)" >> "$GITHUB_OUTPUT" + echo "nlmod_sha=$(git ls-remote https://github.com/gwmod/nlmod.git dev | cut -f1)" >> "$GITHUB_OUTPUT" + + - name: Cache the test environment + uses: actions/cache@v4 + with: + path: ~/.cache/uv + key: uv-${{ hashFiles('pyproject.toml') }}-${{ steps.deps.outputs.data_sha }}-${{ steps.deps.outputs.nlmod_sha }} + + # The MODFLOW binaries are not cached: nlmod downloads them into its own package + # directory, which lives inside the hatch environment and is rebuilt whenever the + # environment cache misses. The single mf6-marked test fetches them via + # nlmod.util.get_exe_path on first use (~30 s) and reuses them afterwards. + - name: Run tests + run: hatch run test:test -- -m "not network" --durations=15 diff --git a/README.md b/README.md index f4658f5..26355ca 100644 --- a/README.md +++ b/README.md @@ -4,4 +4,61 @@ Deze repository bevat de `nhflotools` Python package met tools om NHFLO modellen Note that the content of this repository is available under a not-so-permissive open-source licence: GNU AGPLv3. Please have a look at [choose a license](https://choosealicense.com/licenses/agpl-3.0/) for the key conditions and limitations of this license before getting started. ## Installatie -De `nhflotools` package kan worden geinstallerd door deze repository te clonen en vervolgens `pip install -e . --config-settings editable_mode=strict` te runnen vanuit de package map. \ No newline at end of file +De `nhflotools` package kan worden geinstallerd door deze repository te clonen en vervolgens `pip install -e . --config-settings editable_mode=strict` te runnen vanuit de package map. + +## Overzicht van de modules + +Het modelscript [`09pwnmodel2`](https://github.com/NHFLO/models/tree/main/modelscripts/09pwnmodel2) is +het actieve PWN-model en gebruikt maar een deel van deze package. Die modules vormen de +onderhouden kern: ze worden in CI met pytest getest. De overige modules zijn niet in gebruik +door 09pwnmodel2 en zijn ongetest — behandel ze als legacy of werk-in-uitvoering. + +### Gebruikt door 09pwnmodel2 (getest in CI) + +| Module | Wat het doet | Gebruikt in 09pwnmodel2 | +|---|---|---| +| `major_surface_waters.py` | Grote RWS-oppervlaktewateren: Noordzee als CHD (met zeespiegel-tijdreeks en 18.000 mg/l chloride), IJsselmeer/Markermeer/Noordzeekanaal als GHB | `get_chd_ghb_data_from_major_surface_waters`, `chd_ghb_from_major_surface_waters` | +| `nhi_chloride.py` | NHI-chlorideconcentratie interpoleren naar het modelgrid als begintoestand voor transport | `get_nhi_chloride_concentration` | +| `panden.py` | Infiltratiepanden (ICAS/IKIEF) als RIV-package | `riv_from_oppervlakte_pwn` | +| `polder.py` | Polderpeilgebieden van HHNK als DRN-package, met maaiveld als terugval | `drn_from_waterboard_data` | +| `postprocessing.py` | Waterbalanscontrole, modeluitvoer inlezen, grensvlakken zoet/brak en resultaatkaarten | `check_budget_discrepancy`, `add_output_to_ds`, `plot_result_maps` | +| `pwnlayers/layers.py` | Alleen de maaiveldhoogte uit AHN, inclusief opvullen bij oppervlaktewater en zee | `get_top_from_ahn` | +| `pwnlayers3/layers.py` | Het PWN-lagenmodel (v3): botm, kh en kv uit de bodemlagenkartering, samengevoegd met REGIS | `get_pwn_layer_model` | +| `pwnlayers3/plot.py` | Diagnostische dwarsdoorsneden van het lagenmodel | `plot_diagnostic_cross_sections` | +| `well.py` | Winnings- en infiltratieputten van PWN en Tata Steel | `get_wells_pwn_dataframe`, `get_wells_tata_dataframes` | +| `pwnlayers/merge_layer_models.py` | Twee lagenmodellen samenvoegen via een koppeltabel, met overgangszone | indirect, via `pwnlayers3` | +| `pwnlayers/utils.py` | Ontbrekende en elkaar kruisende laagbodems repareren | indirect, via `merge_layer_models` | + +### Niet gebruikt door 09pwnmodel2 (ongetest) + +| Module | Wat het doet | +|---|---| +| `bergen_utils.py` | Lagenmodel en oppervlaktewater voor het oudere Bergen-model | +| `berging_utils.py` | Bergingscoëfficiënten per perceel, afgeleid van nabijgelegen oppervlaktewater | +| `bofek.py` | BOFEK-bodemprofielen en de bijbehorende berging | +| `cropfactor.py` | Gewasfactoren toepassen op verdamping | +| `geoconverter/` | Command-line tool om geodata naar het NHFLO-dataformaat te converteren | +| `hhnk.py` | Peilbuisreeksen ophalen uit de FEWS-webservice van HHNK | +| `nhflo_utils.py` | Verzameling oudere plot- en gridhulpfuncties | +| `pwnlayers/io.py` | Inlezen van de Mensink- en Bergen-bodemlagen (lagenmodel v1) | +| `utils.py` | Lokale kopie van nlmod's MODFLOW-binaries-afhandeling | + +## Tests + +De suite telt 126 tests en draait in ongeveer 4 seconden. De tests gebruiken uitsluitend +kleine synthetische modellen — geen netwerk, geen gridgen en geen grote datasets — zodat +ze op elke pull request in CI meedraaien. Alle live webservices (HHNK, REGIS, RWS, +achtergrondkaarten) worden per naam gemonkeypatcht; de rest van nlmod draait echt, waardoor +de suite meteen dienstdoet als compatibiliteitscanary voor `nlmod@dev`. + +De enige uitzondering is één test met de marker `mf6`: die draait MODFLOW echt op een +model van 3x3x2 cellen en controleert de waterbalans plus een analytische oplossing. De +binaries worden zo nodig automatisch door nlmod gedownload en daarna hergebruikt. + +```bash +hatch run test:test # de hele suite +pytest -m "not mf6" # zonder de MODFLOW-run +``` + +De opzet en onderbouwing van de suite staan in [`TEST_PLAN.md`](TEST_PLAN.md), inclusief +de mutatietests waarmee per bestand is aangetoond dat de tests echte fouten vangen. \ No newline at end of file diff --git a/TEST_PLAN.md b/TEST_PLAN.md new file mode 100644 index 0000000..4023d53 --- /dev/null +++ b/TEST_PLAN.md @@ -0,0 +1,674 @@ +# Test plan — nhflotools, scoped to the 09pwnmodel2 closure + +Implementation plan only; nothing here is implemented yet. Scope = the 13 nhflotools entry +points `models/modelscripts/09pwnmodel2/01_pwnmodel2.py` imports, plus their transitive +nhflotools closure (`pwnlayers.merge_layer_models`, `pwnlayers.utils`, +`panden.get_oppervlakte_pwn_shapes`). Target: a lean, meaningful pytest suite that runs +green on a fresh ubuntu GitHub runner with zero live network at test time. + +Provenance: produced by a multi-agent review (9 module deep-reads, call-site trace, +git-history mining, existing-test audit, nlmod fixture survey, data/CI feasibility, two +competing designs, adversarial critique). All file:line claims below were verified against +HEAD during that review. + +--- + +## 1. Current state (verified baseline) + +- Existing suite: 18 tests, all passing on a healthy Python 3.14 env in <2 s + (`test_postprocessing.py` 9: `interface_elevation` + one mocked + `check_budget_discrepancy`; `test_well.py` 9: `get_wells_tata_dataframes` only). +- `test_panden.py` and `test_nhflo_utils.py` are 0 bytes; `test_hhnk.py` is a lone + docstring — they collect nothing and misrepresent coverage. +- No `conftest.py` (a previous one plus `test_polder.py`/`test_lakes.py`/`test_pwnlayers.py` + were deleted in b86583e; `__pycache__` remnants confirm). +- CI runs lint only (`.github/workflows/lint.yml`); there is no test workflow. +- `[tool.pytest.ini_options]` has only `testpaths`; pytest-env is installed but unused. +- The hatch-registered test env (`.direnv/test`) is a stale Python 3.11 env that can no + longer resolve deps (`requires-python >=3.14`); it must be removed and recreated once. +- Both NHFLO/tools and NHFLO/data are public → CI needs no secrets. With + `NHFLODATA_LOCATION` unset, `nhflodata.get_abs_data_path` resolves to mockup data + packaged in the wheel; all datasets the script requests have mockup variants. +- **Blocker found: `pyarrow` is missing from dependencies** — `well.py:46` + (`pd.read_feather`) needs it in production and in tests; `import pyarrow` fails in the + tools venv. Add it to the **runtime** deps (not just test extra). + +## 2. Untested-but-used surface + +Zero coverage today for: `major_surface_waters` (both functions), `nhi_chloride`, +`panden` (both), `polder`, `well.get_wells_pwn_dataframe`, +`postprocessing.add_output_to_ds` / `plot_result_maps`, `pwnlayers.layers.get_top_from_ahn`, +all of `pwnlayers3` (layer model + plot). This includes the WEL/RIV/DRN/GHB/CHD boundary +conditions and the layer model — the physics core of 09pwnmodel2. + +--- + +## 3. Architecture + +### 3.1 `tests/util.py` (vendored helpers, ~90 lines) + +nlmod's `tests/util.py` is not importable from an installed nlmod, so vendor equivalents: + +| Helper | Builds | +|---|---| +| `make_rect_vertex_ds(nx=2, ny=2, delr=100.0, botms=(-10.0, -20.0), top=0.0, kh=5.0, transport=0)` | Hand-built rectangular pseudo-vertex ds: dims `(layer, icell2d)`, coords cell-center `x/y`, grid geometry `xv/yv/icvert`, vars `top/botm/kh/kv/area/idomain`, attrs `gridtype='vertex'`, `extent`, `model_name='test'`, `transport`. **No gridgen, no binaries, no network** — the single biggest feasibility lever. | +| `make_structured_ds(...)` | Thin wrapper over `nlmod.get_ds(extent, delr=100, ...)` with `download_exe=False` (param verified present in nlmod/dims/base.py:548). | +| `add_time(ds)` | `nlmod.time.set_ds_time(ds, time=[1], start='2022-01-01', steady=True)` — required before sim/tdis. | +| `make_gwf_disv(ds, tmp_path)` | `nlmod.sim.sim` + `tdis` + `nlmod.gwf.gwf` + `disv`; package construction only, never run. | +| `make_rws_gdf(rows)` | Synthetic EPSG:28992 GeoDataFrame factory, columns `OWMNAAM/peil/bweerstand/geometry`. | +| `write_mf6_listing(path, pct_disc, budgetkey)` | Emits a real MF6 listing **derived from the committed `.lst` of the smoke test's run** (§5.11), editing only the PERCENT DISCREPANCY numbers. Hand-rolled listing text risks `Mf6ListBudget.get_dataframes()` returning None → false-positive "Could not parse" tests. | + +### 3.2 `tests/conftest.py` + +| Fixture | Scope | Content | +|---|---|---| +| `_hygiene` (autouse) | function | Vendored verbatim from nlmod tests/conftest.py:39-46 (`plt.close('all')`, `gc.collect()`, clear xarray FILE_CACHE) + `monkeypatch.delenv('NHFLODATA_LOCATION', raising=False)` so a developer's data mount can never redirect mockup-resolution tests. | +| `vertex_ds` | function | Fresh `make_rect_vertex_ds()` per test — several targets mutate ds in place (`northsea`, `sfw_*`, `drn_*`, `thickness`); no shared mutable state. Build cost sub-millisecond. | +| `gwf_disv` | function | `(ds, sim, gwf)` in `tmp_path`. | +| `pwn_data_tree` | session | Synthetic bodemlagen tree: 7 boundary GeoJSONs (squares buffered +1e-3 m past cell edges to dodge the `min_area_fraction=1.0` float-equality trap), `botm/botm.geojson` point layer (all 14 W/S columns), `conductances/K*/KD*/C*` GeoJSONs, `triwaco_model_nhdz.geojson`, koppeltabel CSV. Immutable, written once. | +| `pwn_layer_model` | module (test_pwnlayers3_layers) | The one expensive result: `get_pwn_layer_model(...)` on a 4×4 synthetic vertex ds with `nlmod.read.regis.get_layer_names` monkeypatched. Read-only; the plot test takes a `.copy()`. | +| `chloride_nc` | session | Tiny `chloride_p50.nc`: 2 source layers (1-D `top/bottom` coords, power-of-two thicknesses), 3×3 grid at 250 m, **descending y** (matches the real file — ascending-y-only fixtures would pass while production fails). | + +### 3.3 pyproject changes + +- Add `pyarrow` to `[project] dependencies` (production requirement of `well.py:46`). +- `[tool.pytest.ini_options]`: `addopts = "--strict-markers"`, + `markers = ["network: live web services, excluded in CI", "mf6: needs MODFLOW binaries"]`, + `env = ["MPLBACKEND=Agg"]` (finally using the declared pytest-env), + a named `filterwarnings` entry for the pyarrow feather FutureWarning. +- Fix the stale `tests/test_pwnlayers.py` reference in the lint script (pyproject.toml:85). +- Delete the empty shells `test_nhflo_utils.py` / `test_hhnk.py` (out of the 09pwnmodel2 + closure; they misrepresent coverage). Local housekeeping (not CI): `hatch env remove test` + once; `tests/model_ws/` (51 MB untracked gridgen leftovers) can be deleted. + +### 3.4 Mocking policy (named seams only) + +Mock only at named network seams, with the reason stated in the test: +`nlmod.read.waterboard.download_data` (live HHNK ArcGIS REST, uncached at polder.py:40); +`nlmod.read.regis.get_layer_names` (OPeNDAP when called without ds); +`nlmod.plot.add_background_map` (contextily tile fetch); +`nlmod.read.rws.get_gdf_surface_water` (synthetic gdf substitution — one patch controls +`discretize_northsea`/`discretize_surface_water` too, verified major_surface_waters.py:42); +`nlmod.gwf.output.get_heads_da` / `nlmod.gwt.output.*` in the transport-derivation tests +(real loader path covered once by the MF6 smoke test). Everything else — `gdf_to_grid`, +`aggregate_vector_per_cell`, `build_spd`, `split_layers_ds`, `discretize_*`, +`calculate_thickness`, `get_isosurface` — runs the **real nlmod code**, so the suite doubles +as an nlmod-@dev compatibility canary (that canary already caught the +`get_isosurface(left=...)` staleness failure). + +### 3.5 Assertion policy + +Exact (`assert_array_equal` / `==`) by default; geometry axis-aligned on grid-multiple +coordinates so shapely areas are exact floats; expected values derived in-test from first +principles, never pasted from implementation output. The only tolerances, each named: +IMS solver `dvclose` (MF6 smoke, atol=1e-8), `rtol=1e-12` on area-weighted means (pandas +groupby accumulation order vs the test's reconstruction), `rtol=1e-12` on linear-griddata +plane reproduction (Qhull barycentric arithmetic). + +**Baseline-validation protocol** (global CLAUDE.md doctrine): every test tagged +`[regression: ]` must be run against the pre-fix baseline +(`git checkout ^ -- src/nhflotools/` in a scratch worktree) and confirmed to +FAIL there, then confirmed to pass on HEAD, before it is trusted. + +--- + +## 4. Test files and per-test specs + +### 4.1 `tests/test_major_surface_waters.py` (new) + +1. **test_get_chd_ghb_data_known_answer_cond_area** — 2×2 vertex ds (100 m cells); + synthetic gdf: sea polygon (`OWMNAAM='Hollandse kust (kustwater)'`) over cell 0, + `'IJsselmeer'` (peil=-0.2, bweerstand=2.0) exactly covering cell 3. Assert + `rws_oppwater_cond[3] == 10000/2.0`, `area[3] == 10000.0`, `stage[3] == -0.2`, dry + cells `== 0.0` (pins the zeros-not-NaN contract downstream code relies on). Catches the + conductance-formula bug class (the panden 83365c7 1e4× error, identical code shape). +2. **test_stage_override_only_where_ijsselmeer_intersects** — Noordzeekanaal wins cond in + cell 2; IJsselmeer (peil=p) also intersects cell 2 but not cell 3. Assert + `stage[2] == p` and `stage[3]` bit-identical to winner-take-all (the + `xr.where(da_peil.isnull())` identity). Catches override-mask inversion. +3. **test_input_ds_mutation_northsea_and_extrapolation** — plant an all-NaN `botm/kh` + column under the sea cell. Afterwards: input ds has `northsea == 1` exactly there, and + the column equals its Euclidean-nearest valid neighbor exactly. Catches the + discarded-`extrapolate_ds`-return hazard (in-place mutation contract, nlmod + dims/base.py:286-301): a copy-semantics change in nlmod leaves NaNs under the sea and + this fails loudly. +4. **test_chd_ghb_float_branch_masks_disjoint_layers** — `chd_ghb_from_major_surface_waters` + on `vertex_ds` + `gwf_disv`; cells: sea / wet / two dry; one wet cell `idomain[0]=0`; + **plus one sea cell `idomain[0]=0` among active sea cells** (must be absent from CHD — + the partial-absence case). `sea_stage=0.04`. Assert: `sfw_stage` NaN exactly where + `northsea==1`, `== rws_oppwater_stage` elsewhere; `sfw_cond == 0.0` on sea; **GHB ∩ CHD + cell sets = ∅** (boundary-type disjointness — a cell that is both GHB and fixed-head is + non-physical); every CHD rec head `== 0.04`, aux `== 18000.0` (`SEA_CHLORIDE_MG_L`); + `ts_sea is None`; the idomain-blocked GHB row sits in layer 1 (first active), the + blocked CHD cell absent (layer-0-only contract). One setup, six bug classes. +5. **test_chd_ts_branch_wiring** — `sea_stage=[(0.0, 0.04), (365.0, 0.04)]`. Assert every + CHD rec head is the **literal string** `'sea_stage'`; `chd.ts` has + `time_series_namerecord == 'sea_stage'`, linear interpolation, timeseries equal to the + input list (a flat series is semantically the constant branch). + `[regression: b42fec1]` — a list-valued sea_stage never worked (NameError + wrong + namerecord). +6. **test_chd_none_guard_returns_three_tuple** — all sea cells `idomain[0]=0`, list-branch + call. Assert `(ghb, None, None)`, no raise. `[regression: 9d6c40c]` + (AttributeError on `chd.ts.initialize`). +7. **test_gdf_permutation_invariance** (lowest priority in file) — permute gdf rows; all + three output vars bit-identical (winner-take-all must be order-free; no area ties by + construction). Catches a last-wins loop rewrite. Cheap. + +*Flagged, deliberately not tested (would entrench defects):* the +`isinstance(sea_stage, float)` int-rejection quirk; GHB silently zeroed where sea and +canal overlap. → §8 issues. + +### 4.2 `tests/test_nhi_chloride.py` (new) + +Uses `chloride_nc`; called without `cachedir` (decorator short-circuits) except test 7. + +1. **test_constant_field_identity** — source ≡ 1000 covering the full vertical extent, + no sea → output == 1000.0 everywhere, layers preserved, `units == 'mg/l'`. Any + non-convex aggregation fails. Exact. +2. **test_weighted_mean_known_answer** (parametrized ×2) — (a) model layer `[0,-8]`, voxels + c=100/300 → exactly `(4*100+4*300)/8 = 200.0`; (b) layer `[-2,-6]` clipping 2 m of each + voxel → `200.0` with different weights; a `[0,-4]`-only layer → `100.0`. Catches + clip/weighting off-by-ones. Exact (integer-valued arithmetic). +3. **test_nan_voxel_excluded_from_denominator** — one voxel NaN → column output `== c1` + exactly, not diluted toward 0. Catches the skipna-numerator/keep-denominator bug (the + exact bug that *does* exist at postprocessing.py:182 — see §8). +4. **test_sea_override_layer0_only** — sea cell: layer 0 `== 18000.0` exactly; deeper + layer unchanged. Catches positional `da.values[0]` misalignment. +5. **test_deep_layer_bfill_ffill** — model layers wholly below and wholly above the source + stack → output NaN-free; deep layer equals last covered value, top layer equals first + covered value (asymmetric column pins bfill-before-ffill order). + `[regression: 83365c7]` — all-NaN deep layers poisoned the transport IC. + **Contingency:** the recon flags a possible empty-valid-points crash in the fillnan + step *before* the bfill/ffill rescue (nhi_chloride.py:44-52 ordering). If this test + crashes rather than fails, that is a live bug find — report it, don't reshape the + fixture to avoid it. +6. **test_vertex_grid_path** — 2×2 pseudo-vertex ds (production grid type): output dims + `(layer, icell2d)`; one cell outside the source's horizontal extent → NaN → vertex + `fillnan` path fills it with its nearest neighbor's value exactly. Covers the + production-only dispatch that structured tests never touch (a constant field would + leave nothing to fill and make this a near-no-op). +7. **test_cache_roundtrip_identical** — attrs `gridtype/extent`, + `cachedir=tmp_path, cachename='chloride'`: `chloride.nc` + `.pklz` written; second call + bit-identical values. Pins the production cachedir path (call site always passes it). + +### 4.3 `tests/test_panden.py` (fill the 0-byte file) + +Real tiny shapefile via `gpd.to_file` (exercises the real read path); structured ds must +include `kh` (build_spd reads it unconditionally). + +1. **test_get_oppervlakte_pwn_shapes_name_filter** — Naam + `['ICAS-noord','IKIEF-3','VIJVER']`: stages exactly `{2.8, 5.8}`, `c == 1.0`, + `rbot == stage - 2.0` elementwise, VIJVER row absent, caplog warning. + `[regression: e16b418 (#41)]` — VIJVER used to abort the model run. +2. **test_riv_known_answer_single_cell** — 50×50 ICAS square inside one cell → exactly one + spd record `[(0,iy,ix), 2.8, 2500.0, 0.8, 0.0, 'ICAS…']` (aux before boundname). + `[regression: 83365c7]` — baseline aggregated resistance `c` instead of `area/c`; cond + would be `1.0`, a silent ~1e4× physics error. +3. **test_riv_cond_conservation_across_cells** — one rectangle straddling two cells (areas + 2500/7500, axis-aligned → exact). Per-cell cond == intersected area and + `sum(cond) == 10000.0` exactly (partition conservation). Same baseline, different + failure mode (per-piece vs per-cell aggregation). +4. **test_riv_mixed_names_cell** — ICAS + IKIEF sharing one cell. Stage + `== (2.8*a1+5.8*a2)/(a1+a2)` (rtol=1e-12, named: groupby accumulation order), + `rbot == 0.8` (min), `cond == a1+a2`; assert build_spd's stage-clip warning did not + fire. Characterizes the mixed-cell aggregation mapping. +5. **test_riv_lay_of_rbot_placement** — 2-layer ds `botm=[1.0, -10]`, rbot=0.8 → record in + layer **1**, not 0. Catches the first-crossing off-by-one. +6. **test_riv_empty_intersection_returns_none** — polygons outside extent → `None`, + `ds.attrs['ssm_sources']` untouched. `[regression: 232de8d]` +7. **test_riv_ssm_registration_idempotent** — `transport=1`, call twice → `'riv'` appears + exactly once; `transport=0` → attrs untouched. `[regression: c613863 (#39)]` +8. **test_riv_pregridded_frame_raises** — feed a frame with duplicate index / already + gridded (the #114 generalization): `gdf_to_grid` raises ValueError + (nlmod grid.py:1996-1997). Negative test documenting the single-intersection contract + so the models-repo crash class is pinned on the tools side. + +### 4.4 `tests/test_polder.py` (recreate deleted file) + +All tests monkeypatch `nlmod.read.waterboard.download_data` (named: live HHNK ArcGIS REST, +uncached at polder.py:40). Grid = `vertex_ds` (+ float `top`, `ahn`, `northsea`) + `gwf_disv`. + +1. **test_drn_known_answer_and_cbot_scaling** (parametrized cbot ∈ {1.0, 2.0}) — polygon + fully covering cell A, `summer==winter==s`: `drn_cond[A] == area/cbot` exactly (1/cbot + proportionality across params), `drn_elev[A] == s`. +2. **test_drn_elev_nan_skipping_mean** — stages (1.0, 3.0) → 2.0; (1.0, NaN) → 1.0, + exactly (one-sided nanmean, not griddata). +3. **test_drn_fallback_partition_and_layer** — 4 cells: covered / uncovered land + (`ahn=4.25`) / `northsea==1` / uncovered land with `idomain[0]=0`. Uncovered land gets + `elev == ahn` and `cond == area/cbot` (**nonzero** — the historical silent failure); + sea cell `drn_cond` NaN and **absent from the DRN reclist**; blocked cell's record in + layer 1. `[regression: 83365c7]` — zero-init gave uncovered land zero-conductance + drains and drained the sea. +4. **test_drn_empty_download_returns_none_ds_untouched** — empty GeoDataFrame → `None`; + `drn_elev`/`drn_cond` **not** in ds. `[regression: a048d98]` +5. **test_drn_duplicate_index_dedup** — download index `['A','A','B']` → no ValueError + from `gdf_to_grid`; assert dedup via the area-weighted **stage** of the covered cell + plus record presence. (Do **not** assert cond from intersected areas: verified + polder.py:59 computes `cond = full-cell area / cbot`, so an intersected-area cond + assertion would fail against correct code.) +6. **test_drn_both_stages_nan_nearest_fill** — polygon in cell A with both stages NaN, + valid stage s in cell B → `drn_elev[A] == s` exactly (nearest donor, member of the + valid set). Documents the all-NaN crash boundary without pinning scipy's error. + +*Flagged, not tested:* full-cell (not intersected-area) conductance for sliver overlaps — +design decision to raise upstream, not entrench. → §8. + +### 4.5 `tests/test_well.py` (extend; keep the existing 9 tata tests, add none for tata) + +Real tiny files: geojson via `gpd.to_file` (deliberately exercises the GeoJSON round-trip +that stringifies `sec_nput` — the reason for the coercion at well.py:59), feather via +`to_feather` (needs pyarrow, §1 blocker). + +1. **test_pwn_q_known_answer_conservation** — 3 wells tag T1, `sec_nput=3` (one given as + string `"3"`), feather `T1 = [-30,-30,-30]` m³/h; second tag with odd asymmetric series + `[-10,-20,-60]` on one well (`sec_nput=1`); an `'ophaal tijdstip'` datetime column. + Assert each T1 well `Q == -240.0` exactly, `sum == -720.0 == 24*median` (the WEL + mass-balance contract the MAW branch scales back up); odd-series well + `Q == 24*(-20.0)`; the datetime column never contaminates the median; `rw == 0.25`, + `CONCENTRATION == 0.0`. `[regression: 83365c7]` (`numeric_only` + factor-24/nput). +2. **test_pwn_drops_bad_rows_warns_keeps_infiltration** — rows: unmapped tag, `sec_nput=0` + (inf Q), zero-median tag, and one **positive**-median infiltration well. Assert exactly + the three bad rows dropped, warning count == dropped count, survivors finite and + nonzero, the positive-Q well survives with positive sign (drop mask is sign-symmetric). + `[regression: 83365c7 + c613863]` +3. **test_pwn_flow_product_error_contract** (parametrized) — `'timeseries'` → + `NotImplementedError`, `'bogus'` → `ValueError` (documents that the feather read + precedes dispatch). + +*(Design A's tata deepest-layer proposal was cut: `well.py:143-146` shows the `lay>0` path +is identical for mid and deepest layers — no new branch.)* + +### 4.6 `tests/test_postprocessing.py` (extend; keep the existing 9) + +1. **test_check_budget_real_listing** (parametrized ×3) — listing text from + `write_mf6_listing` (provenance §3.1), parsed by **real** `flopy.utils.Mf6ListBudget`: + `|disc| = 0.99` passes; `= 1.0` raises (inclusive `>=` boundary); garbage file → + `RuntimeError` "Could not parse". Catches budgetkey-string and flopy API drift the + existing fully-mocked test cannot. +2. **test_add_output_freshwater_head_identity** — monkeypatch the three nlmod output + loaders to synthetic `(time=2, layer=3, icell2d=3)` arrays; conc ≡ 0, `drhodc` attr + set. Assert `freshwater_head == head_filled` **exactly** (ρ = ρ_ref identity). Catches + density-correction wiring (z-term sign, drhodc pickup). +3. **test_add_output_constant_conc_identities** — conc ≡ 3000 (between thresholds): + `concentration_mean == 3000.0` exactly regardless of unequal layer thicknesses; + `dconcentration_mean.isel(time=0) == 0` identically; `grensvlak_zoet == ds['top']` and + `grensvlak_brak == botm[-1]` per the pinning rules; one planted head NaN → + `head_filled == head` wherever head is finite (bfill conservation). +4. **test_add_output_zoet_brak_ordering** (graft from Design A) — monotone-with-depth + synthetic conc: `grensvlak_zoet >= grensvlak_brak` everywhere (the fresh interface is + never deeper than the brackish one). +5. **test_interface_bounds_random_profiles** (graft from Design A) — property test, + `np.random.default_rng(0)`, ~50 random profiles: + `botm[-1] <= interface_elevation(...) <= top` per column — bounds currently unasserted + for interpolated columns. Exact bound check. +6. **test_plot_result_maps_filenames_and_early_return** (parametrized) — monkeypatch + `nlmod.plot.add_background_map` (named: contextily tile fetch, the module's only + network call); tiny vertex ds with `freshwater_head/grensvlak_*` (with `threshold` + attrs) + ctop, nper=2, `iper=-1`. Assert figdir contains **exactly** + `{doorsnedelijnen.png, map_head_L0_t1.png, map_conc_L0_t1.png, grensvlak_zoet_t1.png, + grensvlak_brak_t1.png}` (t1 = the iper-normalization known answer). Variant (b): ds + without `freshwater_head`, `ctop=None` → only `doorsnedelijnen.png`, no raise. Variant + (c): ds with `drn_elev` (the call site always has it — polder.py:81-82 sets it; the + branch is live, contrary to the recon's first read) → `oppervlaktewater.png` in the set. + +### 4.7 `tests/test_pwnlayers_get_top.py` (new — `get_top_from_ahn`) + +1. **test_pure_griddata_nearest_known_answer** — anisotropic 5-cell layout where the + x-nearest and y-nearest donors differ; **query strictly off the q1==q2 diagonal and + donors off p1==p2** (a query at the origin is provably invariant under a one-sided + (y,x)↔(x,y) swap — the geometry must make the swap detectable); flags off (touches + zero nlmod code). Assert the NaN cell receives its true Euclidean-nearest value and all + valid cells are bit-identical. Catches a silent axis-order swap during refactoring. +2. **test_no_nan_identity** — NaN-free ahn → output identical, no raise (pins the + empty-qpoints scipy edge, scipy-version-sensitive). +3. **test_fill_priority_peil_constant_partial** — monkeypatch the `nlmod.read.rws` trio to + hand-built Datasets. Three NaN cells: fully-covered water that is also in the northsea + mask → gets **peil** (priority order); sea-only cell → gets `0.0` with + `replace_northsea_with_constant=0.0` (**the falsy-constant guard** — production passes + 0.0; a refactor of the `is not None` test to truthiness silently disables the sea fill + and this fails); partially-covered cell (area < cell area) → falls through to nearest + interpolation (the isclose full-coverage gate). Highest-value test in the file. +4. **test_missing_ahn_raises_valueerror** — one-line error-surface pin (guards silent key + rename). + +### 4.8 `tests/test_pwnlayers3_layers.py` (new — highest-risk module) + +Pure seams first, then the merge, then one offline integration. + +1. **test_fix_missings_botms_both_copies** (parametrized over the two divergent duplicates + `pwnlayers3.layers` and `pwnlayers.utils` — locks them together; named maintenance + trap) — 3-layer × 4-cell botm, mid-column NaN, one botm crossing above the layer above; + expectation hand-computed from ffill + `minimum.accumulate`. Assert exact equality, + NaN-free, monotone non-increasing, `<= top`, **idempotent** (`f(f(x)) == f(x)`), + **input object unmodified** (purity — `[regression: df20a42]`, the caller assumed + in-place mutation and discarded the fix), ValueError on NaN top. +2. **test_get_thickness_telescoping** — `thickness[k] == botm[k-1] - botm[k]` exact, + `sum('layer') == botm[0] - botm[-1]`, labels are the lower layers, W11 absent. Catches + the sign/labeling class `[regression: 479d673]`. +3. **test_guard_zero_thickness** — thickness `[0.0, 1e-12, 0.5]` → fill_value exactly at + the isclose-zero entries, others untouched; identity when no zeros. +4. **test_kh_kv_harmonic_and_anisotropy** (graft from Design A — the physics core: + resistances → conductivities) — 2×2 vertex ds + `pwn_data_tree` variant with constant-c + conductance polygons. Assert: single constant-c polygon fully covering the grid → + `kv == d/c` and `kh == d·anisotropy/c` exactly; two half-covering polygons c1=2, c2=4 → + `inv_c == 0.5·(1/2 + 1/4) == 0.375` exactly; W layers `kv == kh/anisotropy` exactly on + the mask; S layers on the c-path `kh/kv == anisotropy` exactly where thickness > 0. +5. **test_kh_nhdz_branch_known_answer** (critique gap) — a cell inside the + `triwaco_model_nhdz` region: `kh == KD/d` from the point-griddata path + (layers.py:857-875); also covers the empty-Bergen-region 0-cell GridIntersect hazard + (layers.py:887-889) by keeping one region empty and asserting no crash and no + contamination. +6. **test_interpolate_da_nearest_linear_and_view_mutation** — 4-cell line; + `_interpolate_da` middle cell missing → nearest donor value; linear midpoint == mean of + the two valid values (representable); `ismissing` empty → strict no-op; **assert the + parent Dataset's array actually changed** (the `.loc`-on-`sel`-view contract — an + xarray copy-semantics change turns transition interpolation into a silent no-op). +7. **test_combine_identity_reduction** — synthetic REGIS (NaN-free) + OTHER with all-False + mask/transition, pure-1:1 koppeltabel → output `kh/kv/botm` bit-identical to REGIS, + categories all 1. Baseline sanity for the whole merge. +8. **test_combine_routing_and_split_conservation** — koppeltabel with 1:1, 1:2 (REGIS + split), 2:1 (OTHER split), one NaN-uncoupled deep row; mask True on half the cells. + Assert: category-2 cells carry OTHER values exactly, category-1 REGIS exactly; per + cell, sum of split-sublayer thicknesses == original layer thickness (exact — same-float + subtraction chains); group bottoms preserved; uncoupled layers category 1. Catches the + positional koppeltabel/split-alignment hazard — the scariest silent-corruption vector + in the codebase. +9. **test_combine_transition_convexity** — transition band between mask and REGIS regions: + interpolated transition kh/botm lie within `[min, max]` of surrounding category≠3 + values **and differ from their pre-merge REGIS values** (proves interpolation actually + ran; complements test 6). +10. **test_ratios_forward_inverse** — `_compute_thickness_ratios`: ratios sum to 1 per + cell exactly, reproduce actual fractions, `1/N` at zero-thickness; + `_apply_ratios_to_botm` with equal ratios reproduces `split_layers_ds`'s equal-split + botms exactly (forward+inverse identity). Catches the group-top `first_idx-1` + off-by-one. +11. **test_get_pwn_layer_model_offline_integration** (module fixture `pwn_layer_model`) — + 4×4 synthetic vertex ds_regis (3 fake REGIS layers + `'mv'`), + `nlmod.read.regis.get_layer_names` monkeypatched (named: OPeNDAP in the first line), + `pwn_data_tree` with one conductance polygon deliberately `VALUE=0`. Assert + postconditions: `kh/kv/botm` NaN-free; **`kh, kv > 0` everywhere** (the zero/inf guard + → NaN → fill path; `[regression: df20a42/5daf902]` zero-kh NPF crash); botm monotone + non-increasing; output `top` identical to input; W-layer `kv == kh/10` exactly on + masked cells; diagnostics `cat_botm/botm_method/kh_method/kv_method` present **with + `flag_values`/`flag_meanings` attrs** (the plot module's hard requirement); botm + source points on an affine plane z=ax+by+c reproduced at interior cell centres + (rtol=1e-12, named: Qhull barycentric arithmetic). Parametrized negatives: NaN in top → + ValueError; ds_regis missing a layer → ValueError. Budget ≈ 1–3 s. +12. **test_area_passthrough_values** — ds_regis with an `'area'` var → output area values + identical; without → computed. `[regression: d62ebbf + 83365c7]` (eager-default + KeyError). **Values-passthrough only** — do NOT assert `get_area` is not called: + verified layers.py:232 `ds_regis.get("area", nlmod.dims.get_area(ds_regis))` evaluates + the default eagerly even when `area` exists, so a spy assertion fails on HEAD. The + eager evaluation is filed as a fix candidate instead (§8). + +### 4.9 `tests/test_pwnlayers3_plot.py` (new) + +1. **test_parse_flag_labels_known_answers** (parametrized) — the three **literal** + `flag_meanings` strings copied from layers.py into the test as the independent source → + exact label lists of lengths 5/7/4, each matching its `flag_values` length (the gate at + plot.py:243 silently drops all tick labels on mismatch); adversarial cases: trailing + `;`, entry without `N:` prefix, text after `(` truncated, underscores→spaces. Pure, + instant. +2. **test_load_project_and_overlay** — tmp `botm/botm.geojson` (EPSG:28992), points at + perpendicular distances 30 m and 80 m from line `[(0,0),(100,0)]`, + `buffer_distance=50`, one point at exactly 50 (boundary `<=`); shuffled non-RangeIndex + and non-layer column order. Assert exact `d_along` projections, far-point exclusion, + boundary inclusion, returned keys == layer_names ∩ columns **in layer_names order** + (positional-alignment hazard). Then feed the dict (plus one all-NaN layer and one z + outside `[zmin, zmax]`) to `_overlay_source_botm` on an Agg Axes: number of + PathCollections == layers with ≥1 valid point; no NaN/out-of-range offsets. +3. **test_plot_diagnostic_cross_sections_end_to_end** — reuse `pwn_layer_model` (`.copy()`) + + ds_regis with `xv/yv/icvert` injected (replicating the call-site mutation), midline + through the synthetic extent, `data_path_2024=pwn_data_tree` (the call site always + passes it — 01_pwnmodel2.py:604). Assert `(fig, axes)` with `axes.size == 6` and cat + colorbar ticklabels exactly `['REGIS', 'PWN', 'Transition']`; close fig. Justified + against the no-smoke rule: the named KeyError regressions (missing flag attrs, + `layer_pwn` rename, `return_diagnostics=False` ds) all make precisely this call raise. + +### 4.10 `tests/test_nhflodata_contract.py` (new) + +1. **test_mockup_paths_and_koppeltabel_columns** (parametrized ~10 cases) — with + `NHFLODATA_LOCATION` deleted, `get_abs_data_path(name, 'latest')` + the hardcoded + relative file must exist for every (dataset, file) pair nhflotools reads: + `Panden_ICAS_IKIEF.shp`, `pumping_infiltration_wells.geojson`, `sec_flows.feather`, + both tata geojsons, `chloride_p50.nc`, `bodemlagenvertaaltabelv2.csv` (also read the + 20 KB and assert columns `'Regis II v2.2'` / `'ASSUMPTION1'`), `botm/botm.geojson`, + `boundaries/S11/S11.geojson`, `boundaries/triwaco_model_nhdz.geojson`. Since + `get_abs_data_path` only **warns** on missing paths (get_paths.py:109-110), this is the + sole automated guard against a data-repo restructure silently breaking 09pwnmodel2. + Instant (stat calls; the 31 MB nc is never opened). + +### 4.11 `tests/test_mf6_smoke.py` (new, marker `mf6`) + +1. **test_tiny_run_budget_and_output_loading** — nlmod test_015 recipe: 3×2 cells × 5 + layers structured ds, CHD=1.0 on the edge mask, `write_and_run`. Then, on the **real** + files: `check_budget_discrepancy(ws, name, transport=False)` passes (an all-CHD steady + model closes its budget); `add_output_to_ds(ds, ws, name, transport=False)` returns + `(ds, None)` with `ds['head'] == 1.0` everywhere (atol=1e-8 — named: IMS iterative + solver converged to dvclose, not algebraic exactness). The one end-to-end pipeline test + with a physical invariant, and the only coverage of the real `.hds`/`.grb`/listing + loader paths the monkeypatched tests bypass. Its committed `.lst` (~2 KB) is the + provenance source for `write_mf6_listing` (§3.1). ~3–8 s. + +--- + +## 5. Coverage map (every 09pwnmodel2-used function) + +| Function | Disposition | +|---|---| +| `pwnlayers.layers.get_top_from_ahn` | §4.7 (4 tests) | +| `pwnlayers3.layers.get_pwn_layer_model` (+ get_ds/get_botm/get_kh/get_kv, combine, fix_missings ×2) | §4.8 (12 tests) | +| `major_surface_waters.get_chd_ghb_data_from_major_surface_waters` | §4.1.1–3, 7 | +| `major_surface_waters.chd_ghb_from_major_surface_waters` | §4.1.4–6 | +| `nhi_chloride.get_nhi_chloride_concentration` | §4.2 (7 tests) | +| `well.get_wells_pwn_dataframe` | §4.5 (3 tests) | +| `well.get_wells_tata_dataframes` | existing 9 tests; **no additions** (nearest-cell, kd threshold strictness, chloride-warning boundary, screen offsets, IndexError contract already covered; more violates leanness) | +| `polder.drn_from_waterboard_data` | §4.4 (6 tests) | +| `panden.riv_from_oppervlakte_pwn` + `get_oppervlakte_pwn_shapes` | §4.3 (8 tests) | +| `pwnlayers3.plot.plot_diagnostic_cross_sections` (+ helpers) | §4.9 (3 tests) | +| `postprocessing.check_budget_discrepancy` | existing + §4.6.1 + §4.11 real run | +| `postprocessing.add_output_to_ds` | §4.6.2–5 + §4.11 real run | +| `postprocessing.plot_result_maps` | §4.6.6 | +| `postprocessing.interface_elevation` | existing 8 tests + §4.6.5 bounds property | + +**Named exclusions:** + +- **Full 09pwnmodel2 pipeline / transport run** — needs live REGIS, AHN5, HHNK services + plus a minutes-long MF6 flow+transport run; not CI-realistic. The pipeline obligation is + discharged by §4.11 (real run + budget closure) and §4.8.11 (offline layer-model + integration). A future opt-in `network`-marked nightly is noted, not planned. +- **`pwnlayers3.layers.get_top`** — dead at the call site (script imports + `get_top_from_ahn`); default path hits `download_bathymetry` (network). +- **Legacy `pwnlayers.get_bergen_botm` / `get_mensink_botm`** history candidates (883e517, + 479d673) — outside the 09pwnmodel2 closure; the sign/labeling bug class is covered by + §4.8.2's analog in pwnlayers3. +- **numba-accelerated `get_isosurface` path** — CI env has no numba, so only the numpy + fallback is exercised; numba/numpy divergence is untestable there. Stated, not hidden. +- **`get_transition` monotonicity** — acceptable drop; get_ds's internal validation + (layers.py:400-417) executes inside §4.8.11. +- **`recharge_utils`** — lives in the models repo, zero nhflotools calls, has its own test + file there; port only if it migrates to tools. +- **hhnk.py, nhflo_utils.py, bofek.py, geoconverter, bergen/berging utils** — not in the + 09pwnmodel2 closure; out of scope by the task's restriction. + +--- + +## 6. CI workflow + +New `.github/workflows/test.yml` beside lint.yml: + +```yaml +name: test +on: + push: {branches: [main]} + pull_request: +jobs: + test: + runs-on: ubuntu-latest + env: + MPLBACKEND: Agg + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: {python-version-file: pyproject.toml} # 3.14, same as lint.yml + - name: Resolve moving deps for cache key + id: deps + run: | + echo "data_sha=$(git ls-remote https://github.com/NHFLO/data.git main | cut -f1)" >> "$GITHUB_OUTPUT" + echo "nlmod_sha=$(git ls-remote https://github.com/gwmod/nlmod.git dev | cut -f1)" >> "$GITHUB_OUTPUT" + - uses: actions/cache@v4 # uv download cache incl. the ~235 MB nhflodata build + with: + path: ~/.cache/uv + key: uv-${{ hashFiles('pyproject.toml') }}-${{ steps.deps.outputs.data_sha }}-${{ steps.deps.outputs.nlmod_sha }} + - run: pipx install hatch + - name: Sync test env + run: hatch env create test # uv installer per [tool.hatch.envs.default] + - uses: actions/cache@v4 # MODFLOW binaries + flopy exe metadata + id: mfcache + with: + path: | + ~/mfbin + ~/.local/share/flopy + key: mfbin-${{ steps.deps.outputs.nlmod_sha }} + - name: Download MODFLOW binaries + if: steps.mfcache.outputs.cache-hit != 'true' + run: hatch run test:python -c "import os, nlmod; nlmod.util.download_mfbinaries(bindir=os.path.expanduser('~/mfbin'))" + - name: Run tests + run: hatch run test:test -- -m "not network" --durations=15 +``` + +Notes (each verified during review): + +- `~/.local/share/flopy` (the `get_modflow.json` exe metadata flopy writes outside the + bindir) must be cached alongside `~/mfbin`, otherwise a cache hit skips the download and + exe resolution fails. `bindir` needs `os.path.expanduser` — a literal `'~/mfbin'` is not + tilde-expanded by Python. (Confirm the metadata path on the runner at implementation + time; it is platform-dependent.) +- `NHFLODATA_LOCATION` deliberately unset (and force-deleted in conftest) → mockup + resolution inside the installed wheel is what is exercised. +- The `git ls-remote` SHAs make the unpinned `nlmod@dev` / `nhflodata@main` deps + cache-correct (new upstream commit → fresh env) while caching between pushes. The + observed stale-env failure class (nlmod 0.9.1dev lacking `get_isosurface(left=...)`) is + exactly what this prevents. Longer-term: pin both to SHAs in pyproject. +- Both repos are public — no tokens needed. +- Budget: pytest total ≤ 60 s (estimate ~35 s); job total ≤ 6 min cold, ≤ 2.5 min warm. + +--- + +## 7. Implementation order (value-ranked from the regression analysis) + +1. `test_panden.py` — silent 1e4× physics error class, empty file today + (baselines 83365c7/e16b418/232de8d/c613863) +2. `test_polder.py` — silent zero-cond drains + sea drains (83365c7/a048d98) +3. `test_well.py` PWN additions — silent mass-balance factor (83365c7/c613863); + **prerequisite: add pyarrow to runtime deps** +4. `test_major_surface_waters.py` — crash class + boundary-condition correctness + (9d6c40c/b42fec1) +5. `test_nhi_chloride.py` — transport-IC poisoning (83365c7) +6. `test_pwnlayers3_layers.py` — deepest silent-corruption surface (df20a42/5daf902/479d673) +7. `tests/util.py` + `conftest.py` + CI workflow — prerequisite for 1–6, build alongside 1 +8. `test_pwnlayers_get_top.py`, `test_postprocessing.py` additions, `test_pwnlayers3_plot.py` +9. `test_nhflodata_contract.py`, `test_mf6_smoke.py` (smoke test also generates the + committed `.lst` fixture used by §4.6.1) + +For every `[regression: ]` test: before merging, restore the pre-fix file +(`git checkout ^ -- src/nhflotools/` in a scratch worktree), confirm the +test FAILS, revert, confirm it passes on HEAD. + +## 8. Issues found during review (NOT tests — tests must not entrench them) + +All verified against HEAD and filed on GitHub. + +| # | Finding | Status | +|---|---|---| +| 1 | `pd.read_feather` needs pyarrow, not a declared dependency (well.py:46) | **Fixed in this PR** — added to runtime deps | +| 2 | `concentration_mean` NaN/thickness-denominator bias (postprocessing.py:182) | [#60](https://github.com/NHFLO/tools/issues/60) | +| 3 | Eager `get_area` default at pwnlayers3/layers.py:232 | [#61](https://github.com/NHFLO/tools/issues/61) | +| 4 | `isinstance(sea_stage, float)` mishandles an int stage (major_surface_waters.py:102) | [#62](https://github.com/NHFLO/tools/issues/62) | +| 5 | Polder conductance uses full-cell, not intersected, area (polder.py:59) | Already filed as [#51](https://github.com/NHFLO/tools/issues/51) | +| 6 | Sea override precedes horizontal fill, seeding coastal land with 18000 mg/l | [#63](https://github.com/NHFLO/tools/issues/63) | +| 7 | Possible crash on an all-NaN layer before the bfill/ffill rescue (nhi_chloride.py:44-52) | Contingency — file only if §4.2.5 surfaces it | +| 8 | Housekeeping: stale lint ref, orphaned fixtures, hardcoded local path | [#64](https://github.com/NHFLO/tools/issues/64) | + +--- + +## 9. Implementation progress + +Status legend: ☐ not started · ◐ in progress · ☑ done and green. + +### Infrastructure + +| Item | Status | +|---|---| +| `pyarrow` added to runtime dependencies | ☑ | +| `tests/util.py` — vertex-grid builder etc., verified against real nlmod | ☑ | +| `tests/conftest.py` — hygiene, `vertex_ds`, `gwf_disv` | ☑ | +| pyproject: markers, `MPLBACKEND=Agg`, strict markers, lint-target fix | ☑ | +| GitHub issues for §8 findings | ☑ (#60–#64) | +| Removed empty shells `test_hhnk.py` / `test_nhflo_utils.py` | ☑ | +| `.github/workflows/test.yml` | ☑ | +| README module overview (used-by-09pwnmodel2 vs untested) | ☑ | + +### Test files + +All green. Suite total: **125 passed, 1 xfailed in 4.2 s**; slowest single test 0.59 s. + +| File | Plan § | Tests | Status | +|---|---|---|---| +| `test_panden.py` | 4.3 | 7 | ☑ | +| `test_polder.py` | 4.4 | 4 | ☑ | +| `test_well.py` (extended) | 4.5 | 13 | ☑ | +| `test_major_surface_waters.py` | 4.1 | 7 | ☑ | +| `test_nhi_chloride.py` | 4.2 | 8 | ☑ | +| `test_pwnlayers3_layers.py` | 4.8 | 15 | ☑ (1 xfail, #65) | +| `test_pwnlayers_get_top.py` | 4.7 | 3 | ☑ | +| `test_postprocessing.py` (extended) | 4.6 | 17 | ☑ | +| `test_pwnlayers3_plot.py` | 4.9 | 10 | ☑ | +| `test_nhflodata_contract.py` | 4.10 | 41 | ☑ | +| `test_mf6_smoke.py` | 4.11 | 1 | ☑ | + +Every test file was mutation-checked during implementation: a plausible bug was introduced +into the source, the test confirmed to fail, and the source restored. Highlights — the +panden conductance formula reverted to the pre-83365c7 resistance aggregation (caught by +2 tests), the polder fallback conductance zeroed (caught), the northsea fallback guard +removed (caught), `numeric_only` dropped from the well median (caught), the sea-stage +override mask inverted (caught by 3), `extrapolate_ds`'s in-place contract broken (caught), +the layer-split group top misidentified (caught), the `get_top_from_ahn` nearest-donor +metric axis-swapped (caught), the falsy-constant sea guard introduced (caught), the budget +threshold inverted and the budget key renamed (both caught), and the `.hds` head field +mirrored on the x axis (caught by the MF6 run). + +### Deviations from the plan + +- **§4.3 #8** (pre-gridded frame raises) dropped: unreachable through panden's real code + path, since `get_oppervlakte_pwn_shapes` always builds its frame from `gpd.read_file`. + It would have asserted nlmod's own guard, not nhflotools behaviour. +- **§4.4** planned 6 tests, landed 4: the cbot known-answer merged into the partition test + (parametrized), the NaN-mean and nearest-fill cases merged into one, and the all-NaN + crash-boundary case dropped as untestable without pinning scipy internals. +- **§4.7** planned 4, landed 3: the no-NaN identity folded into the nearest-donor test. +- **Unit tests use the vertex (DISV) grid** from `tests/util.py` rather than the structured + grid the plan sketched, so cellids are `(layer, icell2d)` pairs. This matches production + (09pwnmodel2 is a refined vertex model) and needs no gridgen. +- **§4.2 contingency did not materialise**: the suspected all-NaN-layer crash before the + bfill/ffill rescue does not reproduce on this environment, so §8 item 7 was not filed. +- **§4.8 REGIS-guard test** revealed a real defect (issue #65) and is committed as + `xfail(strict=True)`, so it becomes a visible failure the moment the guard is fixed. + +### Further findings filed during implementation + +| # | Finding | +|---|---| +| [#65](https://github.com/NHFLO/tools/issues/65) | REGIS-completeness guard raises a pandas error in the case it exists for | +| [#66](https://github.com/NHFLO/tools/issues/66) | `chd_ghb_from_major_surface_waters` always returns `ts_sea=None` | +| [#67](https://github.com/NHFLO/tools/issues/67) | Budget check inspects only the incremental budget, never the cumulative | +| [#68](https://github.com/NHFLO/tools/issues/68) | `plot_result_maps` leaks all but the grensvlak figures | +| [#69](https://github.com/NHFLO/tools/issues/69) | Panden SSM guard cannot dedupe on a repeated call against one `gwf` | + +Reported but not filed (lower value, recorded here): `polder.py` `si.griddata` raises when +every cell's elevation is NaN; `well.py` silently maps every well to NaN when the flow-tag +namespaces differ, and does not check `locatie` uniqueness; `panden.py` `str.contains` +yields NaN for a NULL `Naam` and raises on the boolean mask; the two copies of +`fix_missings_botms_and_min_layer_thickness` have already diverged in their logging +arithmetic; `_parse_flag_labels` splits on `;` before stripping parentheticals, so a `;` +inside a description silently splits one flag into two labels. diff --git a/pyproject.toml b/pyproject.toml index 75d4731..9220bb5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,6 +41,7 @@ dependencies = [ "numpy", "pandas", "pedon", + "pyarrow", # pandas feather I/O in well.py is gated behind this optional pandas dep "pykrige", "pyshp", # flopy shapefile I/O in pwnlayers is gated behind this optional flopy dep "scipy", @@ -82,10 +83,10 @@ format = [ "ruff format src/nhflotools/*.py tests", ] lint = [ - "ruff check tests/test_pwnlayers.py src/nhflotools/pwnlayers src/nhflotools/geoconverter src/nhflotools/major_surface_waters.py src/nhflotools/nhi_chloride.py src/nhflotools/panden.py src/nhflotools/polder.py src/nhflotools/well.py", + "ruff check tests src/nhflotools/pwnlayers src/nhflotools/geoconverter src/nhflotools/major_surface_waters.py src/nhflotools/nhi_chloride.py src/nhflotools/panden.py src/nhflotools/polder.py src/nhflotools/well.py", ] lintminimal = [ - "ruff check src/nhflotools/pwnlayers src/nhflotools/geoconverter src/nhflotools/major_surface_waters.py src/nhflotools/nhi_chloride.py src/nhflotools/panden.py src/nhflotools/polder.py src/nhflotools/well.py --config \"lint.select=['E4', 'E7', 'E9', 'F']\"", + "ruff check tests src/nhflotools/pwnlayers src/nhflotools/geoconverter src/nhflotools/major_surface_waters.py src/nhflotools/nhi_chloride.py src/nhflotools/panden.py src/nhflotools/polder.py src/nhflotools/well.py --config \"lint.select=['E4', 'E7', 'E9', 'F']\"", ] [tool.hatch.envs.test] @@ -94,4 +95,10 @@ features = ["test"] test = "pytest -v" [tool.pytest.ini_options] +addopts = "--strict-markers" +env = ["MPLBACKEND=Agg"] +markers = [ + "mf6: requires the MODFLOW 6 executable", + "network: hits live web services; excluded in CI", +] testpaths = ["tests"] diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..eae2b35 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,41 @@ +"""Shared fixtures for the nhflotools test suite.""" + +import matplotlib.pyplot as plt +import pytest +import xarray as xr + +from .util import make_gwf_disv, make_rect_vertex_ds + + +@pytest.fixture(autouse=True) +def _hygiene(monkeypatch): + """Keep tests independent: no leaked figures, open files or data-location env var. + + ``NHFLODATA_LOCATION`` is deleted so a developer's local data mount can never + redirect the tests away from the mockup data that CI resolves. + + ``FILE_CACHE.clear()`` already closes the netCDF handles xarray holds open; an + additional ``gc.collect()`` here cost ~50 ms per test (two thirds of the whole + suite's runtime) and closed nothing further, so it is deliberately absent. + """ + monkeypatch.delenv("NHFLODATA_LOCATION", raising=False) + yield + plt.close("all") + xr.backends.file_manager.FILE_CACHE.clear() + + +@pytest.fixture +def vertex_ds(): + """A fresh 2x2 vertex model dataset. + + Function-scoped on purpose: the functions under test mutate the dataset in place + (adding ``northsea``, ``sfw_*``, ``drn_*``, ``thickness``), so sharing one would + couple tests. Building it costs well under a millisecond. + """ + return make_rect_vertex_ds() + + +@pytest.fixture +def gwf_disv(vertex_ds, tmp_path): + """``(ds, gwf)`` with sim/tdis/gwf/disv built in memory from ``vertex_ds``.""" + return make_gwf_disv(vertex_ds, tmp_path) diff --git a/tests/test_hhnk.py b/tests/test_hhnk.py deleted file mode 100644 index d007b83..0000000 --- a/tests/test_hhnk.py +++ /dev/null @@ -1 +0,0 @@ -"""Test HHNK tools.""" diff --git a/tests/test_major_surface_waters.py b/tests/test_major_surface_waters.py new file mode 100644 index 0000000..1ae0c17 --- /dev/null +++ b/tests/test_major_surface_waters.py @@ -0,0 +1,226 @@ +"""Tests for :mod:`nhflotools.major_surface_waters` (RWS sea/lake CHD + GHB).""" + +import nlmod +import numpy as np +import pytest +from shapely.geometry import box + +from nhflotools.major_surface_waters import ( + chd_ghb_from_major_surface_waters, + get_chd_ghb_data_from_major_surface_waters, +) +from tests.util import cell_polygon, make_gdf, make_gwf_disv, make_rect_vertex_ds + +# One of the five OWMNAAM values nlmod.read.rws.discretize_northsea treats as sea +# (nlmod/read/rws.py:214-224). Any other name is a GHB water body, not sea. +SEA_NAME = "Hollandse kust (kustwater)" + + +def _patch_gdf(monkeypatch, gdf): + """Make ``nlmod.read.rws.get_gdf_surface_water`` return ``gdf``. + + The single patch feeds both ``discretize_northsea`` and + ``discretize_surface_water``, because the function under test threads one + GeoDataFrame into both. + """ + + def _fake_get_gdf_surface_water(*_args, **_kwargs): + return gdf + + monkeypatch.setattr(nlmod.read.rws, "get_gdf_surface_water", _fake_get_gdf_surface_water, raising=True) + + +def _inset_cell(ds, icell2d, inset=10.0): + """Return a box strictly inside one cell, shrunk by ``inset`` on every side. + + A polygon snapped to the cell edges also *touches* its neighbours, and + ``discretize_northsea`` rasterises with a plain ``intersects`` predicate, so an + edge-snapped sea polygon would flag the whole grid as sea. Insetting keeps the + sea confined to one cell while the intersected area stays exact. + """ + xmin, ymin, xmax, ymax = cell_polygon(ds, icell2d).bounds + return box(xmin + inset, ymin + inset, xmax - inset, ymax - inset) + + +def test_discretized_cond_is_area_over_bweerstand_and_dry_cells_are_zero(monkeypatch, vertex_ds): + """Conductance equals intersected area / bweerstand; untouched cells are 0.0.""" + ds = vertex_ds + # Sea over cell 0, inset by 10 m -> 80 x 80 m; IJsselmeer covering cell 3 exactly. + gdf = make_gdf( + [_inset_cell(ds, 0), cell_polygon(ds, 3)], + OWMNAAM=[SEA_NAME, "IJsselmeer"], + peil=[0.0, -0.25], + bweerstand=[1.0, 2.0], + ) + _patch_gdf(monkeypatch, gdf) + + rws_ds = get_chd_ghb_data_from_major_surface_waters(ds, cachedir=None) + + area = rws_ds["rws_oppwater_area"].values + cond = rws_ds["rws_oppwater_cond"].values + stage = rws_ds["rws_oppwater_stage"].values + + # 80 m x 80 m inset box in cell 0; the full 100 m x 100 m cell 3. + np.testing.assert_array_equal(area, [6400.0, 0.0, 0.0, 10000.0]) + # cond = intersected area / bweerstand, per water body. + np.testing.assert_array_equal(cond, [6400.0 / 1.0, 0.0, 0.0, 10000.0 / 2.0]) + np.testing.assert_array_equal(stage, [0.0, 0.0, 0.0, -0.25]) + + # Cells without surface water must be exactly zero, never NaN: downstream + # chd_ghb_from_major_surface_waters masks the GHB on ``cond > 0``. + assert not np.isnan(cond).any() + + +@pytest.mark.parametrize("reverse_rows", [False, True]) +def test_ijsselmeer_stage_override_applies_only_where_it_intersects(monkeypatch, vertex_ds, reverse_rows): + """The IJsselmeer peil overwrites the stage only in cells it intersects. + + Row order is permuted as well: the winner-take-all loop in + ``discretize_surface_water`` compares intersected areas strictly, so the result + must not depend on the order the water bodies appear in. + """ + ds = vertex_ds + # Noordzeekanaal covers the whole southern row (cells 2 and 3) -> 10000 m2 each. + # IJsselmeer covers the western half of cell 2 only -> 5000 m2, so it loses the + # area contest there and never reaches cell 3. + rows = [ + (box(0.0, 0.0, 200.0, 100.0), "Noordzeekanaal", -0.5, 1.0), + (box(0.0, 0.0, 50.0, 100.0), "IJsselmeer", -0.25, 4.0), + ] + if reverse_rows: + rows.reverse() + geoms, names, peilen, bweerstanden = zip(*rows, strict=True) + _patch_gdf(monkeypatch, make_gdf(geoms, OWMNAAM=list(names), peil=list(peilen), bweerstand=list(bweerstanden))) + + rws_ds = get_chd_ghb_data_from_major_surface_waters(ds, cachedir=None) + + # Cell 2: overridden to the IJsselmeer peil. Cell 3: keeps the canal peil. + np.testing.assert_array_equal(rws_ds["rws_oppwater_stage"].values, [0.0, 0.0, -0.25, -0.5]) + # The override touches the stage only; the canal keeps winning area and cond. + np.testing.assert_array_equal(rws_ds["rws_oppwater_area"].values, [0.0, 0.0, 10000.0, 10000.0]) + np.testing.assert_array_equal(rws_ds["rws_oppwater_cond"].values, [0.0, 0.0, 10000.0, 10000.0]) + + +def test_input_ds_is_mutated_with_northsea_and_extrapolated_under_sea(monkeypatch): + """``ds`` gains ``northsea`` and its all-NaN sea column is filled in place. + + ``nlmod.dims.extrapolate_ds(ds)`` is called for its side effect only -- the + return value is discarded -- so this pins the in-place mutation contract. + """ + # 3x1 grid so the nearest valid neighbour of the sea cell is unique: + # cell 0 (x=50) -> cell 1 (x=150, 100 m) beats cell 2 (x=250, 200 m). + ds = make_rect_vertex_ds(nx=3, ny=1) + ds["botm"].values[:, 0] = np.nan + ds["kh"].values[:, 0] = np.nan + ds["botm"].values[:, 1] = [-11.0, -21.0] + ds["botm"].values[:, 2] = [-12.0, -22.0] + ds["kh"].values[:, 1] = 8.0 + ds["kh"].values[:, 2] = 16.0 + + _patch_gdf( + monkeypatch, + make_gdf([_inset_cell(ds, 0)], OWMNAAM=[SEA_NAME], peil=[0.0], bweerstand=[1.0]), + ) + + get_chd_ghb_data_from_major_surface_waters(ds, cachedir=None) + + assert "northsea" in ds + np.testing.assert_array_equal(ds["northsea"].values, [True, False, False]) + # The NaN column is replaced by an exact copy of its nearest neighbour, cell 1. + np.testing.assert_array_equal(ds["botm"].values[:, 0], [-11.0, -21.0]) + np.testing.assert_array_equal(ds["kh"].values[:, 0], [8.0, 8.0]) + # Neighbours are untouched. + np.testing.assert_array_equal(ds["botm"].values[:, 2], [-12.0, -22.0]) + + +def _sea_and_lake_ds(tmp_path, pinch_layer0=()): + """Build a 3x2 vertex ds carrying ``northsea`` and the two ``rws_oppwater`` fields. + + Cells 0 and 1 are sea, cells 3 and 4 carry a GHB water body, cells 2 and 5 are + dry. Cells in ``pinch_layer0`` get a zero-thickness top layer, which makes + nlmod derive ``idomain[0] == 0`` there. + """ + ds = make_rect_vertex_ds(nx=3, ny=2) + for icell2d in pinch_layer0: + ds["botm"].values[0, icell2d] = float(ds["top"].values[icell2d]) + ds["northsea"] = ("icell2d", np.array([1, 1, 0, 0, 0, 0])) + ds["rws_oppwater_stage"] = ("icell2d", np.array([0.5, 0.5, 0.0, -1.0, -2.0, 0.0])) + # Sea cells carry a non-zero cond so the sea-zeroing is observable. + ds["rws_oppwater_cond"] = ("icell2d", np.array([100.0, 100.0, 0.0, 250.0, 500.0, 0.0])) + return make_gwf_disv(ds, tmp_path) + + +def test_chd_ghb_float_branch_masks_layers_and_disjointness(tmp_path): + """Float sea_stage: sfw masks, GHB/CHD disjointness, heads, aux and layers.""" + # Cell 1 (sea) and cell 4 (GHB) have a pinched-out layer 0. + ds, gwf = _sea_and_lake_ds(tmp_path, pinch_layer0=(1, 4)) + sea_stage = 0.0625 # exactly representable + + ghb, chd, ts_sea = chd_ghb_from_major_surface_waters(ds, gwf, sea_stage=sea_stage) + + is_sea = ds["northsea"].values.astype(bool) + sfw_stage = ds["sfw_stage"].values + sfw_cond = ds["sfw_cond"].values + # Stage is NaN exactly on the sea, and the raw RWS stage everywhere else. + np.testing.assert_array_equal(np.isnan(sfw_stage), is_sea) + np.testing.assert_array_equal(sfw_stage[~is_sea], ds["rws_oppwater_stage"].values[~is_sea]) + # Conductance is zeroed on the sea, untouched elsewhere. + np.testing.assert_array_equal(sfw_cond, np.where(is_sea, 0.0, ds["rws_oppwater_cond"].values)) + + ghb_rec = ghb.stress_period_data.get_data(0) + chd_rec = chd.stress_period_data.get_data(0) + ghb_cells = {cellid[1] for cellid in ghb_rec["cellid"]} + chd_cells = {cellid[1] for cellid in chd_rec["cellid"]} + + assert ghb_cells == {3, 4} + # Cell 1 is sea but has no active layer 0; CHD is a layer-0-only package, so it + # drops out while cell 0 remains. + assert chd_cells == {0} + # A cell must never be both a general head boundary and a fixed head. + assert ghb_cells.isdisjoint(chd_cells) + + # GHB goes into the first active layer: 0 for cell 3, 1 for the pinched cell 4. + assert {(cellid[1], cellid[0]) for cellid in ghb_rec["cellid"]} == {(3, 0), (4, 1)} + np.testing.assert_array_equal(ghb_rec["bhead"], [-1.0, -2.0]) + np.testing.assert_array_equal(ghb_rec["cond"], [250.0, 500.0]) + # The GHB water bodies are fresh. + np.testing.assert_array_equal(ghb_rec["CONCENTRATION"], [0.0, 0.0]) + + np.testing.assert_array_equal(chd_rec["head"], [sea_stage]) + # The sea is fixed at 18000 mg Cl-/l in the transport model (SEA_CHLORIDE_MG_L). + np.testing.assert_array_equal(chd_rec["CONCENTRATION"], [18000.0]) + assert ts_sea is None + + +def test_chd_time_series_branch_wiring(tmp_path): + """A list sea_stage wires a linear time series named ``sea_stage`` into the CHD.""" + ds, gwf = _sea_and_lake_ds(tmp_path) + series = [(0.0, 0.25), (365.0, 0.5)] + + _ghb, chd, _ts_sea = chd_ghb_from_major_surface_waters(ds, gwf, sea_stage=series) + + chd_rec = chd.stress_period_data.get_data(0) + # One record per active sea cell, each head the literal time-series name, not a number. + n_sea = int(ds["northsea"].sum()) + assert list(chd_rec["head"]) == ["sea_stage"] * n_sea + assert chd.ts.time_series_namerecord.get_data().tolist() == [("sea_stage",)] + assert chd.ts.interpolation_methodrecord.get_data().tolist() == [("linear",)] + np.testing.assert_array_equal( + np.array(chd.ts.timeseries.get_data().tolist()), + np.array(series), + ) + + +def test_chd_none_guard_returns_three_tuple(tmp_path): + """With no active sea cell the CHD is absent and the call still returns a triple.""" + # Both sea cells lose their layer 0, so nlmod.gwf.chd returns None. + ds, gwf = _sea_and_lake_ds(tmp_path, pinch_layer0=(0, 1)) + + result = chd_ghb_from_major_surface_waters(ds, gwf, sea_stage=[(0.0, 0.25), (365.0, 0.5)]) + + # Unpacking enforces the documented (ghb, chd, ts_sea) arity. + ghb, chd, ts_sea = result + assert chd is None + assert ts_sea is None + # The GHB is unaffected by the missing sea. + assert {cellid[1] for cellid in ghb.stress_period_data.get_data(0)["cellid"]} == {3, 4} diff --git a/tests/test_mf6_smoke.py b/tests/test_mf6_smoke.py new file mode 100644 index 0000000..7540910 --- /dev/null +++ b/tests/test_mf6_smoke.py @@ -0,0 +1,114 @@ +"""End-to-end MODFLOW 6 run behind ``postprocessing``'s output loaders. + +Every other test in the suite monkeypatches ``nlmod.gwf.output`` and hands +``check_budget_discrepancy`` a synthetic listing. This module is the only one that +writes real ``.lst``/``.hds``/``.grb`` files with the MODFLOW 6 executable and reads them +back through the production code path, so a change in nlmod/flopy output handling that +the mocked tests cannot see shows up here. +""" + +import nlmod +import pytest +import xarray as xr + +from nhflotools.postprocessing import add_output_to_ds, check_budget_discrepancy +from tests.util import make_structured_ds + + +pytestmark = pytest.mark.mf6 + + +@pytest.fixture(scope="session") +def mf6_exe(): + """Full path of the MODFLOW 6 executable, downloading it once if it is missing. + + ``nlmod.util.get_exe_path`` searches nlmod's own ``bin`` directory and then flopy's + metadata, so a previously downloaded executable is reused and the download happens at + most once per machine -- including on a fresh CI runner, which is why this test needs + no binaries to be installed beforehand. + """ + return nlmod.util.get_exe_path(exe_name="mf6", download_if_not_found=True) + + +# A 3x3 x 2-layer grid of 100 m cells on extent [0, 300] x [0, 300]. Cell centres are at +# x, y in {50, 150, 250}, so (150, 150) is the one cell not on the perimeter. +_NX = _NY = 3 +_NLAY = 2 +_DELR = 100.0 +_CENTRE = 150.0 +# h = (2*150 + 150) / 50, which is also the mean of the four neighbours: (5 + 13 + 7 + 11) / 4. +_CENTRE_HEAD = 9.0 + + +def _analytic_head(ds): + """Prescribed head field ``h = (2x + y) / 50`` [m NAP] on the cell centres. + + A field linear in x and y is harmonic, so it is the exact steady-state solution of + the finite-volume equations once the whole perimeter is held at its own value. The + coefficients are deliberately unequal in x and y: a swapped or mirrored axis in the + ``.grb``/``.hds`` reading path changes the field, whereas a symmetric field would + survive it. With ``x`` in {50, 150, 250} and ``y`` in {250, 150, 50} every head is a + small integer, so the comparison needs no tolerance for representation error. + """ + return (2.0 * ds["x"] + ds["y"]) / 50.0 + + +def test_tiny_run_budget_and_output_loading(tmp_path, mf6_exe): + """A real MF6 run closes its budget and reloads as the analytic head field. + + Builds the smallest model that still has an unknown: 3x3 cells x 2 layers, the eight + perimeter cells held constant at ``h = (2x + y) / 50`` in both layers, one free + interior column. Then, on the files MODFLOW actually wrote, + :func:`check_budget_discrepancy` must accept the listing and + :func:`add_output_to_ds` must return the analytic field. + """ + ws = str(tmp_path) + model_name = "smoke" + ds = make_structured_ds( + extent=(0.0, _NX * _DELR, 0.0, _NY * _DELR), + delr=_DELR, + top=0.0, + botm=[-10.0, -20.0], + kh=10.0, + kv=1.0, + model_name=model_name, + model_ws=ws, + ) + ds = nlmod.time.set_ds_time(ds, start="2022-01-01", time=[1.0], steady=True) + + analytic = _analytic_head(ds) + ds["chd_head"] = analytic.broadcast_like(ds["botm"]).transpose("layer", "y", "x") + # Everything but the single centre cell is a constant-head cell. + perimeter = (ds["x"] != _CENTRE) | (ds["y"] != _CENTRE) + ds["chd_mask"] = perimeter.broadcast_like(ds["botm"]).transpose("layer", "y", "x").astype(int) + + sim = nlmod.sim.sim(ds, exe_name=mf6_exe) + nlmod.sim.tdis(ds, sim) + # Tighter than the default so the residual is far below the 1e-8 comparison below. + nlmod.sim.ims(sim, complexity="SIMPLE", outer_dvclose=1e-9, inner_dvclose=1e-10) + gwf = nlmod.gwf.gwf(ds, sim) + nlmod.gwf.dis(ds, gwf) + nlmod.gwf.npf(ds, gwf) + # Start far from the answer so the interior column is genuinely solved for. + nlmod.gwf.ic(ds, gwf, starting_head=0.0) + nlmod.gwf.chd(ds, gwf) + nlmod.gwf.oc(ds, gwf) + nlmod.sim.write_and_run(sim, ds, write_ds=False, silent=True) + + # Only constant-head cells exchange water, so the volumetric budget must close; a + # wrong budget key or an inverted threshold turns this into a RuntimeError. + check_budget_discrepancy(ws, model_name, transport=False) + + ds, ctop = add_output_to_ds(ds, ws, model_name, transport=False) + assert ctop is None + assert "concentration" not in ds + + # The constant-head cells must come back at the head they were given, and the free + # centre cell at the mean of its four neighbours. ``expected`` is aligned on the x/y + # coordinates, so a mirrored or swapped axis in the .grb/.hds path is a mismatch. + # atol accommodates the IMS iterative solver, which stops at outer_dvclose rather + # than at algebraic exactness. + expected = analytic.broadcast_like(ds["head"]) + assert float(expected.isel(time=0, layer=0).sel(x=_CENTRE, y=_CENTRE)) == _CENTRE_HEAD + xr.testing.assert_allclose(ds["head"], expected, rtol=0.0, atol=1e-8) + assert ds["head"].sizes == {"time": 1, "layer": _NLAY, "y": _NY, "x": _NX} diff --git a/tests/test_nhflo_utils.py b/tests/test_nhflo_utils.py deleted file mode 100644 index e69de29..0000000 diff --git a/tests/test_nhflodata_contract.py b/tests/test_nhflodata_contract.py new file mode 100644 index 0000000..a0397a5 --- /dev/null +++ b/tests/test_nhflodata_contract.py @@ -0,0 +1,113 @@ +"""Contract tests: every nhflodata file the 09pwnmodel2 closure reads must exist. + +``nhflodata.get_paths.get_abs_data_path`` only *warns* when a dataset path is missing, so a +restructure of the NHFLO/data repository would otherwise surface as a ``FileNotFoundError`` +deep inside a model run. These tests stat every (dataset, relative file) pair that nhflotools +-- or ``modelscripts/09pwnmodel2/01_pwnmodel2.py`` -- hardcodes, against the mockup data +packaged with nhflodata (the autouse hygiene fixture deletes ``NHFLODATA_LOCATION``, so +resolution always lands on the mockup). +""" + +from pathlib import Path + +import pandas as pd +import pytest + +from nhflotools.pwnlayers3.layers import layer_names + +get_paths = pytest.importorskip("nhflodata.get_paths") +get_abs_data_path = get_paths.get_abs_data_path + +# Hydrogeological units of pwnlayers3/layers.py:get_gdf_boundaries. +_BOUNDARY_UNITS = ("11", "12", "13", "21", "22", "31", "32") + +# Aquitards whose kh comes from NHDZ transmissivity points (pwnlayers3/layers.py:get_kh). +_KD_NHDZ_LAYERS = ("S12", "S13", "S21", "S22", "S31") + +_KOPPELTABEL_CSV = "bodemlagenvertaaltabelv2.csv" + + +def _bodemlagen_2024_files(): + """Yield the relative paths read from the ``bodemlagen_pwn_2024`` dataset. + + Derived from the path expressions in ``nhflotools/pwnlayers3/layers.py`` and + ``nhflotools/pwnlayers3/plot.py`` rather than from a directory listing, so a file the + code needs but the dataset stopped shipping is a failure, not a silent omission. + + Yields + ------ + str + Path relative to the dataset root, POSIX-style. + """ + yield "botm/botm.geojson" # get_botm + plot._load_and_project_source_botm + yield "boundaries/triwaco_model_nhdz.geojson" # get_kh + for unit in _BOUNDARY_UNITS: # get_gdf_boundaries + yield f"boundaries/S{unit}/S{unit}.geojson" + for name in layer_names: + # get_kh reads K_combined; aquitards read C_combined (get_kh + get_kv) + prefix = "K" if name.startswith("W") else "C" + yield f"conductances/{prefix}{name}_combined.geojson" + for name in _KD_NHDZ_LAYERS: # get_kh, NHDZ branch + yield f"conductances/KD{name}_NHDZ.geojson" + + +def _data_files(): + """Yield ``(dataset_name, relative_path)`` for every hardcoded read. + + Yields + ------ + tuple of (str, str) + Dataset name as passed to ``get_abs_data_path`` and the relative file joined onto it. + """ + for relative_path in _bodemlagen_2024_files(): + yield "bodemlagen_pwn_2024", relative_path + # panden.py reads the .shp; geopandas additionally requires the .dbf/.shx sidecars and + # needs the .prj to attach EPSG:28992 before the grid intersection. + for suffix in ("shp", "dbf", "shx", "prj"): + yield "oppervlaktewater_pwn_shapes_panden", f"Panden_ICAS_IKIEF.{suffix}" + yield "nhi_chloride_concentration", "chloride_p50.nc" # nhi_chloride.py:32 + yield "wells_pwn", "pumping_infiltration_wells.geojson" # well.py:36 + yield "wells_pwn", "sec_flows.feather" # well.py:46 + yield "wells_tata", "tata_zoutwaterbronnen.geojson" # well.py:108 + yield "wells_tata", "tata_zoetwaterbronnen.geojson" # well.py:120 + # Read by 01_pwnmodel2.py itself (lines 126, 363, 400). + yield "lakes_pwn", "lakes_pwn.geojson" + yield "drains_pwn", "drains_pwn.geojson" + yield "hfb_pwn", "hfb_pwn.geojson" + + +_DATA_FILES = tuple(_data_files()) + + +@pytest.mark.parametrize(("dataset", "relative_path"), _DATA_FILES, ids=[f"{d}:{r}" for d, r in _DATA_FILES]) +def test_mockup_data_file_resolves_and_is_nonempty(dataset, relative_path): + """The resolved path exists and holds bytes. + + ``get_abs_data_path`` warns instead of raising on a missing dataset, and a zero-byte + placeholder would satisfy a bare existence check while breaking every reader, so both + are asserted. + """ + root = Path(get_abs_data_path(name=dataset, version="latest", location="get_from_env")) + path = root / relative_path + + assert path.is_file(), f"{dataset} is missing {relative_path} (resolved to {path})" + assert path.stat().st_size > 0, f"{dataset}/{relative_path} is empty" + + +def test_koppeltabel_columns_and_layer_coverage(): + """The koppeltabel exposes the columns and the layer names the merge indexes by. + + ``merge_layer_models.combine_two_layer_models`` defaults to the column names + ``'Regis II v2.2'`` / ``'ASSUMPTION1'`` and then does + ``layer_model_other.sel(layer=)``. So the non-null ASSUMPTION1 + entries must be exactly the pwnlayers3 layer set: an extra name raises inside ``sel``, + a missing name silently drops that PWN layer from the merged model. Rows without a + REGIS name can never be coupled, so that column must be complete. + """ + root = Path(get_abs_data_path(name="bodemlagen_pwn_regis_koppeltabel", version="latest", location="get_from_env")) + # Same read as pwnlayers3/layers.py:166. + df = pd.read_csv(root / _KOPPELTABEL_CSV, skiprows=0, index_col=0) + + assert {"Regis II v2.2", "ASSUMPTION1"}.issubset(df.columns) + assert set(df["ASSUMPTION1"].dropna()) == set(layer_names) + assert df["Regis II v2.2"].notna().all() diff --git a/tests/test_nhi_chloride.py b/tests/test_nhi_chloride.py new file mode 100644 index 0000000..7b11efa --- /dev/null +++ b/tests/test_nhi_chloride.py @@ -0,0 +1,178 @@ +"""Tests for nhflotools.nhi_chloride.get_nhi_chloride_concentration. + +The source data is a tiny hand-written ``chloride_p50.nc`` with descending y (as the real +NHI file has), so the nearest-neighbour interpolation onto the model grid is exact and the +thickness-weighted mean per model layer has an integer-valued known answer. +""" + +import numpy as np +import pytest +import xarray as xr + +from nhflotools.nhi_chloride import SEA_CHLORIDE_MG_L, get_nhi_chloride_concentration +from tests.util import make_rect_vertex_ds + +# Source voxels of the default fixture: [0, -4] m at 100 mg/l over [-4, -8] m at 300 mg/l. +_VOXEL_TOP = (0.0, -4.0) +_VOXEL_BOTTOM = (-4.0, -8.0) +_VOXEL_C = (100.0, 300.0) + +# Cell centres of the default 2x2 vertex grid (extent 0-200 m, row-major from the NW corner). +_CELL_X = (50.0, 150.0) +_CELL_Y = (150.0, 50.0) + + +def _write_chloride_nc(directory, values, top=_VOXEL_TOP, bottom=_VOXEL_BOTTOM, x=_CELL_X, y=_CELL_Y): + """Write a minimal ``chloride_p50.nc`` and return the directory holding it. + + Parameters + ---------- + directory : pathlib.Path + Directory the file is written into. + values : array_like + Either one concentration per source layer (constant in x and y) or a full + ``(layer, y, x)`` array. + top, bottom : sequence of float + Voxel top/bottom elevations [mNAP], one per source layer. + x, y : sequence of float + Source cell centres; ``y`` is descending, matching the real NHI file. + + Returns + ------- + str + Path of the directory, ready to pass as ``data_path_nhi_chloride``. + """ + values = np.asarray(values, dtype=float) + if values.ndim == 1: + values = np.tile(values[:, None, None], (1, len(y), len(x))) + da = xr.DataArray( + values, + dims=("layer", "y", "x"), + coords={ + "layer": np.arange(len(top)), + "y": list(y), + "x": list(x), + "top": ("layer", list(top)), + "bottom": ("layer", list(bottom)), + }, + ) + xr.Dataset({"chloride_p50": da}).to_netcdf(directory / "chloride_p50.nc") + return str(directory) + + +def _model_ds(top=0.0, botm=(-4.0, -8.0), northsea=None): + """Build a 2x2 vertex model dataset with the ``northsea`` variable the function requires.""" + ds = make_rect_vertex_ds(top=top, botm=botm) + n = ds.sizes["icell2d"] + ds["northsea"] = ("icell2d", np.zeros(n, dtype=int) if northsea is None else np.asarray(northsea, dtype=int)) + return ds + + +@pytest.mark.parametrize( + ("top", "botm", "expected"), + [ + # One model layer spanning both voxels: (4*100 + 4*300) / 8 = 200. + (0.0, (-8.0,), (200.0,)), + # Model layers coincide with the voxels, so each layer reproduces its voxel exactly. + (0.0, (-4.0, -8.0), (100.0, 300.0)), + # Model layer [-1, -6] clips both voxels: (3*100 + 2*300) / 5 = 180. + (-1.0, (-6.0,), (180.0,)), + ], + ids=["spans-both-voxels", "voxel-aligned", "clips-both-voxels"], +) +def test_thickness_weighted_mean_per_model_layer(tmp_path, top, botm, expected): + """The per-layer value is the voxel mean weighted by the clipped overlap thickness.""" + path = _write_chloride_nc(tmp_path, _VOXEL_C) + ds = _model_ds(top=top, botm=botm) + + da = get_nhi_chloride_concentration(ds, path) + + assert da.dims == ("layer", "icell2d") + assert da.attrs["units"] == "mg/l" + # The source is uniform in x and y, so every cell of a layer carries the same value. + np.testing.assert_array_equal(da.values, np.repeat(np.array(expected)[:, None], ds.sizes["icell2d"], axis=1)) + + +def test_nan_voxel_leaves_denominator_untouched(tmp_path): + """A NaN voxel drops out of numerator *and* denominator, so the answer is undiluted.""" + path = _write_chloride_nc(tmp_path, [np.nan, 300.0]) + ds = _model_ds(botm=(-8.0,)) + + da = get_nhi_chloride_concentration(ds, path) + + # Both voxels are 4 m thick; keeping the NaN voxel's thickness in the denominator + # would halve the result to 150.0. + np.testing.assert_array_equal(da.values, np.full((1, 4), 300.0)) + + +def test_sea_override_applies_to_layer_zero_only(tmp_path): + """``northsea == 1`` forces the sea concentration in layer 0 and nowhere else.""" + path = _write_chloride_nc(tmp_path, _VOXEL_C) + ds = _model_ds(northsea=[1, 0, 0, 0]) + + da = get_nhi_chloride_concentration(ds, path) + + expected = np.array([ + [SEA_CHLORIDE_MG_L, 100.0, 100.0, 100.0], + [300.0, 300.0, 300.0, 300.0], + ]) + np.testing.assert_array_equal(da.values, expected) + + +def test_model_layers_outside_source_stack_are_filled_vertically(tmp_path): + """Layers above and below the voxel stack are rescued by the bfill/ffill pass.""" + path = _write_chloride_nc(tmp_path, _VOXEL_C) + # Layer 0 = [4, 0] sits entirely above the voxels and layer 3 = [-8, -12] entirely + # below, so both aggregate to 0/0 = NaN before the fill. + ds = _model_ds(top=4.0, botm=(0.0, -4.0, -8.0, -12.0)) + + da = get_nhi_chloride_concentration(ds, path) + + assert not da.isnull().any() + # bfill carries layer 1 up into layer 0; ffill carries layer 2 down into layer 3. + expected = np.repeat(np.array([100.0, 100.0, 300.0, 300.0])[:, None], 4, axis=1) + np.testing.assert_array_equal(da.values, expected) + + +def test_cells_outside_source_extent_take_the_nearest_neighbour(tmp_path): + """Uncovered cells are filled horizontally from the nearest covered cell, not extrapolated.""" + # The source only covers y = 150 and y = 100, and its two rows differ, so the southern + # model cells (y = 50) fall outside its extent and interpolate to NaN. + values = np.array([ + [[100.0, 900.0], [500.0, 700.0]], + [[300.0, 1100.0], [600.0, 800.0]], + ]) + path = _write_chloride_nc(tmp_path, values, y=(150.0, 100.0)) + ds = _model_ds() + + da = get_nhi_chloride_concentration(ds, path) + + # Cell 2 (50, 50) is 100 m from cell 0 and 141 m from cell 1, so it inherits cell 0; + # cell 3 likewise inherits cell 1. Extrapolating the source instead would yield the + # y = 100 row, i.e. 500/700 and 600/800. + expected = np.array([ + [100.0, 900.0, 100.0, 900.0], + [300.0, 1100.0, 300.0, 1100.0], + ]) + np.testing.assert_array_equal(da.values, expected) + + +def test_cached_call_reproduces_the_uncached_result(tmp_path): + """The cache_netcdf round trip returns the same values and attrs as a direct call.""" + path = _write_chloride_nc(tmp_path, _VOXEL_C) + cachedir = tmp_path / "cache" + cachedir.mkdir() + + reference = get_nhi_chloride_concentration(_model_ds(northsea=[1, 0, 0, 0]), path) + first = get_nhi_chloride_concentration( + _model_ds(northsea=[1, 0, 0, 0]), path, cachedir=str(cachedir), cachename="chloride" + ) + second = get_nhi_chloride_concentration( + _model_ds(northsea=[1, 0, 0, 0]), path, cachedir=str(cachedir), cachename="chloride" + ) + + assert (cachedir / "chloride.nc").exists() + assert (cachedir / "chloride.pklz").exists() + np.testing.assert_array_equal(first.values, reference.values) + np.testing.assert_array_equal(second.values, reference.values) + assert second.attrs["units"] == "mg/l" diff --git a/tests/test_panden.py b/tests/test_panden.py index e69de29..0792e16 100644 --- a/tests/test_panden.py +++ b/tests/test_panden.py @@ -0,0 +1,184 @@ +"""Tests for :mod:`nhflotools.panden` (the PWN infiltration-pond RIV package).""" + +import logging + +import geopandas as gpd +import numpy as np +import pytest +from shapely.geometry import box + +from nhflotools.panden import get_oppervlakte_pwn_shapes, riv_from_oppervlakte_pwn +from tests.util import make_gwf_disv, make_rect_vertex_ds + +# Stages hardcoded in panden.py; rbot is stage - 2.0 m. +ICAS_STAGE = 2.8 +IKIEF_STAGE = 5.8 +DEPTH = 2.0 + + +def _write_panden_shp(directory, geometries, names): + """Write the real ``Panden_ICAS_IKIEF.shp`` the reader expects. + + A genuine shapefile (rather than a stubbed GeoDataFrame) is written so the + ``gpd.read_file`` path, the DBF round-trip of ``Naam`` and ``make_valid`` are all + exercised. + + Parameters + ---------- + directory : pathlib.Path + Directory that becomes ``data_path_panden``. + geometries : sequence of shapely.geometry.Polygon + One polygon per pand. + names : sequence of str + Value of the ``Naam`` attribute of each pand. + + Returns + ------- + str + ``directory`` as a string, ready to pass as ``data_path_panden``. + """ + gdf = gpd.GeoDataFrame({"Naam": list(names)}, geometry=list(geometries), crs="EPSG:28992") + gdf.to_file(directory / "Panden_ICAS_IKIEF.shp") + return str(directory) + + +@pytest.fixture +def panden_ds_gwf(tmp_path): + """``(ds, gwf)`` on the default 2x2 100 m grid, transport enabled. + + Cell 0 spans x in [0, 100], y in [100, 200]; cell 1 spans x in [100, 200], y in + [100, 200]. Transport is on so the SSM registration branch is reachable. + """ + ds = make_rect_vertex_ds(transport=1) + return make_gwf_disv(ds, tmp_path / "model") + + +def test_get_oppervlakte_pwn_shapes_assigns_stages_and_drops_other_names(tmp_path, caplog): + """Only ICAS/IKIEF panden survive, with their stage, resistance and rbot set. + + A "VIJVER" pond carries neither name: it must be dropped with a warning rather than + reaching the RIV package with a NaN stage (regression, see #41). + """ + path = _write_panden_shp( + tmp_path, + [box(0, 0, 10, 10), box(20, 0, 30, 10), box(40, 0, 50, 10)], + ["ICAS-noord", "IKIEF-3", "VIJVER"], + ) + + with caplog.at_level(logging.WARNING, logger="nhflotools.panden"): + shapes = get_oppervlakte_pwn_shapes(data_path_panden=path) + + assert list(shapes["Naam"]) == ["ICAS-noord", "IKIEF-3"] + np.testing.assert_array_equal(shapes["stage"].to_numpy(), [ICAS_STAGE, IKIEF_STAGE]) + np.testing.assert_array_equal(shapes["c"].to_numpy(), [1.0, 1.0]) + # rbot is 2 m below the stage of that same pand. + np.testing.assert_array_equal(shapes["rbot"].to_numpy(), [ICAS_STAGE - DEPTH, IKIEF_STAGE - DEPTH]) + assert any("VIJVER" in rec.getMessage() for rec in caplog.records) + + +def test_riv_single_cell_known_answer(tmp_path, panden_ds_gwf): + """One 50x50 m pand inside one cell yields one fully derived RIV record.""" + ds, gwf = panden_ds_gwf + # 50 m x 50 m, wholly inside cell 0 (x in [0, 100], y in [100, 200]). + path = _write_panden_shp(tmp_path, [box(10, 110, 60, 160)], ["ICAS-noord"]) + + riv = riv_from_oppervlakte_pwn(ds, gwf, data_path_panden=path) + spd = riv.stress_period_data.get_data(0) + + # cond = intersected area / c = (50 * 50) / 1.0. rbot = 2.8 - 2.0 is above botm[0] + # = -10, so lay_of_rbot puts the record in layer 0. The aux (CONCENTRATION) value + # precedes the boundname, which is the order flopy reads them in. + assert len(spd) == 1 + cellid, stage, cond, rbot, aux, boundname = tuple(spd[0]) + assert tuple(cellid) == (0, 0) + assert stage == ICAS_STAGE + assert cond == 50.0 * 50.0 + assert rbot == ICAS_STAGE - DEPTH + assert aux == 0.0 + assert boundname.strip() == "ICAS-noord" + + +def test_riv_conductance_is_conserved_and_stage_area_weighted(tmp_path, panden_ds_gwf): + """Splitting a pand over cells partitions its conductance; a shared cell blends stages. + + The pre-fix code aggregated the resistance ``c`` instead of ``area / c``, a silent + ~1e4x conductance error that this partition identity pins down. + """ + ds, gwf = panden_ds_gwf + path = _write_panden_shp( + tmp_path, + # Straddles the cell 0 / cell 1 boundary at x = 100: 25 m x 50 m in cell 0 and + # 75 m x 50 m in cell 1. The second pand sits wholly in cell 0. + [box(75, 110, 175, 160), box(10, 110, 60, 160)], + ["ICAS-noord", "IKIEF-3"], + ) + + riv = riv_from_oppervlakte_pwn(ds, gwf, data_path_panden=path) + records = {tuple(rec[0]): rec for rec in riv.stress_period_data.get_data(0)} + + icas_cell0 = 25.0 * 50.0 + icas_cell1 = 75.0 * 50.0 + ikief_cell0 = 50.0 * 50.0 + assert set(records) == {(0, 0), (0, 1)} + + # c == 1.0 everywhere, so cond per cell is exactly the intersected area, and the + # two pieces of the straddling pand must add back up to its full area. + assert records[(0, 0)][2] == icas_cell0 + ikief_cell0 + assert records[(0, 1)][2] == icas_cell1 + total_cond = records[(0, 0)][2] + records[(0, 1)][2] + assert total_cond == icas_cell0 + icas_cell1 + ikief_cell0 + + # Cell 0 holds both panden: area-weighted stage, minimum rbot. + expected_stage = (ICAS_STAGE * icas_cell0 + IKIEF_STAGE * ikief_cell0) / (icas_cell0 + ikief_cell0) + # rel tolerance accommodates the pandas groupby accumulation order of the + # area * stage products, which need not match the order used here. + assert records[(0, 0)][1] == pytest.approx(expected_stage, rel=1e-12) + assert records[(0, 0)][3] == ICAS_STAGE - DEPTH + assert records[(0, 1)][1] == ICAS_STAGE + assert records[(0, 1)][3] == ICAS_STAGE - DEPTH + + +def test_riv_rbot_selects_first_layer_it_reaches_into(tmp_path): + """A pand bottom below the first layer bottom is placed in the layer that holds it.""" + # Layer 0 spans [5, 1], layer 1 spans [1, -10]; rbot = 0.8 lies inside layer 1. + ds = make_rect_vertex_ds(top=5.0, botm=(1.0, -10.0)) + ds, gwf = make_gwf_disv(ds, tmp_path / "model") + path = _write_panden_shp(tmp_path, [box(10, 110, 60, 160)], ["ICAS-noord"]) + + riv = riv_from_oppervlakte_pwn(ds, gwf, data_path_panden=path) + spd = riv.stress_period_data.get_data(0) + + assert len(spd) == 1 + assert tuple(spd[0][0]) == (1, 0) + + +def test_riv_outside_extent_returns_none_and_leaves_ds_untouched(tmp_path, panden_ds_gwf): + """No intersection means no RIV package and no SSM source registered.""" + ds, gwf = panden_ds_gwf + # Model extent is [0, 200] x [0, 200]. + path = _write_panden_shp(tmp_path, [box(1000, 1000, 1100, 1100)], ["ICAS-noord"]) + + assert riv_from_oppervlakte_pwn(ds, gwf, data_path_panden=path) is None + assert "ssm_sources" not in ds.attrs + assert not [p for p in gwf.get_package_list() if p.lower().startswith("riv")] + + +@pytest.mark.parametrize("transport", [0, 1]) +def test_riv_ssm_registration_is_transport_gated_and_idempotent(tmp_path, transport): + """The RIV aux is registered as an SSM source once, and only under transport. + + Rebuilding the flow model on an existing dataset (the model script re-run against a + cached ``ds``) must not append the same package a second time. + """ + ds = make_rect_vertex_ds(transport=transport) + path = _write_panden_shp(tmp_path, [box(10, 110, 60, 160)], ["ICAS-noord"]) + + ds, gwf = make_gwf_disv(ds, tmp_path / "run1") + riv = riv_from_oppervlakte_pwn(ds, gwf, data_path_panden=path) + ds, gwf2 = make_gwf_disv(ds, tmp_path / "run2") + riv_from_oppervlakte_pwn(ds, gwf2, data_path_panden=path) + + if transport: + assert ds.attrs["ssm_sources"] == [riv.package_name] + else: + assert "ssm_sources" not in ds.attrs diff --git a/tests/test_polder.py b/tests/test_polder.py new file mode 100644 index 0000000..25e4e1e --- /dev/null +++ b/tests/test_polder.py @@ -0,0 +1,136 @@ +"""Tests for nhflotools.polder.drn_from_waterboard_data. + +Only the download is faked (``nlmod.read.waterboard.download_data`` hits a live, uncached +and flaky HHNK ArcGIS REST service at polder.py:40); the grid intersection, the +area-weighted aggregation, the nearest-neighbour stage fill and the flopy DRN package are +all built for real. + +Note on conductance: polder.py:59 deliberately reads the FULL cell area, not the +intersected area (issue #51), so nothing here asserts an intersected-area conductance. +""" + +import numpy as np +import xarray as xr +from shapely.geometry import box + +from nhflotools.polder import drn_from_waterboard_data +from tests.util import cell_polygon, make_gdf, make_gwf_disv, make_rect_vertex_ds + +# make_rect_vertex_ds() default: 2x2 cells of 100 m, row-major from the north-west corner, +# so cell 0=(50,150), 1=(150,150), 2=(50,50), 3=(150,50) and every cell area is 100*100. +CELL_AREA = 100.0 * 100.0 + + +def _patch_download(monkeypatch, gdf): + """Replace the live waterboard download with ``gdf`` and record the call kwargs.""" + calls = {} + + def _fake(**kwargs): + calls.update(kwargs) + return gdf + + monkeypatch.setattr("nlmod.read.waterboard.download_data", _fake) + return calls + + +def _reclist(drn): + """Return the DRN stress period data as sorted ``(layer, icell2d, elev, cond)`` tuples.""" + return sorted((int(r[0][0]), int(r[0][1]), float(r[1]), float(r[2])) for r in drn.stress_period_data.get_data(0)) + + +def test_drn_partition_over_polder_land_and_sea(tmp_path, monkeypatch): + """Polder cells, uncovered land, sea cells and blocked layers each land in the right bin. + + Regression for the zero-initialised conductance array, which gave uncovered land + zero-conductance (i.e. inert) drains and simultaneously drained the North Sea. + + ``cbot`` is 2.0, not 1.0: dividing the cell area by 1.0 is an identity, so a + conductance that ignored ``cbot`` altogether would pass unnoticed. + """ + cbot = 2.0 + ds = make_rect_vertex_ds() + # Zero-thickness layer 0 in cell 3 -> nlmod derives idomain 0 there, first active layer 1. + ds["botm"][0, 3] = float(ds["top"][3]) + ds, gwf = make_gwf_disv(ds, tmp_path) + ds["ahn"] = xr.DataArray(np.array([0.0, 4.25, 1.0, 2.5]), dims="icell2d") + ds["northsea"] = xr.DataArray(np.array([0, 0, 1, 0]), dims="icell2d") + + gdf = make_gdf([cell_polygon(ds, 0)], summer_stage=1.0, winter_stage=3.0) + _patch_download(monkeypatch, gdf) + + drn = drn_from_waterboard_data(ds, gwf, cbot=cbot) + + # cell 0 covered -> mean(1.0, 3.0); cells 1-3 uncovered -> ahn. + np.testing.assert_array_equal(ds["drn_elev"].values, [2.0, 4.25, 1.0, 2.5]) + # Every non-sea cell gets area/cbot, whether or not a polder polygon covers it. + np.testing.assert_array_equal(ds["drn_cond"].values, [CELL_AREA / cbot, CELL_AREA / cbot, np.nan, CELL_AREA / cbot]) + # The historical failure mode was a *zero* fallback conductance, not a NaN one. + assert ds["drn_cond"].values[1] > 0.0 + + # The sea cell is dropped; the blocked cell drains from its first active layer. + assert _reclist(drn) == [ + (0, 0, 2.0, CELL_AREA / cbot), + (0, 1, 4.25, CELL_AREA / cbot), + (1, 3, 2.5, CELL_AREA / cbot), + ] + + +def test_drn_elev_skips_nan_stages_and_fills_all_nan_cells_from_nearest(gwf_disv, monkeypatch): + """A missing stage is skipped in the mean; a cell missing both is filled by its nearest peer.""" + ds, gwf = gwf_disv + # ahn differs from every stage, so an elevation equal to a stage proves the stage path won. + ds["ahn"] = xr.DataArray(np.array([9.0, 7.0, 8.0, 5.0]), dims="icell2d") + + gdf = make_gdf( + [cell_polygon(ds, 0), cell_polygon(ds, 1), cell_polygon(ds, 2)], + summer_stage=[1.0, np.nan, 2.0], + winter_stage=[np.nan, np.nan, 6.0], + ) + _patch_download(monkeypatch, gdf) + + drn_from_waterboard_data(ds, gwf) + + # cell 0: one-sided mean of (1.0, NaN) is 1.0, not 0.5. + # cell 1 (150,150): both stages NaN -> nearest valid donor is cell 0 at (50,150), 100 m + # away, versus cell 2 at (50,50), 100*sqrt(2) m away; so 1.0, not 4.0 and not ahn. + # cell 2: mean(2.0, 6.0) = 4.0. cell 3: no polygon -> ahn. + np.testing.assert_array_equal(ds["drn_elev"].values, [1.0, 1.0, 4.0, 5.0]) + + +def test_drn_empty_download_returns_none_and_leaves_ds_untouched(gwf_disv, monkeypatch): + """No level areas in the extent short-circuits before any variable is written to ds.""" + ds, gwf = gwf_disv + ds["ahn"] = xr.DataArray(np.zeros(4), dims="icell2d") + calls = _patch_download(monkeypatch, make_gdf([], summer_stage=[], winter_stage=[])) + + assert drn_from_waterboard_data(ds, gwf) is None + assert "drn_elev" not in ds + assert "drn_cond" not in ds + # The download must be scoped to the model, not to the whole waterboard. + assert calls["extent"] == ds.extent + + +def test_drn_duplicate_download_index_is_deduplicated(gwf_disv, monkeypatch): + """Repeated level-area identifiers are made unique so gdf_to_grid accepts the frame. + + ``nlmod.grid.gdf_to_grid`` raises ``ValueError: gdf should not have duplicate columns + or index``, so reaching an area-weighted elevation at all proves the renaming ran. + """ + ds, gwf = gwf_disv + ds["ahn"] = xr.DataArray(np.zeros(4), dims="icell2d") + + # Two features share the id "A": one fills cell 0, one fills the west half of cell 1. + # "B" fills the east half of cell 1. + gdf = make_gdf( + [cell_polygon(ds, 0), box(100.0, 100.0, 150.0, 200.0), box(150.0, 100.0, 200.0, 200.0)], + summer_stage=[2.0, 6.0, 10.0], + winter_stage=[2.0, 6.0, 10.0], + ) + gdf.index = ["A", "A", "B"] + _patch_download(monkeypatch, gdf) + + drn_from_waterboard_data(ds, gwf) + + # cell 1 keeps both halves: (5000*6.0 + 5000*10.0) / 10000 = 8.0. Losing either + # duplicate would give 6.0 or 10.0 instead. + np.testing.assert_array_equal(ds["drn_elev"].values, [2.0, 8.0, 0.0, 0.0]) diff --git a/tests/test_postprocessing.py b/tests/test_postprocessing.py index 681077a..12589d0 100644 --- a/tests/test_postprocessing.py +++ b/tests/test_postprocessing.py @@ -2,17 +2,21 @@ The interface-elevation helper supersedes the old layer-discretised ``grensvlak`` with ``nlmod.dims.get_isosurface``. These tests pin the intended behaviour: it agrees with the old -logic at the fresh/salt limits (per column) and interpolates the crossing in between. +logic at the fresh/salt limits (per column) and interpolates the crossing in between. The +remaining tests cover the budget guard against a real MODFLOW listing, the derived fields +``add_output_to_ds`` builds, and the exact set of figures ``plot_result_maps`` writes. """ from unittest import mock import numpy as np import pandas as pd +import pytest import xarray as xr import nhflotools.postprocessing as pp from nhflotools.postprocessing import interface_elevation +from tests.util import make_rect_vertex_ds, write_mf6_listing # col0 crosses the fresh threshold, col1 stays fresh throughout, col2 stays salt throughout. _CONC = [[100, 500, 2000, 9000], [50, 80, 90, 120], [3000, 4000, 5000, 6000]] @@ -140,40 +144,214 @@ def test_interface_deepens_with_threshold(): assert depths[0] > depths[1] > depths[2] -def _run_budget_check(values, max_pct): - """Call check_budget_discrepancy with Mf6ListBudget stubbed to a known discrepancy column.""" - fake = mock.Mock() - fake.get_dataframes.return_value = [pd.DataFrame({"PERCENT_DISCREPANCY": values})] - with mock.patch.object(pp.flopy.utils, "Mf6ListBudget", return_value=fake): - pp.check_budget_discrepancy("ws", "m", transport=False, max_gwf_pct=max_pct) +# --- check_budget_discrepancy against a real MODFLOW listing ------------------------------- -def test_check_budget_discrepancy_uses_absolute_value(): - """A large negative discrepancy trips the guard (absolute value), a small one passes.""" - raised = False - try: - _run_budget_check([0.1, -5.0, 0.2], max_pct=1.0) - except RuntimeError: - raised = True - assert raised, "expected RuntimeError when |discrepancy| exceeds the threshold" +@pytest.mark.parametrize( + ("pct_disc", "max_gwf_pct", "expect_raise"), + [ + (2.0, 4.0, False), # comfortably below the threshold + (2.0, 2.0, True), # exactly at it: the guard is inclusive (>=) + (-2.0, 2.0, True), # a negative discrepancy is compared on its absolute value + ], +) +def test_check_budget_real_listing(tmp_path, pct_disc, max_gwf_pct, expect_raise): + """The guard trips at (not just above) the threshold when flopy really parses the listing. - _run_budget_check([0.1, -0.2, 0.3], max_pct=1.0) # below threshold -> must not raise + Values are integral so the float32 column flopy produces holds them exactly and the + ``>=`` boundary is unambiguous. + """ + write_mf6_listing(tmp_path / "test.lst", pct_disc) + def run(): + pp.check_budget_discrepancy(str(tmp_path), "test", transport=False, max_gwf_pct=max_gwf_pct) + + if expect_raise: + with pytest.raises(RuntimeError, match="discrepancy too large"): + run() + else: + run() + + +def test_check_budget_unparsable_listing(tmp_path): + """A listing flopy cannot parse raises rather than silently passing the guard.""" + (tmp_path / "test.lst").write_text("this is not a MODFLOW listing\n") + with pytest.raises(RuntimeError, match="Could not parse"): + pp.check_budget_discrepancy(str(tmp_path), "test", transport=False) + + +# --- add_output_to_ds ---------------------------------------------------------------------- + +# Three layers of thickness 4, 4 and 8 m below top=0: cell centres at -2, -6 and -12 m NAP and +# a thickness sum of 16, so every weighted mean below is exactly representable. +_BOTM = (-4.0, -8.0, -16.0) +_Z_CENTRE = np.array([-2.0, -6.0, -12.0]) +_THICKNESS = np.array([4.0, 4.0, 8.0]) +_THRESHOLD_BRAK = 8000.0 +# chloride 500/1500/14500 over thicknesses 4/4/8: (4*500 + 4*1500 + 8*14500) / 16 +_PROFILE = np.array([500.0, 1500.0, 14500.0]) +_MEAN_T0 = 7750.0 + + +def _output_ds(ntime=2): + """Two cells x three layers, with the density attrs ``freshwater_head`` reads from ds.""" + ds = make_rect_vertex_ds(nx=2, ny=1, botm=_BOTM, transport=1) + ds.attrs["drhodc"] = 0.5 + return ds, pd.to_datetime([f"2022-{m + 1:02d}-01" for m in range(ntime)]) + + +def _da(values, time, layer): + """Broadcast a per-layer sequence to a ``(time, layer, icell2d)`` DataArray of 2 cells.""" + arr = np.broadcast_to(np.asarray(values, float)[None, :, None], (len(time), len(layer), 2)) + return xr.DataArray(arr.copy(), dims=("time", "layer", "icell2d"), coords={"time": time, "layer": layer}) -def test_interface_time_dim_uses_2d_slices(): - """With a time dim, get_isosurface is called per time step (its numpy fallback rejects 3D input).""" - ds = _synthetic_ds() - base = _conc(_CONC) - conc = xr.concat([base, base], dim="time").assign_coords(time=[0, 1]) - seen_dims = [] - real = pp.nlmod.dims.get_isosurface - def spy(da, *args, **kwargs): - seen_dims.append(tuple(da.dims)) - return real(da, *args, **kwargs) +def _patch_loaders(monkeypatch, head, conc): + """Stub the three nlmod readers of the .hds/.ucn binaries (named seam, no MODFLOW run). - with mock.patch.object(pp.nlmod.dims, "get_isosurface", side_effect=spy): - pp.interface_elevation(ds, conc, _THRESHOLD) + Returns the list the transport readers append to, so a test can assert they stayed unused. + """ + calls = [] + monkeypatch.setattr(pp.nlmod.gwf.output, "get_heads_da", lambda *_a, **_kw: head) + + def _conc(*_a, **_kw): + calls.append("conc") + return conc - assert seen_dims, "get_isosurface was never called" - assert all("time" not in dims for dims in seen_dims), f"get_isosurface received a time dim: {seen_dims}" + def _ctop(c, **_kw): + calls.append("ctop") + return c.isel(layer=0) + + monkeypatch.setattr(pp.nlmod.gwt.output, "get_concentration_da", _conc) + monkeypatch.setattr(pp.nlmod.gwt.output, "get_concentration_at_gw_surface", _ctop) + return calls + + +def test_add_output_freshwater_head_and_head_fill(monkeypatch): + """The density correction is an identity at zero chloride and exact and signed below it. + + With denseref=1024 and drhodc=0.5, chloride 1024 gives density 1536, so the two factors + density/denseref = 1.5 and (density - denseref)/denseref = 0.5 are exact binary fractions. + """ + ds, time = _output_ds() + layer = ds.layer.values + # layer 0 is chloride-free (rho == rho_ref -> hf must equal the head itself) + conc = _da([0.0, 1024.0, 1024.0], time, layer) + head = _da([2.0, 2.0, 2.0], time, layer) + head[0, 0, 0] = np.nan # one dry cell, filled from the layer below by bfill + _patch_loaders(monkeypatch, head, conc) + + ds, _ = pp.add_output_to_ds(ds, "ws", "test", denseref=1024.0) + + # bfill must leave every finite head untouched and fill the planted gap from layer 1 + finite = np.isfinite(head.values) + np.testing.assert_array_equal(ds["head_filled"].values[finite], head.values[finite]) + assert ds["head_filled"].values[0, 0, 0] == head.values[0, 1, 0] + + # hf = rho/rho_ref * h - (rho - rho_ref)/rho_ref * z, with h = 2 everywhere after the fill + expected = np.array([2.0, 1.5 * 2.0 - 0.5 * _Z_CENTRE[1], 1.5 * 2.0 - 0.5 * _Z_CENTRE[2]]) + np.testing.assert_array_equal(ds["freshwater_head"].values, _da(expected, time, layer).values) + # the zero-chloride layer is the pure identity: any scaling of the density term breaks it + np.testing.assert_array_equal(ds["freshwater_head"].isel(layer=0).values, ds["head_filled"].isel(layer=0).values) + + +def test_add_output_concentration_mean_and_grensvlak(monkeypatch): + """The mean is thickness-weighted and the two thresholds reach the right grensvlak. + + Chloride 500/1500/14500 crosses 1000 exactly halfway between the layer-0 and layer-1 + centres and 8000 exactly halfway between the layer-1 and layer-2 centres. + """ + ds, time = _output_ds() + layer = ds.layer.values + conc = _da(_PROFILE, time, layer) + conc[1] *= 2.0 # a second time step so the change vs. t=0 is not trivially zero + _patch_loaders(monkeypatch, _da([1.0, 1.0, 1.0], time, layer), conc) + + ds, ctop = pp.add_output_to_ds(ds, "ws", "test") + + # thickness-weighted: 7750; an unweighted mean of the same profile would give 5500 instead + mean0 = float(_PROFILE @ _THICKNESS) / _THICKNESS.sum() + assert mean0 == _MEAN_T0 + expected_mean = np.full((2, 2), mean0) + expected_mean[1] = 2.0 * mean0 + np.testing.assert_array_equal(ds["concentration_mean"].values, expected_mean) + np.testing.assert_array_equal(ds["dconcentration_mean"].values, expected_mean - mean0) + + # midpoint crossings: fresh between centres -2 and -6, brackish between -6 and -12 + zoet = (_Z_CENTRE[0] + _Z_CENTRE[1]) / 2.0 + brak = (_Z_CENTRE[1] + _Z_CENTRE[2]) / 2.0 + assert (zoet, brak) == (-4.0, -9.0) + np.testing.assert_array_equal(ds["grensvlak_zoet"].isel(time=0).values, np.full(2, zoet)) + np.testing.assert_array_equal(ds["grensvlak_brak"].isel(time=0).values, np.full(2, brak)) + # the fresh interface is never deeper than the brackish one, and the thresholds are not swapped + assert (ds["grensvlak_zoet"] >= ds["grensvlak_brak"]).all() + assert ds["grensvlak_zoet"].attrs["threshold"] == _THRESHOLD + assert ds["grensvlak_brak"].attrs["threshold"] == _THRESHOLD_BRAK + + np.testing.assert_array_equal(ctop.values, conc.isel(layer=0).values) + + +def test_add_output_without_transport_skips_the_ucn(monkeypatch): + """``transport=False`` returns ``ctop=None`` and never reads the transport output.""" + ds, time = _output_ds() + layer = ds.layer.values + calls = _patch_loaders(monkeypatch, _da([1.0, 1.0, 1.0], time, layer), _da([0.0, 0.0, 0.0], time, layer)) + + ds, ctop = pp.add_output_to_ds(ds, "ws", "test", transport=False) + + assert ctop is None + assert calls == [] + assert "concentration" not in ds + assert "freshwater_head" not in ds + + +# --- plot_result_maps ---------------------------------------------------------------------- + +_GRID_MAP = {"doorsnedelijnen.png"} +_DRN_MAP = {"oppervlaktewater.png"} +# nper=3 with iper=-1 must normalise to 2; a sign flip instead of nper+iper would give 1 +_TRANSPORT_MAPS = { + "map_head_L0_t2.png", + "map_conc_L0_t2.png", + "grensvlak_zoet_t2.png", + "grensvlak_brak_t2.png", +} + + +def _plot_ds(*, transport, drn): + """Build a 2x2 vertex ds with three stress periods and the fields plot_result_maps consumes.""" + ds = make_rect_vertex_ds() + time = pd.to_datetime(["2022-01-01", "2022-02-01", "2022-03-01"]) + ds = ds.assign_coords(time=time) + ncell = ds.sizes["icell2d"] + ctop = None + if drn: + ds["drn_elev"] = xr.DataArray(np.full(ncell, -1.0), dims=("icell2d",)) + if transport: + ds["freshwater_head"] = xr.DataArray( + np.zeros((3, ds.sizes["layer"], ncell)), dims=("time", "layer", "icell2d"), coords={"time": time} + ) + for name, threshold in (("zoet", _THRESHOLD), ("brak", _THRESHOLD_BRAK)): + da = xr.DataArray(np.full((3, ncell), -10.0), dims=("time", "icell2d"), coords={"time": time}) + da.attrs["threshold"] = threshold + ds[f"grensvlak_{name}"] = da + ctop = xr.DataArray(np.full((3, ncell), 100.0), dims=("time", "icell2d"), coords={"time": time}) + return ds, ctop + + +@pytest.mark.parametrize( + ("transport", "drn", "expected"), + [ + (True, False, _GRID_MAP | _TRANSPORT_MAPS), + (True, True, _GRID_MAP | _TRANSPORT_MAPS | _DRN_MAP), + (False, False, _GRID_MAP), # early return: no freshwater_head, ctop=None + ], +) +def test_plot_result_maps_filenames(tmp_path, transport, drn, expected): + """Exactly the expected figures are written, with iper=-1 normalised against nper=3.""" + ds, ctop = _plot_ds(transport=transport, drn=drn) + # named seam: add_background_map fetches contextily tiles, the module's only network call + with mock.patch.object(pp.nlmod.plot, "add_background_map"): + pp.plot_result_maps(ds, ctop, str(tmp_path)) + + assert {p.name for p in tmp_path.iterdir()} == expected diff --git a/tests/test_pwnlayers3_layers.py b/tests/test_pwnlayers3_layers.py new file mode 100644 index 0000000..4a499db --- /dev/null +++ b/tests/test_pwnlayers3_layers.py @@ -0,0 +1,586 @@ +"""Tests for the PWN layer-model seams of ``nhflotools.pwnlayers3.layers``. + +The pure helpers are pinned with hand-derived arithmetic: the botm repair (which lives as +two divergent copies), thickness telescoping, the zero-thickness guard, the in-place +griddata fill, and the REGIS/PWN merge that splits, routes and interpolates layers. +""" + +import geopandas as gpd +import nlmod +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from shapely.geometry import Point, box + +from nhflotools.pwnlayers.merge_layer_models import ( + _apply_ratios_to_botm, + _compute_thickness_ratios, + _interpolate_da, + combine_two_layer_models, +) +from nhflotools.pwnlayers.utils import ( + fix_missings_botms_and_min_layer_thickness as fix_botms_utils, +) +from nhflotools.pwnlayers3.layers import _guard_zero_thickness, get_kv, get_pwn_layer_model, get_thickness +from nhflotools.pwnlayers3.layers import ( + fix_missings_botms_and_min_layer_thickness as fix_botms_pwnlayers3, +) +from nhflotools.pwnlayers3.layers import layer_names as PWN_LAYER_NAMES + +from .util import make_rect_vertex_ds + +# Header names combine_two_layer_models defaults to; the koppeltabel columns must match. +_H_REGIS = "Regis II v2.2" +_H_OTHER = "ASSUMPTION1" +_CRS = "EPSG:28992" + + +def _cell_coords(ds): + """Coordinate mapping shared by every DataArray built on the grid of ``ds``.""" + return {"icell2d": ds.icell2d, "x": ds.x, "y": ds.y} + + +def _spread(values, ncell): + """Turn a per-layer sequence of scalars into a ``(layer, icell2d)`` array.""" + return np.asarray(values, dtype=float)[:, None] * np.ones((1, ncell)) + + +def _layer_model(ds, layers, botm, kh, kv): + """Build a REGIS-shaped layer model with uniform values per layer.""" + ncell = ds.sizes["icell2d"] + return xr.Dataset( + { + "botm": (("layer", "icell2d"), _spread(botm, ncell)), + "kh": (("layer", "icell2d"), _spread(kh, ncell)), + "kv": (("layer", "icell2d"), _spread(kv, ncell)), + }, + coords={"layer": list(layers), **_cell_coords(ds)}, + attrs={"extent": ds.attrs["extent"], "gridtype": ds.attrs["gridtype"]}, + ) + + +def _flags(ds, layers, cells_true): + """Boolean mask/transition Dataset that is True on ``cells_true`` for every layer.""" + data = np.zeros((len(layers), ds.sizes["icell2d"]), dtype=bool) + data[:, list(cells_true)] = True + return xr.Dataset( + {var: (("layer", "icell2d"), data.copy()) for var in ("botm", "kh", "kv")}, + coords={"layer": list(layers), **_cell_coords(ds)}, + ) + + +def _top(ds, value=0.0): + return xr.DataArray(np.full(ds.sizes["icell2d"], float(value)), dims="icell2d", coords=_cell_coords(ds), name="top") + + +def _koppeltabel(rows): + return pd.DataFrame(list(rows), columns=[_H_REGIS, _H_OTHER]) + + +@pytest.mark.parametrize( + "fix_botms", + [fix_botms_pwnlayers3, fix_botms_utils], + ids=["pwnlayers3.layers", "pwnlayers.utils"], +) +def test_fix_missings_botms_is_pure_monotone_and_idempotent(fix_botms): + """Both divergent copies fill NaNs downward, clip crossings and leave the input alone. + + Parametrized over the two copies on purpose: editing one and not the other is the + standing maintenance trap in this codebase. + """ + top = xr.DataArray(np.zeros(3), dims="icell2d", coords={"icell2d": [0, 1, 2]}) + botm = xr.DataArray( + np.array([ + [-10.0, -10.0, np.nan], + [np.nan, -5.0, np.nan], + [-30.0, -20.0, -8.0], + ]), + dims=("layer", "icell2d"), + coords={"layer": ["a", "b", "c"], "icell2d": [0, 1, 2]}, + ) + # ffill down the column with `top` prepended, then a running minimum: + # cell 0: [0, -10, nan, -30] -> ffill [0, -10, -10, -30] -> already decreasing. + # cell 1: [0, -10, -5, -20] -> -5 lies above -10, so it is pulled down to -10. + # cell 2: [0, nan, nan, -8] -> leading NaNs are filled from the top elevation, 0. + expected = np.array([ + [-10.0, -10.0, 0.0], + [-10.0, -10.0, 0.0], + [-30.0, -20.0, -8.0], + ]) + before = botm.copy(deep=True) + + out = fix_botms(top=top, botm=botm) + + np.testing.assert_array_equal(out.values, expected) + assert out.dims == ("layer", "icell2d") + assert not out.isnull().any() + assert (out.diff(dim="layer") <= 0.0).all(), "layer bottoms must be non-increasing downward" + assert (out <= top).all() + # Purity: the caller keeps only the return value, so mutating the argument would corrupt + # the source model silently. + np.testing.assert_array_equal(botm.values, before.values) + # Idempotence: re-running the repair on repaired botms must be a no-op. + np.testing.assert_array_equal(fix_botms(top=top, botm=out).values, expected) + + with pytest.raises(ValueError, match="nan"): + fix_botms(top=top.where(top.icell2d != 1), botm=botm) + + +def test_get_thickness_telescopes_and_labels_the_lower_layer(): + """thickness[k] == botm[k-1] - botm[k], labelled with the lower layer; the top layer drops out.""" + botm = xr.DataArray( + np.array([[-2.0, -2.0], [-6.0, -10.0], [-14.0, -12.0]]), + dims=("layer", "icell2d"), + coords={"layer": ["W11", "S11", "W12"], "icell2d": [0, 1]}, + ) + thickness = get_thickness(botm=botm) + + # W11 needs the model top and is therefore absent; the label is the *lower* of each pair. + assert list(thickness.layer.values) == ["S11", "W12"] + np.testing.assert_array_equal(thickness.values, np.array([[4.0, 8.0], [8.0, 2.0]])) + # Telescoping: the column sums back to the distance between the first and last bottom. + np.testing.assert_array_equal( + thickness.sum(dim="layer").values, botm.isel(layer=0).values - botm.isel(layer=-1).values + ) + + +def test_guard_zero_thickness_replaces_only_the_isclose_zero_cells(): + """Cells within np.isclose of zero thickness take the fill value; the rest are untouched.""" + values = np.array([1.0, 2.0, 3.0, 4.0]) + # 1e-12 is inside np.isclose's default atol of 1e-8, 0.5 and 2.0 are not. + guarded = _guard_zero_thickness(values.copy(), np.array([0.0, 1e-12, 0.5, 2.0]), 7.0, "S11") + np.testing.assert_array_equal(guarded, [7.0, 7.0, 3.0, 4.0]) + + identity = _guard_zero_thickness(values.copy(), np.array([0.5, 1.0, 2.0, 4.0]), 7.0, "S11") + np.testing.assert_array_equal(identity, values) + + +def test_interpolate_da_writes_through_to_the_parent_dataset(): + """The fill reaches the parent Dataset, honours the method, and no-ops when nothing is missing. + + ``_interpolate_ds`` hands ``_interpolate_da`` a ``.sel(layer=...)`` view and relies on the + ``.loc`` assignment propagating back; an xarray copy-semantics change would silently turn + the whole transition interpolation into a no-op. + """ + # Four corners of a 4x4 square carry data; the cell at (1, 1) is missing. + x = np.array([0.0, 4.0, 0.0, 4.0, 1.0]) + y = np.array([0.0, 0.0, 4.0, 4.0, 1.0]) + plane = x + 2.0 * y # any linear interpolant reproduces an affine field exactly + sentinel = -999.0 + values = np.stack([plane, plane]) + values[:, 4] = sentinel + ds = xr.Dataset( + {"kh": (("layer", "icell2d"), values)}, + coords={"layer": ["A", "B"], "icell2d": np.arange(5), "x": ("icell2d", x), "y": ("icell2d", y)}, + ) + isvalid = xr.DataArray([True] * 4 + [False], dims="icell2d", coords={"icell2d": np.arange(5)}) + ismissing = ~isvalid + nothing_missing = xr.zeros_like(isvalid, dtype=bool) + + _interpolate_da(ds["kh"].sel(layer="A"), isvalid=isvalid, ismissing=ismissing, method="linear") + _interpolate_da(ds["kh"].sel(layer="B"), isvalid=isvalid, ismissing=ismissing, method="nearest") + + # rtol accommodates Qhull barycentric arithmetic; the field is affine so the plane value + # 1 + 2*1 = 3 is triangulation-independent. + np.testing.assert_allclose(ds["kh"].sel(layer="A").values[4], 3.0, rtol=1e-12) + # Nearest donor of (1, 1) is the corner (0, 0), at distance sqrt(2) against 3.16 for the others. + assert ds["kh"].sel(layer="B").values[4] == 0.0 + # Valid cells are never rewritten. + np.testing.assert_array_equal(ds["kh"].values[:, :4], np.stack([plane[:4], plane[:4]])) + + ds["kh"].loc[{"layer": "A", "icell2d": 4}] = sentinel + _interpolate_da(ds["kh"].sel(layer="A"), isvalid=isvalid, ismissing=nothing_missing, method="linear") + assert ds["kh"].sel(layer="A").values[4] == sentinel + + +def test_combine_two_layer_models_reduces_to_regis_when_pwn_is_absent(): + """An all-False mask with a 1:1 koppeltabel must hand back REGIS untouched.""" + ds = make_rect_vertex_ds(nx=2, ny=2) + top = _top(ds) + regis = _layer_model(ds, ["A", "B", "C"], botm=[-8.0, -16.0, -24.0], kh=[1.0, 2.0, 4.0], kv=[0.1, 0.2, 0.4]) + other = _layer_model(ds, ["p", "q", "r"], botm=[-2.0, -6.0, -12.0], kh=[16.0, 32.0, 64.0], kv=[1.6, 3.2, 6.4]) + + out, cat = combine_two_layer_models( + layer_model_regis=regis, + layer_model_other=other, + mask_model_other=_flags(ds, ["p", "q", "r"], cells_true=[]), + transition_model=_flags(ds, ["p", "q", "r"], cells_true=[]), + top=top, + df_koppeltabel=_koppeltabel([("A", "p"), ("B", "q"), ("C", "r")]), + split_method="nearest_ratio", + ) + + for var in ("botm", "kh", "kv"): + np.testing.assert_array_equal(out[var].values, regis[var].values) + for var in ("botm", "kh", "kv"): + np.testing.assert_array_equal(cat[var].values, np.ones_like(regis[var].values, dtype=int)) + + +@pytest.mark.parametrize( + ("split_method", "botm_a1_regis_cells"), + [ + # 'equal' halves the 8 m REGIS layer A; 'nearest_ratio' takes the 2 m : 6 m split of the + # PWN pair (p, q) and applies it to A, giving 0 - (2/8)*8 = -2. + ("equal", -4.0), + ("nearest_ratio", -2.0), + ], +) +def test_combine_two_layer_models_routes_and_conserves_split_thickness(split_method, botm_a1_regis_cells): + """Split layers take their values from the right source model and preserve group bottoms. + + The koppeltabel mixes a 1:2 REGIS split (A -> p, q), a 2:1 PWN split (B, C -> r) and one + uncoupled deep REGIS layer (D). Cells 0 and 1 are PWN, cells 2 and 3 are REGIS. + """ + ds = make_rect_vertex_ds(nx=2, ny=2) + top = _top(ds) + regis = _layer_model( + ds, ["A", "B", "C", "D"], botm=[-8.0, -16.0, -24.0, -40.0], kh=[1.0, 2.0, 4.0, 8.0], kv=[0.1, 0.2, 0.4, 0.8] + ) + other = _layer_model(ds, ["p", "q", "r"], botm=[-2.0, -8.0, -20.0], kh=[16.0, 32.0, 64.0], kv=[1.6, 3.2, 6.4]) + pwn_cells, regis_cells = [0, 1], [2, 3] + + out, cat = combine_two_layer_models( + layer_model_regis=regis, + layer_model_other=other, + mask_model_other=_flags(ds, ["p", "q", "r"], cells_true=pwn_cells), + transition_model=_flags(ds, ["p", "q", "r"], cells_true=[]), + top=top, + df_koppeltabel=_koppeltabel([("A", "p"), ("A", "q"), ("B", "r"), ("C", "r"), ("D", np.nan)]), + split_method=split_method, + ) + + assert list(out.layer.values) == ["A_1", "A_2", "B_1", "C_1", "D"] + # Category 2 = PWN, 1 = REGIS; the uncoupled layer D is REGIS everywhere. + np.testing.assert_array_equal(cat["kh"].values, np.array([[2, 2, 1, 1]] * 4 + [[1, 1, 1, 1]]).reshape(5, 4)) + + # Routing: sublayers inherit the conductivity of the layer they were split from. + kh_pwn = np.array([16.0, 32.0, 64.0, 64.0, 8.0]) # p, q, r, r, REGIS D + kh_regis = np.array([1.0, 1.0, 2.0, 4.0, 8.0]) # A, A, B, C, D + np.testing.assert_array_equal(out["kh"].sel(icell2d=pwn_cells).values, np.tile(kh_pwn[:, None], (1, 2))) + np.testing.assert_array_equal(out["kh"].sel(icell2d=regis_cells).values, np.tile(kh_regis[:, None], (1, 2))) + np.testing.assert_array_equal(out["kv"].sel(icell2d=pwn_cells).values, np.tile(kh_pwn[:, None] / 10.0, (1, 2))) + + # Bottoms: PWN cells get p, q and the equal halves of r between -8 and -20; REGIS cells get + # the A split (method dependent) plus the untouched B, C. D is REGIS in every cell. + botm_pwn = np.array([-2.0, -8.0, -14.0, -20.0, -40.0]) + botm_regis = np.array([botm_a1_regis_cells, -8.0, -16.0, -24.0, -40.0]) + np.testing.assert_array_equal(out["botm"].sel(icell2d=pwn_cells).values, np.tile(botm_pwn[:, None], (1, 2))) + np.testing.assert_array_equal(out["botm"].sel(icell2d=regis_cells).values, np.tile(botm_regis[:, None], (1, 2))) + + # Conservation: splitting redistributes thickness inside a group but never changes the + # group's total thickness, i.e. the group bottom equals the source layer's bottom. + thickness = get_thickness(botm=xr.concat([top.expand_dims(layer=["mv"]), out["botm"]], dim="layer")) + group_a = thickness.sel(layer=["A_1", "A_2"]).sum(dim="layer") + group_r = thickness.sel(layer=["B_1", "C_1"]).sum(dim="layer") + np.testing.assert_array_equal(group_a.sel(icell2d=regis_cells).values, [8.0, 8.0]) # REGIS A: 0 - -8 + np.testing.assert_array_equal(group_a.sel(icell2d=pwn_cells).values, [8.0, 8.0]) # PWN p+q: 0 - -8 + np.testing.assert_array_equal(group_r.sel(icell2d=pwn_cells).values, [12.0, 12.0]) # PWN r: -8 - -20 + np.testing.assert_array_equal(group_r.sel(icell2d=regis_cells).values, [16.0, 16.0]) # REGIS B+C + + +def test_combine_two_layer_models_interpolates_the_transition_band(): + """Transition cells are interpolated between the PWN and REGIS values, not copied from either. + + The 3x3 grid is banded north to south: PWN row, transition row, REGIS row. Because each band + carries one constant value, the valid data is an affine function of y and the interpolated + middle row is the exact midpoint whatever triangulation scipy picks. + """ + ds = make_rect_vertex_ds(nx=3, ny=3) + top = _top(ds) + regis = _layer_model(ds, ["A"], botm=[-16.0], kh=[2.0], kv=[0.5]) + other = _layer_model(ds, ["p"], botm=[-4.0], kh=[8.0], kv=[2.0]) + pwn_cells, transition_cells, regis_cells = [0, 1, 2], [3, 4, 5], [6, 7, 8] + + out, cat = combine_two_layer_models( + layer_model_regis=regis, + layer_model_other=other, + mask_model_other=_flags(ds, ["p"], cells_true=pwn_cells), + transition_model=_flags(ds, ["p"], cells_true=transition_cells), + top=top, + df_koppeltabel=_koppeltabel([("A", "p")]), + split_method="nearest_ratio", + ) + + np.testing.assert_array_equal(cat["kh"].values, [[2, 2, 2, 3, 3, 3, 1, 1, 1]]) + # rtol accommodates Qhull barycentric arithmetic on the affine field. + for var, pwn_value, regis_value in [("kh", 8.0, 2.0), ("kv", 2.0, 0.5), ("botm", -4.0, -16.0)]: + values = out[var].sel(layer="A_1").values + np.testing.assert_array_equal(values[pwn_cells], [pwn_value] * 3) + np.testing.assert_array_equal(values[regis_cells], [regis_value] * 3) + np.testing.assert_allclose(values[transition_cells], [(pwn_value + regis_value) / 2.0] * 3, rtol=1e-12) + + +def test_thickness_ratios_round_trip_through_apply_ratios_to_botm(): + """Ratios reproduce the source split exactly, fall back to 1/N and spread by nearest neighbour. + + The group starts at the second layer, so ``_compute_thickness_ratios`` and + ``_apply_ratios_to_botm`` must both take the group top from the layer above (the + ``first_idx - 1`` lookup) rather than from the model top. + """ + ds = make_rect_vertex_ds(nx=3, ny=1) + top = _top(ds) + layers = ["X", "Y", "Z"] + botm = np.array([ + [-4.0, -4.0, -4.0], # X + [-8.0, -4.0, -8.0], # Y: cell 1 is a collapsed, zero-thickness group + [-20.0, -4.0, -20.0], # Z + ]) + source = xr.Dataset( + {"botm": (("layer", "icell2d"), botm), "top": top}, + coords={"layer": layers, **_cell_coords(ds)}, + ) + # Cell 2 has no source data and must inherit from its nearest valid neighbour, cell 1 + # (100 m away) rather than cell 0 (200 m away). + mask_valid = xr.DataArray([True, True, False], dims="icell2d", coords=_cell_coords(ds)) + + ratios = _compute_thickness_ratios(source, ["Y", "Z"], mask_valid) + + # cell 0: group top is botm X = -4, so Y is 4 m of the 16 m group and Z is 12 m -> 1/4, 3/4. + # cell 1: zero group thickness -> equal ratios 1/2. cell 2: nearest copy of cell 1. + np.testing.assert_array_equal( + ratios.transpose("layer", "icell2d").values, np.array([[0.25, 0.5, 0.5], [0.75, 0.5, 0.5]]) + ) + np.testing.assert_array_equal(ratios.sum(dim="layer").values, np.ones(3)) + + target = source.copy(deep=True) + _apply_ratios_to_botm(target, top, ["Y", "Z"], ratios) + + # cell 0 round-trips to its own botm; cell 1 stays collapsed; cell 2 gets the halved split + # of its 16 m group: -4 - 0.5*16 = -12. + np.testing.assert_array_equal(target["botm"].sel(layer="Y").values, [-8.0, -4.0, -12.0]) + # The group bottom is preserved in every cell, so total group thickness is unchanged. + np.testing.assert_array_equal(target["botm"].sel(layer="Z").values, botm[2]) + + +def _write_resistance_polygons(data_path, values, boxes): + """Write one C_combined.geojson per aquitard, holding constant-c polygons.""" + gdf = gpd.GeoDataFrame({"VALUE": list(values)}, geometry=[box(*b) for b in boxes], crs=_CRS) + conductances = data_path / "conductances" + conductances.mkdir(parents=True, exist_ok=True) + for name in PWN_LAYER_NAMES: + if name.startswith("S"): + gdf.to_file(conductances / f"C{name}_combined.geojson", driver="GeoJSON") + + +def test_get_kv_uses_harmonic_area_weighting_and_anisotropy(tmp_path): + """Aquitards get kv = d/c with 1/c averaged by area; aquifers get kv = kh/anisotropy.""" + ds = make_rect_vertex_ds(nx=2, ny=1) + # Cell 0 spans x in [0, 100] and is halved by the polygon boundary at x = 50; cell 1 + # (x in [100, 200]) lies wholly inside the second polygon. + _write_resistance_polygons(tmp_path, values=[2.0, 4.0], boxes=[(0, 0, 50, 200), (50, 0, 200, 200)]) + + nlay, ncell = len(PWN_LAYER_NAMES), ds.sizes["icell2d"] + thickness = xr.DataArray( + np.full((nlay, ncell), 8.0), + dims=("layer", "icell2d"), + coords={"layer": PWN_LAYER_NAMES, "icell2d": ds.icell2d}, + ) + thickness.loc[{"layer": "S11", "icell2d": 1}] = 0.0 # collapsed cell -> fill value + kh = xr.DataArray( + np.full((nlay, ncell), 30.0), + dims=("layer", "icell2d"), + coords={"layer": PWN_LAYER_NAMES, "icell2d": ds.icell2d}, + ) + + kv = get_kv( + ds=ds, + data_path_2024=tmp_path, + kh=kh, + thickness=thickness, + anisotropy=10.0, + fill_value_kv=1.0, + isin_bounds=np.ones((nlay, ncell), dtype=bool), + ) + + # Aquifers: 30 / 10. + for name in [n for n in PWN_LAYER_NAMES if n.startswith("W")]: + np.testing.assert_array_equal(kv.sel(layer=name).values, [3.0, 3.0]) + # Aquitards: 1/c averaged over the cell, then kv = d * (1/c). + # cell 0: 0.5*(1/2) + 0.5*(1/4) = 0.375 -> 8 * 0.375 = 3.0. cell 1: 1/4 -> 8 * 0.25 = 2.0. + for name in [n for n in PWN_LAYER_NAMES if n.startswith("S") and n != "S11"]: + np.testing.assert_array_equal(kv.sel(layer=name).values, [3.0, 2.0]) + np.testing.assert_array_equal(kv.sel(layer="S11").values, [3.0, 1.0]) + + +# ── Offline integration of get_pwn_layer_model ────────────────────────────────────────── +# A synthetic 4x4 grid of 100 m cells over [0, 400]^2, with a hand-written PWN data tree. +# The PWN boundary covers x < 200 (columns 0 and 1), so with a 150 m transition buffer the +# grid splits into PWN columns 0-1, a transition column 2 and a REGIS-only column 3. The +# NHDZ region covers y > 200, splitting the PWN block into a KD half and a Bergen half. +_MASKED_NHDZ, _MASKED_BERGEN = [0, 1, 4, 5], [8, 9, 12, 13] +_TRANSITION_CELLS, _REGIS_CELLS = [2, 6, 10, 14], [3, 7, 11, 15] +_REGIS_LAYERS = ["mv", "A", "B", "C", "D", "E"] +_KOPPELTABEL_GROUPS = [("A", 4), ("B", 4), ("C", 4), ("D", 2)] +# Layer k of the PWN model has its bottom at -2*(k+1) minus a x/200 tilt, so every PWN layer +# is exactly 2 m thick and linear interpolation of the source points is exact. +_PWN_BASE_BOTM = -2.0 * (np.arange(len(PWN_LAYER_NAMES)) + 1.0) +_C_VALUE, _KD_VALUE, _K_VALUE, _ANISOTROPY, _REGIS_KH = 4.0, 8.0, 20.0, 10.0, 5.0 + + +def _write_pwn_data_tree(path): + """Write the boundary, botm, conductance and NHDZ GeoJSONs get_ds reads.""" + points = [Point(x, y) for x in (0.0, 100.0, 200.0, 300.0, 400.0) for y in (0.0, 100.0, 200.0, 300.0, 400.0)] + xs = np.array([p.x for p in points]) + + (path / "botm").mkdir(parents=True, exist_ok=True) + botm_columns = {name: _PWN_BASE_BOTM[i] - xs / 200.0 for i, name in enumerate(PWN_LAYER_NAMES)} + gpd.GeoDataFrame(botm_columns, geometry=points, crs=_CRS).to_file(path / "botm" / "botm.geojson", driver="GeoJSON") + + for unit in ["11", "12", "13", "21", "22", "31", "32"]: + (path / "boundaries" / f"S{unit}").mkdir(parents=True, exist_ok=True) + gpd.GeoDataFrame(geometry=[box(-1.0, -1.0, 201.0, 401.0)], crs=_CRS).to_file( + path / "boundaries" / f"S{unit}" / f"S{unit}.geojson", driver="GeoJSON" + ) + gpd.GeoDataFrame(geometry=[box(-1.0, 199.0, 401.0, 401.0)], crs=_CRS).to_file( + path / "boundaries" / "triwaco_model_nhdz.geojson", driver="GeoJSON" + ) + + _write_resistance_polygons(path, values=[_C_VALUE], boxes=[(-1.0, -1.0, 401.0, 401.0)]) + conductances = path / "conductances" + for name in PWN_LAYER_NAMES: + if name.startswith("W"): + gpd.GeoDataFrame({"VALUE": [_K_VALUE]}, geometry=[box(-1.0, -1.0, 401.0, 401.0)], crs=_CRS).to_file( + conductances / f"K{name}_combined.geojson", driver="GeoJSON" + ) + for name in ["S12", "S13", "S21", "S22", "S31"]: + gpd.GeoDataFrame({"VALUE": [_KD_VALUE] * len(points)}, geometry=points, crs=_CRS).to_file( + conductances / f"KD{name}_NHDZ.geojson", driver="GeoJSON" + ) + + +def _regis_input(): + """4x4 vertex ds standing in for REGIS, plus the model top and the koppeltabel rows.""" + ds = make_rect_vertex_ds(nx=4, ny=4, botm=(0.0, -10.0, -20.0, -30.0, -40.0, -120.0), kh=_REGIS_KH) + ds = ds.assign_coords(layer=_REGIS_LAYERS) + regis_per_pwn_layer = [group for group, n in _KOPPELTABEL_GROUPS for _ in range(n)] + rows = list(zip(regis_per_pwn_layer, PWN_LAYER_NAMES, strict=True)) + return ds, _top(ds), _koppeltabel([*rows, ("E", np.nan)]) + + +def _expected_pwn_kh(in_nhdz): + """Kh the PWN model should produce per layer: K polygon, KD/d, or d*anisotropy/c.""" + d, inv_c = 2.0, 1.0 / _C_VALUE + kh_from_c = d * _ANISOTROPY * inv_c + return np.array([ + _K_VALUE + if name.startswith("W") + else (kh_from_c if name in {"S11", "S32"} else (_KD_VALUE / d if in_nhdz else kh_from_c)) + for name in PWN_LAYER_NAMES + ]) + + +def test_get_pwn_layer_model_offline_known_answers(tmp_path, monkeypatch): + """The merged model carries the PWN values, in koppeltabel order, everywhere PWN is valid. + + ``nlmod.read.regis.get_layer_names`` is patched because it opens an OPeNDAP endpoint when + called without a dataset; it is the only network seam in this pipeline. + """ + monkeypatch.setattr(nlmod.read.regis, "get_layer_names", lambda: pd.Index(_REGIS_LAYERS, name="layer")) + _write_pwn_data_tree(tmp_path) + ds_regis, top, koppeltabel = _regis_input() + fname_koppeltabel = tmp_path / "koppeltabel.csv" + koppeltabel.to_csv(fname_koppeltabel) + + out = get_pwn_layer_model( + ds_regis=ds_regis, + data_path_2024=tmp_path, + fname_koppeltabel=fname_koppeltabel, + top=top, + anisotropy=_ANISOTROPY, + distance_transition=150.0, + return_diagnostics=True, + ) + + # Split names follow the koppeltabel: A splits into 4, ..., D into 2, E stays uncoupled. + expected_layers = [f"{g}_{i + 1}" for g, n in _KOPPELTABEL_GROUPS for i in range(n)] + ["E"] + assert list(out.layer.values) == expected_layers + + coupled = out.sel(layer=expected_layers[:-1]) + np.testing.assert_array_equal(coupled["cat_kh"].values[:, _MASKED_NHDZ + _MASKED_BERGEN], 2) + np.testing.assert_array_equal(coupled["cat_kh"].values[:, _TRANSITION_CELLS], 3) + np.testing.assert_array_equal(coupled["cat_kh"].values[:, _REGIS_CELLS], 1) + np.testing.assert_array_equal(out["cat_kh"].sel(layer="E").values, 1) + + # rtol accommodates Qhull barycentric arithmetic in the botm/KD interpolations, which + # propagates into the thickness the conductivities are derived from. + for cells, in_nhdz in [(_MASKED_NHDZ, True), (_MASKED_BERGEN, False)]: + kh_expected = _expected_pwn_kh(in_nhdz) + np.testing.assert_allclose( + coupled["kh"].values[:, cells], np.tile(kh_expected[:, None], (1, len(cells))), rtol=1e-12 + ) + # Aquifers keep kv = kh/anisotropy; aquitards get kv = d/c = 2/4. + kv_expected = np.where([n.startswith("W") for n in PWN_LAYER_NAMES], kh_expected / _ANISOTROPY, 2.0 / _C_VALUE) + np.testing.assert_allclose( + coupled["kv"].values[:, cells], np.tile(kv_expected[:, None], (1, len(cells))), rtol=1e-12 + ) + # Bottoms come from the PWN point cloud: -2*(k+1) tilted by -x/200. + x = out.x.values[cells] + np.testing.assert_allclose( + coupled["botm"].values[:, cells], _PWN_BASE_BOTM[:, None] - x[None, :] / 200.0, rtol=1e-12 + ) + + # Transition cells mix both models: strictly between REGIS kh (5) and PWN aquifer kh (20). + aquifer_layers = [ + layer for layer, name in zip(expected_layers[:-1], PWN_LAYER_NAMES, strict=True) if name.startswith("W") + ] + kh_transition = out["kh"].sel(layer=aquifer_layers).values[:, _TRANSITION_CELLS] + assert np.all(kh_transition > _REGIS_KH) + assert np.all(kh_transition < _K_VALUE) + + # Postconditions the NPF package depends on: strictly positive, finite, monotone. + assert np.all(out["kh"].values > 0.0) + assert np.all(out["kv"].values > 0.0) + assert np.all(np.isfinite(out["kh"].values)) + assert np.all(np.isfinite(out["kv"].values)) + assert (out["botm"].diff(dim="layer") <= 0.0).all() + assert (out["botm"] <= out["top"]).all() + np.testing.assert_array_equal(out["top"].values, top.values) + # Values passthrough only: layers.py evaluates the get_area default eagerly (issue #61). + np.testing.assert_array_equal(out["area"].values, ds_regis["area"].values) + + # The plot module ticks its colourbars off these attributes, and every code that actually + # occurs must be declared. Seeing 1/2/4/5 also proves all three kh branches ran. + assert set(np.unique(out["kh_method"].values)) == {0, 1, 2, 4, 5} + for name in ["botm_method", "kh_method", "kv_method"]: + flag_values = set(out[name].attrs["flag_values"]) + assert set(np.unique(out[name].values)) <= flag_values + assert len(out[name].attrs["flag_meanings"].split(";")) == len(flag_values) + + +@pytest.mark.parametrize( + ("bad", "match"), + [("top", "should not contain nan"), ("layers", "All REGIS layers should be present")], +) +def test_get_pwn_layer_model_rejects_nan_top_and_wrong_regis_layers(bad, match, monkeypatch): + """The two fail-fast guards trip before any data is read.""" + monkeypatch.setattr(nlmod.read.regis, "get_layer_names", lambda: pd.Index(_REGIS_LAYERS, name="layer")) + ds_regis, top, _ = _regis_input() + if bad == "top": + top = top.where(top.icell2d != 0) + else: + # Same number of layers, one renamed: the merge would otherwise silently couple the + # koppeltabel to the wrong REGIS unit. + ds_regis = ds_regis.assign_coords(layer=[*_REGIS_LAYERS[:-1], "X"]) + + with pytest.raises(ValueError, match=match): + get_pwn_layer_model(ds_regis=ds_regis, data_path_2024=None, fname_koppeltabel=None, top=top) + + +@pytest.mark.xfail( + strict=True, + reason="The guard compares layer names elementwise, so a missing layer raises pandas' " + "'Lengths must match to compare' before the actionable message is reached. NHFLO/tools#65", +) +def test_get_pwn_layer_model_reports_a_missing_regis_layer(monkeypatch): + """A dropped REGIS layer must be reported with the message that names the fix. + + This is what ``get_regis(.., remove_nan_layers=True)`` hands over, so it is the case + the guard exists for -- and the one it currently fails to report usefully. + """ + monkeypatch.setattr(nlmod.read.regis, "get_layer_names", lambda: pd.Index(_REGIS_LAYERS, name="layer")) + ds_regis, top, _ = _regis_input() + ds_regis = ds_regis.sel(layer=_REGIS_LAYERS[:-1]) + + with pytest.raises(ValueError, match="All REGIS layers should be present"): + get_pwn_layer_model(ds_regis=ds_regis, data_path_2024=None, fname_koppeltabel=None, top=top) diff --git a/tests/test_pwnlayers3_plot.py b/tests/test_pwnlayers3_plot.py new file mode 100644 index 0000000..315d2ce --- /dev/null +++ b/tests/test_pwnlayers3_plot.py @@ -0,0 +1,182 @@ +"""Tests for the PWN layer model cross-section plotting helpers.""" + +import geopandas +import numpy as np +import pytest +from matplotlib.figure import Figure +from shapely.geometry import Point + +from nhflotools.pwnlayers3.layers import layer_names +from nhflotools.pwnlayers3.plot import ( + _load_and_project_source_botm, # noqa: PLC2701 + _overlay_source_botm, # noqa: PLC2701 + _parse_flag_labels, # noqa: PLC2701 +) + +# The three flag_meanings strings that pwnlayers3.layers attaches to the diagnostic +# method arrays, copied verbatim from layers.py (botm_method, kh_method, kv_method). +# They are duplicated here on purpose: the test is only meaningful if its input is an +# independent transcription of what production attaches, not an import of it. +BOTM_METHOD_MEANINGS = ( + "0: no_data (outside boundary polygon); " + "1: linear_interpolation (griddata linear from botm.geojson point data); " + "2: nearest_interpolation (griddata nearest fallback where linear produced NaN); " + "3: forward_fill (missing botm filled from layer above by fix_missings_botms_and_min_layer_thickness); " + "4: shifted_for_min_thickness (botm shifted downward to enforce monotonically decreasing sequence)" +) +BOTM_METHOD_FLAG_VALUES = [0, 1, 2, 3, 4] + +KH_METHOD_MEANINGS = ( + "0: no_data (outside boundary polygon); " + "1: W_layer_polygon_kh (direct kh from area-weighted K polygon data); " + "2: S_layer_NHDZ_KD_linear (kh = KD/d, KD from linear interpolation of point data); " + "3: S_layer_NHDZ_KD_nearest (kh = KD/d, KD from nearest-neighbor interpolation fallback); " + "4: S_layer_Bergen_c_to_kh (kh = d*anisotropy/c, c from harmonic area-weighted polygon data); " + "5: S_layer_c_to_kh (kh = d*anisotropy/c, c from harmonic area-weighted polygon data, full extent); " + "6: fill_value (zero-thickness cell, set to fill_value_kh)" +) +KH_METHOD_FLAG_VALUES = [0, 1, 2, 3, 4, 5, 6] + +KV_METHOD_MEANINGS = ( + "0: no_data (outside boundary polygon); " + "1: W_layer_kh_anisotropy (kv = kh/anisotropy); " + "2: S_layer_d_over_c (kv = d/c, c from harmonic area-weighted polygon data); " + "3: fill_value (zero-thickness cell, set to fill_value_kv)" +) +KV_METHOD_FLAG_VALUES = [0, 1, 2, 3] + + +@pytest.mark.parametrize( + ("meanings", "flag_values", "expected"), + [ + pytest.param( + BOTM_METHOD_MEANINGS, + BOTM_METHOD_FLAG_VALUES, + [ + "no data", + "linear interpolation", + "nearest interpolation", + "forward fill", + "shifted for min thickness", + ], + id="botm_method", + ), + pytest.param( + KH_METHOD_MEANINGS, + KH_METHOD_FLAG_VALUES, + [ + "no data", + "W layer polygon kh", + "S layer NHDZ KD linear", + "S layer NHDZ KD nearest", + "S layer Bergen c to kh", + "S layer c to kh", + "fill value", + ], + id="kh_method", + ), + pytest.param( + KV_METHOD_MEANINGS, + KV_METHOD_FLAG_VALUES, + ["no data", "W layer kh anisotropy", "S layer d over c", "fill value"], + id="kv_method", + ), + ], +) +def test_parse_flag_labels_matches_flag_values_of_production_strings(meanings, flag_values, expected): + """Every production flag_meanings string yields one label per flag value. + + ``_plot_method_cross_section`` only applies the tick labels when + ``len(labels) == len(flag_values)``; a parser that produces one label too many or + too few silently leaves the colorbar showing raw integers instead of failing. + """ + labels = _parse_flag_labels(meanings) + + assert labels == expected + assert len(labels) == len(flag_values) + + +@pytest.mark.parametrize( + ("meanings", "expected"), + [ + # A trailing separator must not produce a phantom empty label, which would + # break the len(labels) == n_flags gate. + pytest.param("0: a; 1: b;", ["a", "b"], id="trailing-semicolon"), + pytest.param("0: a;;1: b", ["a", "b"], id="empty-entry"), + # No "N:" prefix: the whole entry is kept as the label. + pytest.param("plain_label", ["plain label"], id="no-flag-prefix"), + # Everything from the first "(" on is descriptive detail and is dropped. + pytest.param("7: short (long (nested) tail)", ["short"], id="parenthetical-dropped"), + # Multi-digit flag numbers are prefixes too, not part of the label. + pytest.param("10: ten_th_flag", ["ten th flag"], id="multi-digit-prefix"), + ], +) +def test_parse_flag_labels_adversarial_inputs(meanings, expected): + """Malformed or unusual flag_meanings entries are normalised, not mis-split.""" + assert _parse_flag_labels(meanings) == expected + + +# Cross-section line along the x-axis; ``project`` of a point therefore returns its +# x-coordinate exactly and ``distance`` returns |y|. +LINE = [(0.0, 0.0), (100.0, 0.0)] +BUFFER = 50.0 + + +@pytest.fixture +def botm_geojson(tmp_path): + """Write a four-point ``botm/botm.geojson`` with a shuffled, non-Range index. + + Perpendicular distances to ``LINE`` are 30, 80, 50 and 50 m, so with + ``BUFFER = 50`` the second point is excluded and the two points sitting exactly on + the buffer boundary are kept (the comparison is ``<=``). + """ + gdf = geopandas.GeoDataFrame( + { + # Column order deliberately differs from layer_names order (W11 precedes + # S13 there) and includes a column that is not a layer at all. + "S13": [100.0, 0.0, 4.0, -200.0], + "OBJECTID": [1, 2, 3, 4], + "W11": [10.0, 999.0, 5.0, -5.0], + }, + geometry=[Point(20.0, 30.0), Point(40.0, 80.0), Point(60.0, -50.0), Point(80.0, 50.0)], + index=[7, 3, 9, 1], + crs="EPSG:28992", + ) + fp = tmp_path / "botm" / "botm.geojson" + fp.parent.mkdir() + gdf.to_file(fp, driver="GeoJSON") + return tmp_path + + +def test_load_and_project_source_botm_selects_by_buffer_and_orders_by_layer_names(botm_geojson): + """Points are filtered by perpendicular distance and keyed in layer_names order.""" + result = _load_and_project_source_botm(data_path_2024=botm_geojson, line=LINE, buffer_distance=BUFFER) + + # Points 1, 3 and 4 survive; ``project`` onto the x-axis returns their x-coordinate. + np.testing.assert_array_equal(result["d_along"], [20.0, 60.0, 80.0]) + + # Keys follow layer_names order (W11 is index 0, S13 index 5), not the column order + # in the file. The values are positionally aligned with d_along, so any reordering + # of the keys without reordering the values would silently mis-pair z with x. + assert list(result["layers"]) == ["W11", "S13"] + assert layer_names.get_loc("W11") < layer_names.get_loc("S13") + np.testing.assert_array_equal(result["layers"]["W11"], [10.0, 5.0, -5.0]) + np.testing.assert_array_equal(result["layers"]["S13"], [100.0, 4.0, -200.0]) + + +def test_overlay_source_botm_draws_only_layers_with_in_window_points(botm_geojson): + """Only layers with at least one finite, in-window z become a scatter collection.""" + source_points = _load_and_project_source_botm(data_path_2024=botm_geojson, line=LINE, buffer_distance=BUFFER) + # An all-NaN layer is the case a layer model produces outside the PWN boundary. + source_points["layers"]["S11"] = np.full(3, np.nan) + + zmin, zmax = -120.0, 25.0 + ax = Figure().subplots() + _overlay_source_botm(ax, source_points, zmin=zmin, zmax=zmax) + + # W11 = [10, 5, -5] -> all three inside [-120, 25]; S13 = [100, 4, -200] -> only + # the middle point survives; S11 is all NaN -> no collection at all. + n_drawn_layers = 2 + assert len(ax.collections) == n_drawn_layers + np.testing.assert_array_equal(ax.collections[0].get_offsets(), [[20.0, 10.0], [60.0, 5.0], [80.0, -5.0]]) + np.testing.assert_array_equal(ax.collections[1].get_offsets(), [[60.0, 4.0]]) diff --git a/tests/test_pwnlayers_get_top.py b/tests/test_pwnlayers_get_top.py new file mode 100644 index 0000000..8807458 --- /dev/null +++ b/tests/test_pwnlayers_get_top.py @@ -0,0 +1,152 @@ +"""Tests for nhflotools.pwnlayers.layers.get_top_from_ahn. + +The function fills the gaps in the AHN surface in a fixed priority order: surface-water +peil first, then a constant over the North Sea, then interpolation from the remaining +valid cells. These tests pin that order, the full-coverage gate on the peil fill and the +(y, x) axis order handed to ``scipy.interpolate.griddata``. +""" + +import nlmod +import numpy as np +import pytest +import xarray as xr + +from nhflotools.pwnlayers.layers import get_top_from_ahn +from tests.util import make_rect_vertex_ds + + +def _ds_with_ahn(ahn, x=None, y=None): + """Build a one-row vertex ds carrying an ``ahn`` variable. + + With both replacement flags off, ``get_top_from_ahn`` reads nothing but ``ds['ahn']`` + and its ``x``/``y`` coordinates, so the cell centres may be moved off the regular + lattice to create a layout with unique nearest neighbours. + + Parameters + ---------- + ahn : sequence of float + AHN value per cell; NaN marks a gap to be filled. + x, y : sequence of float, optional + Replacement cell centres. Defaults to the regular grid of ``make_rect_vertex_ds``. + + Returns + ------- + xarray.Dataset + Vertex dataset with ``ahn`` on the ``icell2d`` dimension. + """ + ahn = np.asarray(ahn, dtype=float) + ds = make_rect_vertex_ds(nx=ahn.size, ny=1) + if x is not None: + ds = ds.assign_coords( + x=("icell2d", np.asarray(x, dtype=float)), + y=("icell2d", np.asarray(y, dtype=float)), + ) + ds["ahn"] = ("icell2d", ahn) + return ds + + +# Four donors and one gap, all off the x == y diagonal so that reversing one of the two +# coordinate tuples handed to griddata is detectable. Distances from the gap at +# (200, 600), in units of m**2: +# A (0, 800): 200**2 + 200**2 = 80000 <- true nearest, value 1.0 +# B (800, 0): 600**2 + 600**2 = 720000 +# C (1000, 600): 800**2 + 0**2 = 640000 +# D (200, 0): 0**2 + 600**2 = 360000 +# The implementation builds both donor and query tuples as (y, x). Reversing only the +# donors compares (x_d, y_d) against (600, 200) instead, giving +# A: 600**2 + 600**2 = 720000 B: 200**2 + 200**2 = 80000 <- swapped nearest, 2.0 +# C: 400**2 + 400**2 = 320000 D: 400**2 + 200**2 = 200000 +# so a one-sided axis swap moves the answer from 1.0 to 2.0. +_X = [0.0, 800.0, 1000.0, 200.0, 200.0] +_Y = [800.0, 0.0, 600.0, 0.0, 600.0] +_VALUES = [1.0, 2.0, 4.0, 8.0, 16.0] + + +def test_nearest_fill_is_euclidean_and_valid_cells_are_untouched(): + """The gap takes its Euclidean-nearest donor and valid cells pass through unchanged.""" + ahn = np.array(_VALUES, dtype=float) + ahn[4] = np.nan + + top = get_top_from_ahn( + _ds_with_ahn(ahn, x=_X, y=_Y), + replace_surface_water_with_peil=False, + replace_northsea_with_constant=None, + ) + + # Donor A at (0, 800) is nearest; the swapped metric would pick donor B at 2.0. + expected = np.array([1.0, 2.0, 4.0, 8.0, 1.0]) + np.testing.assert_array_equal(top.values, expected) + + # No gaps at all: griddata is handed an empty query set and the field is returned + # bit-identical (pins the scipy empty-xi edge). + full = get_top_from_ahn( + _ds_with_ahn(_VALUES, x=_X, y=_Y), + replace_surface_water_with_peil=False, + replace_northsea_with_constant=None, + ) + np.testing.assert_array_equal(full.values, np.array(_VALUES)) + + +def test_fill_priority_peil_then_sea_constant_then_nearest(monkeypatch): + """Peil wins over the sea constant, 0.0 still fills, partial cover falls through. + + Five cells in a row at x = 50, 150, 250, 350, 450 (100 m cells, area 10000 m2): + + ===== ======================================== ================================ + cell input expected + ===== ======================================== ================================ + 0 NaN, water over the whole cell + sea 2.0 (peil, not the sea constant) + 1 NaN, sea only 0.0 (falsy constant still fills) + 2 4.0 4.0 + 3 8.0 8.0 + 4 NaN, water over half the cell 8.0 (nearest donor, cell 3) + ===== ======================================== ================================ + + Cell 4 lies at the end of the row, so cell 3 at 100 m is its unique nearest donor + (cell 2 is 200 m away); its stage of 99.0 would surface instead if the full-coverage + gate were dropped. + """ + ds = _ds_with_ahn([np.nan, np.nan, 4.0, 8.0, np.nan]) + icell2d = ds["icell2d"] + seen = {} + + def fake_get_gdf_surface_water(extent=None, **_kwargs): + seen["extent"] = extent + return "gdf-sentinel" + + def fake_discretize_surface_water(_ds, gdf=None, **_kwargs): + seen["gdf"] = gdf + return xr.Dataset( + { + "rws_oppwater_area": ("icell2d", np.array([10000.0, 0.0, 0.0, 0.0, 5000.0])), + "rws_oppwater_stage": ("icell2d", np.array([2.0, np.nan, np.nan, np.nan, 99.0])), + }, + coords={"icell2d": icell2d}, + ) + + def fake_discretize_northsea(_ds, **_kwargs): + return xr.Dataset( + {"northsea": ("icell2d", np.array([True, True, False, False, False]))}, + coords={"icell2d": icell2d}, + ) + + # The rws readers download from live web services; they are the only seam mocked here. + monkeypatch.setattr(nlmod.read.rws, "get_gdf_surface_water", fake_get_gdf_surface_water) + monkeypatch.setattr(nlmod.read.rws, "discretize_surface_water", fake_discretize_surface_water) + monkeypatch.setattr(nlmod.read.rws, "discretize_northsea", fake_discretize_northsea) + + top = get_top_from_ahn( + ds, + replace_surface_water_with_peil=True, + replace_northsea_with_constant=0.0, + ) + + np.testing.assert_array_equal(top.values, np.array([2.0, 0.0, 4.0, 8.0, 8.0])) + assert seen["extent"] == ds.extent + assert seen["gdf"] == "gdf-sentinel" + + +def test_missing_ahn_raises_valueerror(): + """A ds without AHN is rejected up front rather than silently substituted.""" + with pytest.raises(ValueError, match="AHN"): + get_top_from_ahn(make_rect_vertex_ds(), replace_surface_water_with_peil=False) diff --git a/tests/test_well.py b/tests/test_well.py index 7b06408..cf2199c 100644 --- a/tests/test_well.py +++ b/tests/test_well.py @@ -4,6 +4,7 @@ import geopandas as gpd import numpy as np +import pandas as pd import pytest import xarray as xr from shapely.geometry import Point @@ -319,3 +320,131 @@ def test_tata_fresh_wells_preserve_no_qualifying_kd_failure(monkeypatch, tmp_pat with pytest.raises(IndexError): well.get_wells_tata_dataframes(tmp_path, ds) + + +# --- get_wells_pwn_dataframe ------------------------------------------------ +# +# One synthetic secundair bookkeeping, shared by every PWN test below. +# +# Wells (``sec_nput`` = number of wells the secundair flow is spread over): +# W1/W2/W3 tag T1 sec_nput 3 -> extraction, split three ways +# W4 tag T2 sec_nput 1 -> extraction, asymmetric series +# W5 tag TX sec_nput 2 -> tag absent from the feather (drop) +# W6 tag T1 sec_nput 0 -> division by zero -> +/-inf (drop) +# W7 tag T0 sec_nput 2 -> secundair median is zero (drop) +# W8 tag TINF sec_nput 2 -> infiltration, positive median (keep) +# +# Series medians, by hand (median, not mean -- T2 and TINF are asymmetric so a +# mean would give a different answer): +# T1 [-30, -30, -30] -> -30 (mean -30, not a discriminator) +# T2 [-10, -20, -60] -> -20 (mean -30) +# T0 [ -4, 0, 4] -> 0 (drops on the ``Q != 0`` mask) +# TINF [ 4, 8, 20] -> 8 (mean 32/3) +# The 'ophaal tijdstip' column is datetime: it must be excluded by +# ``numeric_only=True`` and must not become a mappable secundair tag. +_PWN_WELLS = ( + # locatie, sec_flow_tag, sec_nput, x, y + ("W1", "T1", 3, 0.0, 1.0), + ("W2", "T1", "3", 100.0, 2.0), # string on purpose; GeoJSON stringifies the column anyway + ("W3", "T1", 3, 200.0, 3.0), + ("W4", "T2", 1, 300.0, 4.0), + ("W5", "TX", 2, 400.0, 5.0), + ("W6", "T1", 0, 500.0, 6.0), + ("W7", "T0", 2, 600.0, 7.0), + ("W8", "TINF", 2, 700.0, 8.0), +) +_PWN_FLOWS = { + "T1": [-30.0, -30.0, -30.0], + "T2": [-10.0, -20.0, -60.0], + "T0": [-4.0, 0.0, 4.0], + "TINF": [4.0, 8.0, 20.0], +} +_PWN_MEDIANS = {"T1": -30.0, "T2": -20.0, "T0": 0.0, "TINF": 8.0} + + +@pytest.fixture +def pwn_data_path(tmp_path): + """Write a real ``pumping_infiltration_wells.geojson`` and ``sec_flows.feather``. + + Real files rather than patched readers: the GeoJSON round-trip is what turns + ``sec_nput`` into strings, which is exactly the input ``well.py`` has to coerce. + + Returns + ------- + pathlib.Path + Directory holding both input files. + """ + wells = gpd.GeoDataFrame( + { + "locatie": [row[0] for row in _PWN_WELLS], + "sec_flow_tag": [row[1] for row in _PWN_WELLS], + "sec_nput": [row[2] for row in _PWN_WELLS], + }, + geometry=[Point(row[3], row[4]) for row in _PWN_WELLS], + crs="EPSG:28992", + ) + wells.to_file(tmp_path / "pumping_infiltration_wells.geojson", driver="GeoJSON") + + flows = pd.DataFrame({ + **_PWN_FLOWS, + "ophaal tijdstip": pd.to_datetime(["2024-01-01", "2024-01-02", "2024-01-03"]), + }) + flows.to_feather(tmp_path / "sec_flows.feather") + return tmp_path + + +def test_pwn_q_splits_secundair_median_and_conserves_the_total(pwn_data_path): + """Per-well Q is 24 * median / sec_nput, so a secundair sums back to 24 * median.""" + wdf = well.get_wells_pwn_dataframe(pwn_data_path) + + # Surviving wells, in input order; W5/W6/W7 have no usable flow. + assert wdf.index.to_list() == ["W1", "W2", "W3", "W4", "W8"] + assert wdf.index.name == "locatie" + + # m3/h -> m3/day is a factor 24; the secundair flow is split over sec_nput wells. + assert wdf.loc[["W1", "W2", "W3"], "Q"].to_list() == [-240.0] * 3 # 24 * -30 / 3 + assert wdf.loc["W4", "Q"] == 24.0 * _PWN_MEDIANS["T2"] / 1 # -480.0, median not mean + assert wdf.loc["W8", "Q"] == 24.0 * _PWN_MEDIANS["TINF"] / 2 # +96.0 + + # The WEL mass-balance contract: the three T1 wells reconstruct the whole secundair. + assert wdf.loc[["W1", "W2", "W3"], "Q"].sum() == 24.0 * _PWN_MEDIANS["T1"] + + # "3" arrived as a string from the GeoJSON and must have been coerced to a number. + assert pd.api.types.is_numeric_dtype(wdf["sec_nput"]) + assert wdf["sec_nput"].to_list() == [3, 3, 3, 1, 2] + + # Geometry-derived and constant MAW/transport columns. + assert wdf["x"].to_list() == [0.0, 100.0, 200.0, 300.0, 700.0] + assert wdf["y"].to_list() == [1.0, 2.0, 3.0, 4.0, 8.0] + assert wdf["rw"].to_list() == [0.25] * 5 + assert wdf["CONCENTRATION"].to_list() == [0.0] * 5 + + +def test_pwn_drops_unusable_wells_warns_once_and_keeps_infiltration(pwn_data_path, caplog): + """Unmapped, zero-nput and zero-median wells drop; a positive infiltration survives.""" + with caplog.at_level("WARNING", logger=well.logger.name): + wdf = well.get_wells_pwn_dataframe(pwn_data_path) + + assert set(wdf.index) == {"W1", "W2", "W3", "W4", "W8"} + + # One summary warning carrying (n_dropped, n_total): 3 of the 8 input wells. + warnings = [record for record in caplog.records if "without a nonzero secundair flow" in record.getMessage()] + assert len(warnings) == 1 + assert warnings[0].args == (3, len(_PWN_WELLS)) + + # The drop mask is sign-symmetric: it removes non-finite and zero, never negatives, + # and never the sole infiltration well. + assert np.isfinite(wdf["Q"]).all() + assert (wdf["Q"] != 0.0).all() + assert wdf.loc["W8", "Q"] > 0.0 + assert (wdf.loc[["W1", "W2", "W3", "W4"], "Q"] < 0.0).all() + + +@pytest.mark.parametrize( + ("flow_product", "expected_error"), + [("timeseries", NotImplementedError), ("bogus", ValueError)], +) +def test_pwn_unsupported_flow_product_raises(pwn_data_path, flow_product, expected_error): + """Only the median product is implemented; other products fail loudly.""" + with pytest.raises(expected_error): + well.get_wells_pwn_dataframe(pwn_data_path, flow_product=flow_product) diff --git a/tests/util.py b/tests/util.py new file mode 100644 index 0000000..5b62dc8 --- /dev/null +++ b/tests/util.py @@ -0,0 +1,203 @@ +"""Builders for tiny synthetic model datasets. + +nlmod's own ``tests/util.py`` is not importable from an installed nlmod, so the few +helpers the nhflotools suite needs are built here. Everything is hand-built: no gridgen, +no MODFLOW binaries and no network, which keeps the unit tests to milliseconds. +""" + +import flopy +import geopandas as gpd +import nlmod +import numpy as np +import xarray as xr +from shapely.geometry import box + + +def make_rect_vertex_ds( + nx=2, + ny=2, + delr=100.0, + botm=(-10.0, -20.0), + top=0.0, + kh=5.0, + transport=0, +): + """Build a rectangular vertex (DISV) model dataset. + + Cells are numbered row-major from the north-west corner, matching MODFLOW's + convention. The grid is regular, so every cell area is ``delr**2`` exactly. + + Parameters + ---------- + nx, ny : int + Number of cells in the x and y direction. + delr : float + Cell size [m]; cells are square. + botm : sequence of float + Layer bottom elevations [mNAP], one per layer. + top : float + Surface elevation [mNAP]. + kh : float + Horizontal conductivity [m/day]. + transport : int + Value of the ``transport`` attribute nlmod's package builders read. + + Returns + ------- + xarray.Dataset + Vertex dataset with ``top``, ``botm``, ``kh``, ``kv``, ``area`` and ``idomain``, + the ``xv``/``yv``/``icvert`` grid geometry and the attributes nlmod requires. + """ + botm = np.asarray(botm, dtype=float) + nlay = botm.size + ncell = nx * ny + extent = [0.0, nx * delr, 0.0, ny * delr] + + # Cell centres, row-major from the north-west corner. + ix, iy = np.meshgrid(np.arange(nx), np.arange(ny)) + x = (ix.ravel() + 0.5) * delr + y = extent[3] - (iy.ravel() + 0.5) * delr + + # Vertices of the (nx+1) x (ny+1) lattice, numbered the same way. + vx, vy = np.meshgrid(np.arange(nx + 1) * delr, extent[3] - np.arange(ny + 1) * delr) + vx, vy = vx.ravel(), vy.ravel() + + def vertex_id(row, col): + return row * (nx + 1) + col + + icvert = np.array( + [ + [ + vertex_id(r, c), + vertex_id(r, c + 1), + vertex_id(r + 1, c + 1), + vertex_id(r + 1, c), + ] + for r in range(ny) + for c in range(nx) + ], + dtype=int, + ) + + ds = xr.Dataset( + data_vars={ + "top": ("icell2d", np.full(ncell, float(top))), + "botm": (("layer", "icell2d"), np.tile(botm[:, None], (1, ncell))), + "kh": (("layer", "icell2d"), np.full((nlay, ncell), float(kh))), + "kv": (("layer", "icell2d"), np.full((nlay, ncell), float(kh) / 10.0)), + "area": ("icell2d", np.full(ncell, delr**2)), + "idomain": (("layer", "icell2d"), np.ones((nlay, ncell), dtype=int)), + "xv": ("iv", vx), + "yv": ("iv", vy), + "icvert": (("icell2d", "nvert"), icvert), + }, + coords={ + "layer": np.arange(nlay), + "icell2d": np.arange(ncell), + "iv": np.arange(vx.size), + "x": ("icell2d", x), + "y": ("icell2d", y), + }, + attrs={ + "gridtype": "vertex", + "extent": extent, + "model_name": "test", + "mfversion": "mf6", + "model_ws": ".", + "transport": transport, + }, + ) + ds["icvert"].attrs["nodata"] = -1 + return ds + + +def add_time(ds, start="2022-01-01"): + """Add a single steady-state stress period, required before building sim/tdis.""" + return nlmod.time.set_ds_time(ds, start=start, time=[1.0], steady=True) + + +def make_gwf_disv(ds, model_ws): + """Build an in-memory sim/gwf/disv from a vertex ds. Never written, never run.""" + ds = add_time(ds) + ds.attrs["model_ws"] = str(model_ws) + # Packages are only built in memory, never written or run, so the exe need not exist. + sim = nlmod.sim.sim(ds, exe_name="mf6") + nlmod.sim.tdis(ds, sim) + gwf = nlmod.gwf.gwf(ds, sim) + nlmod.gwf.disv(ds, gwf) + return ds, gwf + + +def cell_polygon(ds, icell2d): + """Return the square polygon of one cell of a ``make_rect_vertex_ds`` grid.""" + delr = np.sqrt(float(ds["area"].isel(icell2d=icell2d))) + x = float(ds["x"].isel(icell2d=icell2d)) + y = float(ds["y"].isel(icell2d=icell2d)) + return box(x - delr / 2, y - delr / 2, x + delr / 2, y + delr / 2) + + +def make_gdf(geometries, crs="EPSG:28992", **columns): + """Build a GeoDataFrame from geometries plus scalar-or-sequence columns.""" + n = len(geometries) + data = {k: (v if isinstance(v, (list, tuple, np.ndarray)) else [v] * n) for k, v in columns.items()} + return gpd.GeoDataFrame(data, geometry=list(geometries), crs=crs) + + +def write_mf6_listing(path, pct_disc, budgetkey="VOLUME BUDGET FOR ENTIRE MODEL"): + """Write a minimal MF6 listing file that flopy's Mf6ListBudget can parse. + + The layout mirrors a real MF6 listing: a volume budget block with matching IN/OUT + entries, the percent-discrepancy line, and the time-summary block flopy needs to + attach a time index to the budget. + """ + inflow, outflow = 100.0, 100.0 * (1.0 - pct_disc / 100.0) + text = f""" + {budgetkey} AT END OF TIME STEP 1, STRESS PERIOD 1 + ------------------------------------------------------------------------------ + + CUMULATIVE VOLUME L**3 RATES FOR THIS TIME STEP L**3/T + ------------------ ------------------------ + + IN: IN: + --- --- + CHD = 0.0000 CHD ={inflow:>17.4f} + + TOTAL IN = 0.0000 TOTAL IN ={inflow:>17.4f} + + OUT: OUT: + ---- ---- + CHD = 0.0000 CHD ={outflow:>17.4f} + + TOTAL OUT = 0.0000 TOTAL OUT ={outflow:>17.4f} + + IN - OUT = 0.0000 IN - OUT ={inflow - outflow:>17.4f} + + PERCENT DISCREPANCY = 0.00 PERCENT DISCREPANCY ={pct_disc:>17.2f} + + + TIME SUMMARY AT END OF TIME STEP 1 IN STRESS PERIOD 1 + SECONDS MINUTES HOURS DAYS YEARS + ----------------------------------------------------------- + TIME STEP LENGTH 86400. 1440.0 24.000 1.0000 2.73785E-03 + STRESS PERIOD TIME 86400. 1440.0 24.000 1.0000 2.73785E-03 + TOTAL TIME 86400. 1440.0 24.000 1.0000 2.73785E-03 +""" + path.write_text(text) + return path + + +def make_structured_ds(extent=(0.0, 200.0, 0.0, 200.0), delr=100.0, **kwargs): + """Tiny structured ds via nlmod's own builder, with executable download suppressed.""" + return nlmod.get_ds(list(extent), delr=delr, download_exe=False, **kwargs) + + +__all__ = [ + "add_time", + "cell_polygon", + "flopy", + "make_gdf", + "make_gwf_disv", + "make_rect_vertex_ds", + "make_structured_ds", + "write_mf6_listing", +] From 9f93ef43868b8558b45806b70f59b00b9dd0e7f9 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 20 Jul 2026 14:15:59 +0200 Subject: [PATCH 2/5] Suppress NHFLODATA_LOCATION for the whole session, not per test The variable was unset by the autouse hygiene fixture, which is function scoped and therefore runs after collection and after any session- or module-scoped fixture. Nothing resolves data paths that early today, so the suite did pass with the variable set, but the guard sat one scope too low to be relied on: a parametrisation or a session fixture that resolved a dataset path would have read the developer's mount instead of the packaged mockup. Moving it to pytest_configure closes that: it fires before collection, so every later mechanism sees the variable unset. pytest_unconfigure puts it back for in-process runners that outlive the session. get_abs_data_path resolves against the variable whenever it is set and only warns when the result is missing, so without this the module would pass or fail on whichever datasets happen to be mounted locally rather than on what NHFLO/data ships. Adds test_data_location_env_is_suppressed_for_the_session so removing the suppression fails with the cause named, instead of 41 opaque path failures. Verified by running the suite with the variable unset, set to an existing empty directory, and set to the real mockup root, and by removing the suppression with the variable set: 42 fail, led by the new test. --- README.md | 6 +++++- TEST_PLAN.md | 12 +++++++++-- tests/conftest.py | 36 ++++++++++++++++++++++++++------ tests/test_nhflodata_contract.py | 17 +++++++++++++-- 4 files changed, 60 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index 26355ca..7b6e926 100644 --- a/README.md +++ b/README.md @@ -45,12 +45,16 @@ door 09pwnmodel2 en zijn ongetest — behandel ze als legacy of werk-in-uitvoeri ## Tests -De suite telt 126 tests en draait in ongeveer 4 seconden. De tests gebruiken uitsluitend +De suite telt 127 tests en draait in ongeveer 4 seconden. De tests gebruiken uitsluitend kleine synthetische modellen — geen netwerk, geen gridgen en geen grote datasets — zodat ze op elke pull request in CI meedraaien. Alle live webservices (HHNK, REGIS, RWS, achtergrondkaarten) worden per naam gemonkeypatcht; de rest van nlmod draait echt, waardoor de suite meteen dienstdoet als compatibiliteitscanary voor `nlmod@dev`. +`NHFLODATA_LOCATION` wordt voor de duur van de testsessie uitgezet en daarna hersteld, zodat +de tests altijd tegen de meegeleverde mockup-data draaien — ook op een machine waar een +echte datamap is gekoppeld. + De enige uitzondering is één test met de marker `mf6`: die draait MODFLOW echt op een model van 3x3x2 cellen en controleert de waterbalans plus een analytische oplossing. De binaries worden zo nodig automatisch door nlmod gedownload en daarna hergebruikt. diff --git a/TEST_PLAN.md b/TEST_PLAN.md index 4023d53..6c173cd 100644 --- a/TEST_PLAN.md +++ b/TEST_PLAN.md @@ -611,7 +611,15 @@ Status legend: ☐ not started · ◐ in progress · ☑ done and green. ### Test files -All green. Suite total: **125 passed, 1 xfailed in 4.2 s**; slowest single test 0.59 s. +All green. Suite total: **126 passed, 1 xfailed in ~4 s**; slowest single test 0.59 s. + +The suite is independent of the developer's environment: `conftest.pytest_configure` unsets +`NHFLODATA_LOCATION` before collection (so parametrisation and any higher-scoped fixture +also see it unset) and `pytest_unconfigure` restores it. Verified by running the whole suite +three ways — variable unset, set to an existing empty directory, and set to the real mockup +root — all green, and by mutation: with the suppression removed and the variable set, 42 +tests fail, led by `test_data_location_env_is_suppressed_for_the_session`, which names the +cause rather than leaving 41 opaque path failures. | File | Plan § | Tests | Status | |---|---|---|---| @@ -624,7 +632,7 @@ All green. Suite total: **125 passed, 1 xfailed in 4.2 s**; slowest single test | `test_pwnlayers_get_top.py` | 4.7 | 3 | ☑ | | `test_postprocessing.py` (extended) | 4.6 | 17 | ☑ | | `test_pwnlayers3_plot.py` | 4.9 | 10 | ☑ | -| `test_nhflodata_contract.py` | 4.10 | 41 | ☑ | +| `test_nhflodata_contract.py` | 4.10 | 42 | ☑ | | `test_mf6_smoke.py` | 4.11 | 1 | ☑ | Every test file was mutation-checked during implementation: a plausible bug was introduced diff --git a/tests/conftest.py b/tests/conftest.py index eae2b35..2adc8df 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,24 +1,48 @@ """Shared fixtures for the nhflotools test suite.""" +import os + import matplotlib.pyplot as plt import pytest import xarray as xr from .util import make_gwf_disv, make_rect_vertex_ds +_DATA_LOCATION_ENV = "NHFLODATA_LOCATION" +_saved_data_location = None -@pytest.fixture(autouse=True) -def _hygiene(monkeypatch): - """Keep tests independent: no leaked figures, open files or data-location env var. - ``NHFLODATA_LOCATION`` is deleted so a developer's local data mount can never - redirect the tests away from the mockup data that CI resolves. +def pytest_configure(config): # noqa: ARG001 + """Suppress ``NHFLODATA_LOCATION`` for the whole session, restoring it afterwards. + + ``nhflodata.get_abs_data_path`` resolves against that variable when it is set and only + *warns* when the resulting path is missing, so a developer with a real data mount would + otherwise run these tests against their own data instead of the mockup data shipped in + the wheel -- the same data CI resolves. The tests would then pass or fail depending on + which datasets happen to be mounted locally. + + This runs in ``pytest_configure`` rather than in a fixture deliberately: it fires before + collection, so parametrisation and any session- or module-scoped fixture also see the + variable unset. A function-scoped fixture would run too late for both. + """ + global _saved_data_location # noqa: PLW0603 + _saved_data_location = os.environ.pop(_DATA_LOCATION_ENV, None) + + +def pytest_unconfigure(config): # noqa: ARG001 + """Put ``NHFLODATA_LOCATION`` back, for in-process runners that outlive the session.""" + if _saved_data_location is not None: + os.environ[_DATA_LOCATION_ENV] = _saved_data_location + + +@pytest.fixture(autouse=True) +def _hygiene(): + """Keep tests independent: no figures or open netCDF handles leak between tests. ``FILE_CACHE.clear()`` already closes the netCDF handles xarray holds open; an additional ``gc.collect()`` here cost ~50 ms per test (two thirds of the whole suite's runtime) and closed nothing further, so it is deliberately absent. """ - monkeypatch.delenv("NHFLODATA_LOCATION", raising=False) yield plt.close("all") xr.backends.file_manager.FILE_CACHE.clear() diff --git a/tests/test_nhflodata_contract.py b/tests/test_nhflodata_contract.py index a0397a5..2777132 100644 --- a/tests/test_nhflodata_contract.py +++ b/tests/test_nhflodata_contract.py @@ -4,10 +4,11 @@ restructure of the NHFLO/data repository would otherwise surface as a ``FileNotFoundError`` deep inside a model run. These tests stat every (dataset, relative file) pair that nhflotools -- or ``modelscripts/09pwnmodel2/01_pwnmodel2.py`` -- hardcodes, against the mockup data -packaged with nhflodata (the autouse hygiene fixture deletes ``NHFLODATA_LOCATION``, so -resolution always lands on the mockup). +packaged with nhflodata (``conftest.pytest_configure`` unsets ``NHFLODATA_LOCATION`` for the +whole session, so resolution always lands on the mockup). """ +import os from pathlib import Path import pandas as pd @@ -94,6 +95,18 @@ def test_mockup_data_file_resolves_and_is_nonempty(dataset, relative_path): assert path.stat().st_size > 0, f"{dataset}/{relative_path} is empty" +def test_data_location_env_is_suppressed_for_the_session(): + """Every path check below assumes resolution lands on the packaged mockup data. + + ``get_abs_data_path`` resolves against ``NHFLODATA_LOCATION`` whenever it is set, and + only warns when the result is missing. On a machine with a real data mount that would + silently point this whole module at local data, so it would pass or fail on whichever + datasets happen to be mounted rather than on what NHFLO/data ships. This asserts the + precondition ``conftest.pytest_configure`` establishes. + """ + assert os.environ.get("NHFLODATA_LOCATION") is None + + def test_koppeltabel_columns_and_layer_coverage(): """The koppeltabel exposes the columns and the layer names the merge indexes by. From 269ac23b5e55b8c42eb3205371daa2b60c1b2ec8 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Mon, 20 Jul 2026 14:32:38 +0200 Subject: [PATCH 3/5] Set NHFLODATA_LOCATION empty via pytest-env instead of a conftest hook pytest-env was already a test dependency and already carried MPLBACKEND, so the data location belongs there too: env = ["MPLBACKEND=Agg", "NHFLODATA_LOCATION="] pytest-env overrides an inherited value with the empty string, and get_paths reads that as "use mockup" (get_paths.py:69 defaults the variable to "" and :96 branches on it being falsy), so empty resolves exactly as unset does. Being an ini setting it also applies before collection, which is what the previous pytest_configure hook was there to guarantee -- one declarative line replaces thirty of conftest. The guard test now asserts the resolved path lies under the installed nhflodata/data/mockup rather than asserting the variable is None. That checks the outcome instead of the mechanism, so it keeps holding whichever way the suppression is implemented, and it also fails if nhflodata ever stops treating an empty value as "use mockup" -- which the previous test could not detect. Verified with the variable unset, set to an existing empty directory, and set to a nonexistent path: green in all three. With the pytest-env entry removed and the variable set, 42 tests fail, led by the guard test. --- README.md | 7 ++++--- TEST_PLAN.md | 27 ++++++++++++++++++++------- pyproject.toml | 2 +- tests/conftest.py | 28 ---------------------------- tests/test_nhflodata_contract.py | 26 ++++++++++++++++---------- 5 files changed, 41 insertions(+), 49 deletions(-) diff --git a/README.md b/README.md index 7b6e926..eab39fb 100644 --- a/README.md +++ b/README.md @@ -51,9 +51,10 @@ ze op elke pull request in CI meedraaien. Alle live webservices (HHNK, REGIS, RW achtergrondkaarten) worden per naam gemonkeypatcht; de rest van nlmod draait echt, waardoor de suite meteen dienstdoet als compatibiliteitscanary voor `nlmod@dev`. -`NHFLODATA_LOCATION` wordt voor de duur van de testsessie uitgezet en daarna hersteld, zodat -de tests altijd tegen de meegeleverde mockup-data draaien — ook op een machine waar een -echte datamap is gekoppeld. +`NHFLODATA_LOCATION` wordt in `pyproject.toml` via pytest-env op leeg gezet, zodat de tests +altijd tegen de meegeleverde mockup-data draaien — ook op een machine waar een echte +datamap is gekoppeld. Dat geldt alleen binnen pytest; een modelscript gebruikt gewoon jouw +eigen datamap. De enige uitzondering is één test met de marker `mf6`: die draait MODFLOW echt op een model van 3x3x2 cellen en controleert de waterbalans plus een analytische oplossing. De diff --git a/TEST_PLAN.md b/TEST_PLAN.md index 6c173cd..80d6eca 100644 --- a/TEST_PLAN.md +++ b/TEST_PLAN.md @@ -613,13 +613,26 @@ Status legend: ☐ not started · ◐ in progress · ☑ done and green. All green. Suite total: **126 passed, 1 xfailed in ~4 s**; slowest single test 0.59 s. -The suite is independent of the developer's environment: `conftest.pytest_configure` unsets -`NHFLODATA_LOCATION` before collection (so parametrisation and any higher-scoped fixture -also see it unset) and `pytest_unconfigure` restores it. Verified by running the whole suite -three ways — variable unset, set to an existing empty directory, and set to the real mockup -root — all green, and by mutation: with the suppression removed and the variable set, 42 -tests fail, led by `test_data_location_env_is_suppressed_for_the_session`, which names the -cause rather than leaving 41 opaque path failures. +The suite is independent of the developer's environment. `pyproject.toml` sets +`NHFLODATA_LOCATION=` through pytest-env, alongside `MPLBACKEND=Agg`: + +```toml +env = ["MPLBACKEND=Agg", "NHFLODATA_LOCATION="] +``` + +pytest-env overrides an inherited value with the empty string, and `get_paths.py:69` +defaults that variable to `""` while `:96` branches on it being falsy — so empty resolves to +the packaged mockup exactly as unset does. Being an ini setting it applies before collection, +so parametrisation and every fixture scope see it too. + +`test_paths_resolve_to_the_packaged_mockup` asserts the resolved path lies under the +installed `nhflodata/data/mockup`, i.e. the outcome rather than the variable, so it stays +honest if that empty-means-mockup contract ever changes. + +Verified by running the whole suite with the variable unset, set to an existing empty +directory, and set to a nonexistent path — all green — and by mutation: with the pytest-env +entry removed and the variable set, 42 tests fail, led by that guard test, which names the +cause instead of leaving 41 opaque path failures. | File | Plan § | Tests | Status | |---|---|---|---| diff --git a/pyproject.toml b/pyproject.toml index 9220bb5..ba46ae7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -96,7 +96,7 @@ test = "pytest -v" [tool.pytest.ini_options] addopts = "--strict-markers" -env = ["MPLBACKEND=Agg"] +env = ["MPLBACKEND=Agg", "NHFLODATA_LOCATION="] markers = [ "mf6: requires the MODFLOW 6 executable", "network: hits live web services; excluded in CI", diff --git a/tests/conftest.py b/tests/conftest.py index 2adc8df..7d1f65b 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,39 +1,11 @@ """Shared fixtures for the nhflotools test suite.""" -import os - import matplotlib.pyplot as plt import pytest import xarray as xr from .util import make_gwf_disv, make_rect_vertex_ds -_DATA_LOCATION_ENV = "NHFLODATA_LOCATION" -_saved_data_location = None - - -def pytest_configure(config): # noqa: ARG001 - """Suppress ``NHFLODATA_LOCATION`` for the whole session, restoring it afterwards. - - ``nhflodata.get_abs_data_path`` resolves against that variable when it is set and only - *warns* when the resulting path is missing, so a developer with a real data mount would - otherwise run these tests against their own data instead of the mockup data shipped in - the wheel -- the same data CI resolves. The tests would then pass or fail depending on - which datasets happen to be mounted locally. - - This runs in ``pytest_configure`` rather than in a fixture deliberately: it fires before - collection, so parametrisation and any session- or module-scoped fixture also see the - variable unset. A function-scoped fixture would run too late for both. - """ - global _saved_data_location # noqa: PLW0603 - _saved_data_location = os.environ.pop(_DATA_LOCATION_ENV, None) - - -def pytest_unconfigure(config): # noqa: ARG001 - """Put ``NHFLODATA_LOCATION`` back, for in-process runners that outlive the session.""" - if _saved_data_location is not None: - os.environ[_DATA_LOCATION_ENV] = _saved_data_location - @pytest.fixture(autouse=True) def _hygiene(): diff --git a/tests/test_nhflodata_contract.py b/tests/test_nhflodata_contract.py index 2777132..f7e4dfa 100644 --- a/tests/test_nhflodata_contract.py +++ b/tests/test_nhflodata_contract.py @@ -4,11 +4,10 @@ restructure of the NHFLO/data repository would otherwise surface as a ``FileNotFoundError`` deep inside a model run. These tests stat every (dataset, relative file) pair that nhflotools -- or ``modelscripts/09pwnmodel2/01_pwnmodel2.py`` -- hardcodes, against the mockup data -packaged with nhflodata (``conftest.pytest_configure`` unsets ``NHFLODATA_LOCATION`` for the -whole session, so resolution always lands on the mockup). +packaged with nhflodata (``pyproject.toml`` sets ``NHFLODATA_LOCATION=`` through pytest-env, +so resolution always lands on the mockup). """ -import os from pathlib import Path import pandas as pd @@ -95,16 +94,23 @@ def test_mockup_data_file_resolves_and_is_nonempty(dataset, relative_path): assert path.stat().st_size > 0, f"{dataset}/{relative_path} is empty" -def test_data_location_env_is_suppressed_for_the_session(): +def test_paths_resolve_to_the_packaged_mockup(): """Every path check below assumes resolution lands on the packaged mockup data. - ``get_abs_data_path`` resolves against ``NHFLODATA_LOCATION`` whenever it is set, and - only warns when the result is missing. On a machine with a real data mount that would - silently point this whole module at local data, so it would pass or fail on whichever - datasets happen to be mounted rather than on what NHFLO/data ships. This asserts the - precondition ``conftest.pytest_configure`` establishes. + ``get_abs_data_path`` resolves against ``NHFLODATA_LOCATION`` whenever it holds a + non-empty value, and only warns when the result is missing. On a machine with a real + data mount that would silently point this whole module at local data, so it would pass + or fail on whichever datasets happen to be mounted rather than on what NHFLO/data ships. + + ``pyproject.toml`` therefore sets ``NHFLODATA_LOCATION=`` through pytest-env, which + overrides any inherited value with the empty string that ``get_paths`` reads as "use + mockup" (``get_paths.py:69`` defaults to ``""``, ``:96`` branches on it being falsy). + Asserting the resolved location rather than the variable keeps this honest if that + empty-means-mockup contract ever changes. """ - assert os.environ.get("NHFLODATA_LOCATION") is None + mockup_root = Path(get_paths.__file__).parent / "data" / "mockup" + resolved = Path(get_abs_data_path(name="bodemlagen_pwn_regis_koppeltabel", version="latest")) + assert mockup_root in resolved.parents def test_koppeltabel_columns_and_layer_coverage(): From 1258d8770cf6c0ed90af2aaafe688590d5b0938a Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Wed, 22 Jul 2026 12:11:27 +0200 Subject: [PATCH 4/5] Pass source_botm_name='bottom' now nlmod defaults it to 'botm' nlmod@dev gave aggregate_by_weighted_mean_to_ds configurable source layer names, defaulting the floor to 'botm'; the NHI chloride file names it 'bottom'. Caught by the eight test_nhi_chloride tests on CI. --- src/nhflotools/nhi_chloride.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/nhflotools/nhi_chloride.py b/src/nhflotools/nhi_chloride.py index fb1c049..53b3aca 100644 --- a/src/nhflotools/nhi_chloride.py +++ b/src/nhflotools/nhi_chloride.py @@ -37,7 +37,8 @@ def get_nhi_chloride_concentration(ds, data_path_nhi_chloride): # cli has x and y of ds but layer of cl cli = cl.interp(x=ds.x, y=ds.y, method="nearest").drop_vars(["dy", "dx", "percentile"], errors="ignore") - da = nlmod.layers.aggregate_by_weighted_mean_to_ds(ds, xr.Dataset({"p50": cli}), "p50") + # The NHI file names its voxel floor "bottom"; nlmod's default is "botm". + da = nlmod.layers.aggregate_by_weighted_mean_to_ds(ds, xr.Dataset({"p50": cli}), "p50", source_botm_name="bottom") da.values[0] = xr.where(ds["northsea"] == 1, SEA_CHLORIDE_MG_L, da.values[0]) From 8caf0e32e7f5a3e9afb1bb78be71d0e7683930db Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Wed, 22 Jul 2026 12:11:27 +0200 Subject: [PATCH 5/5] =?UTF-8?q?Correct=20TEST=5FPLAN's=20stale=20header:?= =?UTF-8?q?=20the=20plan=20is=20implemented=20(=C2=A79)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- TEST_PLAN.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TEST_PLAN.md b/TEST_PLAN.md index 80d6eca..6c54b78 100644 --- a/TEST_PLAN.md +++ b/TEST_PLAN.md @@ -1,6 +1,6 @@ # Test plan — nhflotools, scoped to the 09pwnmodel2 closure -Implementation plan only; nothing here is implemented yet. Scope = the 13 nhflotools entry +Implemented in full — §9 records the progress and the deviations. Scope = the 13 nhflotools entry points `models/modelscripts/09pwnmodel2/01_pwnmodel2.py` imports, plus their transitive nhflotools closure (`pwnlayers.merge_layer_models`, `pwnlayers.utils`, `panden.get_oppervlakte_pwn_shapes`). Target: a lean, meaningful pytest suite that runs