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..eab39fb 100644 --- a/README.md +++ b/README.md @@ -4,4 +4,66 @@ 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 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 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 +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..6c54b78 --- /dev/null +++ b/TEST_PLAN.md @@ -0,0 +1,695 @@ +# Test plan — nhflotools, scoped to the 09pwnmodel2 closure + +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 +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: **126 passed, 1 xfailed in ~4 s**; slowest single test 0.59 s. + +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 | +|---|---|---|---| +| `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 | 42 | ☑ | +| `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..ba46ae7 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", "NHFLODATA_LOCATION="] +markers = [ + "mf6: requires the MODFLOW 6 executable", + "network: hits live web services; excluded in CI", +] testpaths = ["tests"] 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]) diff --git a/tests/conftest.py b/tests/conftest.py index e92f4f9..9102ce6 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,16 +1,44 @@ -"""Shared synthetic disv fixtures for the lake/polder tests. - -These build a tiny vertex (disv) MODFLOW 6 model plus a matching xarray dataset, so the -lake carve/stage/mask helpers can be exercised without a solver, real PWN data, or a -network connection. -""" +"""Shared fixtures for the nhflotools test suite.""" import flopy +import matplotlib.pyplot as plt import numpy as np import pytest import xarray as xr from shapely.geometry import box +from tests.util import make_gwf_disv, make_rect_vertex_ds + + +@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. + """ + yield + plt.close("all") + xr.backends.file_manager.FILE_CACHE.clear() + + +@pytest.fixture +def vertex_ds(): + """Build 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) + def _build_disv(nrow=3, ncol=3, top=5.0, dz=10.0, transport=0): """Build an independent (ds, gwf, cell_geometries) synthetic disv model. 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..f7e4dfa --- /dev/null +++ b/tests/test_nhflodata_contract.py @@ -0,0 +1,132 @@ +"""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 (``pyproject.toml`` sets ``NHFLODATA_LOCATION=`` through pytest-env, +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_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 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. + """ + 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(): + """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 index ba4876f..9373eef 100644 --- a/tests/test_polder.py +++ b/tests/test_polder.py @@ -1,23 +1,153 @@ -"""Tests for nhflotools.polder. +"""Tests for nhflotools.polder.drn_from_waterboard_data. -The polder DRN gains an optional ``exclude`` mask so carved lake cells can be handed over -to a dedicated stage boundary. The reach must be dropped whether its stage comes from real -HHNK peilgebied data or from the maaiveld fallback, which is why the exclusion nulls the -conductance after the celldata assignment rather than only editing the fallback mask. +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. + +The polder DRN also gains an optional ``exclude`` mask so carved lake cells can be handed +over to a dedicated stage boundary. The reach must be dropped whether its stage comes from +real HHNK peilgebied data or from the maaiveld fallback, which is why the exclusion nulls +the conductance after the celldata assignment rather than only editing the fallback mask. """ import geopandas as gpd import nlmod +import numpy as np import xarray as xr from shapely.geometry import box from nhflotools import polder +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 -# Cell ids of the shared 3x3 disv grid. +# Cell ids of the 3x3 disv_grid used by the exclude tests. CELL_REAL_STAGE = 4 # centre cell, covered by the mocked HHNK level area (real peilgebied stage) CELL_FALLBACK = 0 # a cell with no HHNK stage -> maaiveld fallback drain +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]) + + def _reach_cells(drn): """Return the set of icell2d values that carry a DRN reach.""" if drn is None: 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", +]