From c821ae0bc3abdf850e6633f66bb82f634fd89f57 Mon Sep 17 00:00:00 2001 From: Bas des Tombe Date: Wed, 22 Jul 2026 16:10:26 +0200 Subject: [PATCH] Fractional open-water coverage for the recharge mask and polder DRN - riv_from_oppervlakte_pwn stores ds['panden_coverage']: the panden-covered fraction of each cell from the intersected piece areas. - recharge_pond_mask(fractional=True) uses that fraction at panden cells (summed with the lake coverage and clipped at 1) instead of a hard 1.0; without the variable the previous 1.0-at-panden-cells behavior is kept. - drn_from_waterboard_data's exclude generalizes from a boolean drop mask to an open-water fraction: the conductance is scaled by the drainable land share 1 - exclude, dropping the reach where the share is zero. Boolean masks behave exactly as before. - Tests: known-answer panden coverage, combined lake+panden fraction with clipping, and land-share-scaled DRN conductance. --- src/nhflotools/geoconverter/utils.py | 2 +- src/nhflotools/lakes.py | 18 +++++++++++++----- src/nhflotools/panden.py | 13 ++++++++++++- src/nhflotools/polder.py | 24 +++++++++++++----------- tests/test_lakes.py | 20 ++++++++++++++++++++ tests/test_panden.py | 11 +++++++++++ tests/test_polder.py | 18 ++++++++++++++++++ 7 files changed, 88 insertions(+), 18 deletions(-) diff --git a/src/nhflotools/geoconverter/utils.py b/src/nhflotools/geoconverter/utils.py index 88a863a..139d234 100644 --- a/src/nhflotools/geoconverter/utils.py +++ b/src/nhflotools/geoconverter/utils.py @@ -119,7 +119,7 @@ def read_tabular( y_column: str | None, wkt_column: str | None, target_crs: str = "EPSG:28992", - **kwargs + **kwargs, ) -> gpd.GeoDataFrame: """Read tabular data with geometry information. diff --git a/src/nhflotools/lakes.py b/src/nhflotools/lakes.py index 33584d1..c0a38a1 100644 --- a/src/nhflotools/lakes.py +++ b/src/nhflotools/lakes.py @@ -314,10 +314,12 @@ def recharge_pond_mask(ds, panden_riv=None, *, fractional=False): fractional : bool, optional When False (default) return the boolean mask: True at carved lake cells and panden RIV cells. When True return a float fraction per cell instead: 1.0 at carved lake - cells (the whole cell is modeled as lake bed held by the RIV) and at panden RIV - cells, and the open-water coverage fraction at cells below the carve threshold — - precipitation on such a partial lake sliver feeds the lake (whose balance the - prescribed stage absorbs), not the aquifer, so the caller can scale recharge by + cells (the whole cell is modeled as lake bed held by the RIV), the panden-covered + fraction from ``ds['panden_coverage']`` (set by + :func:`nhflotools.panden.riv_from_oppervlakte_pwn`; 1.0 at panden RIV cells when + that variable is absent), and the open-water coverage fraction at cells below the + carve threshold — precipitation on such a partial lake or pand sliver feeds the + stage boundary, not the aquifer, so the caller can scale recharge by ``1 - fraction`` there instead of keeping the full land P-E. Returns @@ -332,4 +334,10 @@ def recharge_pond_mask(ds, panden_riv=None, *, fractional=False): mask.loc[{"icell2d": cells}] = True if not fractional: return mask - return xr.where(mask, 1.0, ds["lake_coverage"], keep_attrs=True) + if "panden_coverage" in ds: + coverage = ds["lake_coverage"] + ds["panden_coverage"] + elif panden_riv is not None: + coverage = xr.where(mask & ~ds["lake_cell"], 1.0, ds["lake_coverage"]) + else: + coverage = ds["lake_coverage"] + return xr.where(ds["lake_cell"], 1.0, coverage.clip(max=1.0), keep_attrs=True) diff --git a/src/nhflotools/panden.py b/src/nhflotools/panden.py index 9bcbdb0..11d999a 100644 --- a/src/nhflotools/panden.py +++ b/src/nhflotools/panden.py @@ -5,9 +5,11 @@ import geopandas as gpd import nlmod import numpy as np +import xarray as xr logger = logging.getLogger(__name__) + def get_oppervlakte_pwn_shapes(data_path_panden): """Get oppervlakte shapes for PWN area. @@ -59,7 +61,9 @@ def riv_from_oppervlakte_pwn(ds, gwf, data_path_panden): Returns ------- flopy.mf6.ModflowGwfriv or None - RIV package, or None when no data intersects the model extent. + RIV package, or None when no data intersects the model extent. As a side effect + ``ds['panden_coverage']`` is set: the panden-covered fraction of each cell, used by + :func:`nhflotools.lakes.recharge_pond_mask` for the fractional recharge mask. """ panden_shp = get_oppervlakte_pwn_shapes(data_path_panden=data_path_panden) rivdata = nlmod.grid.gdf_to_grid(panden_shp, ds) @@ -67,6 +71,13 @@ def riv_from_oppervlakte_pwn(ds, gwf, data_path_panden): logger.warning("No RIV data found within the provided extent.") return None + piece_area = rivdata.groupby("cellid")["area"].sum() + coverage = xr.zeros_like(ds["area"], dtype=float) + coverage.loc[{"icell2d": piece_area.index.to_numpy()}] = ( + piece_area.to_numpy() / ds["area"].sel(icell2d=piece_area.index.to_numpy()).to_numpy() + ).clip(max=1.0) + ds["panden_coverage"] = coverage + rivdata["cond"] = rivdata["area"] / rivdata["c"] agg = nlmod.grid.aggregate_vector_per_cell( rivdata, diff --git a/src/nhflotools/polder.py b/src/nhflotools/polder.py index e163305..dc12884 100644 --- a/src/nhflotools/polder.py +++ b/src/nhflotools/polder.py @@ -34,12 +34,13 @@ def drn_from_waterboard_data(ds, gwf, wb="Hollands Noorderkwartier", cbot=1.0, e Bottom resistance of the drains [days], by default 1.0. The per-cell conductance is ``cell_area / cbot``. exclude : xarray.DataArray or None, optional - Boolean mask over ``icell2d`` marking cells at which no DRN reach should be - emitted. Both the conductance and the elevation are nulled at these cells, so a - reach is dropped whether its stage comes from HHNK peilgebied data or from the - maaiveld fallback (``nlmod.gwf.drn`` masks on ``cond > 0``). The default None - applies no exclusion and is fully backward-compatible. Used to hand carved lake - cells over to a dedicated stage boundary. + Boolean mask or open-water fraction (0-1) over ``icell2d``. The conductance is + scaled by the drainable land share ``1 - exclude`` (a boolean ``True`` scales by + 0), so a reach is dropped at fully open-water cells — whether its stage comes from + HHNK peilgebied data or from the maaiveld fallback (``nlmod.gwf.drn`` masks on + ``cond > 0``) — and only drains the land share elsewhere. The default None applies + no exclusion and is fully backward-compatible. Used to hand carved lake cells over + to a dedicated stage boundary and to keep open water out of the drained area. Returns ------- @@ -91,10 +92,11 @@ def drn_from_waterboard_data(ds, gwf, wb="Hollands Noorderkwartier", cbot=1.0, e ds["drn_cond"] = drn_cond if exclude is not None: - # Null both fields at excluded cells. nlmod.gwf.drn builds a reach wherever - # cond > 0, so nulling drn_cond drops the reach whether its stage came from the - # celldata assignment or the maaiveld fallback above. - ds["drn_cond"] = ds["drn_cond"].where(~exclude) - ds["drn_elev"] = ds["drn_elev"].where(~exclude) + # Scale the conductance by the drainable land share; nlmod.gwf.drn builds a reach + # wherever cond > 0, so a zero land share drops the reach whether its stage came + # from the celldata assignment or the maaiveld fallback above. + land = (1.0 - exclude.astype(float)).clip(min=0.0) + ds["drn_cond"] *= land + ds["drn_elev"] = ds["drn_elev"].where(land > 0) return nlmod.gwf.drn(ds, gwf, elev="drn_elev", cond="drn_cond") diff --git a/tests/test_lakes.py b/tests/test_lakes.py index 996aa8f..4a339cb 100644 --- a/tests/test_lakes.py +++ b/tests/test_lakes.py @@ -384,3 +384,23 @@ def test_lak_gdf_from_lakes_pwn_returns_none_when_no_lake_cells(disv_grid): sliver = _lake_gdf([_piece(CELL_40, box(0, 400, 20, 420), 2.0, 1.0)]) ds_carved, _ = lakes.carve_lake_cells(ds, sliver, min_area_fraction=0.5) assert lakes.lak_gdf_from_lakes_pwn(ds_carved, sliver, min_area_fraction=0.5) is None + + +def test_pond_mask_fractional_uses_panden_coverage(disv_grid): + """With ``ds['panden_coverage']`` present the panden share is fractional and sums with the lake share.""" + ds, _gwf, geoms = disv_grid() + rows = [ + _piece(CELL_40, box(0, 400, 80, 600), 2.0, 1.0), # 40% lake -> not carved + _piece(CELL_FULL, geoms[CELL_FULL], 2.0, 1.0), # carved + ] + ds_carved, _ = lakes.carve_lake_cells(ds, _lake_gdf(rows), min_area_fraction=0.5) + panden_coverage = xr.zeros_like(ds_carved["lake_coverage"]) + panden_coverage.loc[{"icell2d": PANDEN_A}] = 0.3 + panden_coverage.loc[{"icell2d": CELL_40}] = 0.9 + ds_carved["panden_coverage"] = panden_coverage + + frac = lakes.recharge_pond_mask(ds_carved, None, fractional=True) + + assert frac.sel(icell2d=PANDEN_A).item() == pytest.approx(0.3) # fraction, not 1.0 + assert frac.sel(icell2d=CELL_40).item() == pytest.approx(1.0) # 0.4 lake + 0.9 panden, clipped + assert frac.sel(icell2d=CELL_FULL).item() == pytest.approx(1.0) # carved stays fully covered diff --git a/tests/test_panden.py b/tests/test_panden.py index 0792e16..3be4e5d 100644 --- a/tests/test_panden.py +++ b/tests/test_panden.py @@ -182,3 +182,14 @@ def test_riv_ssm_registration_is_transport_gated_and_idempotent(tmp_path, transp assert ds.attrs["ssm_sources"] == [riv.package_name] else: assert "ssm_sources" not in ds.attrs + + +def test_riv_stores_panden_coverage(tmp_path, panden_ds_gwf): + """The panden-covered fraction of each cell is stored for the fractional recharge mask.""" + ds, gwf = panden_ds_gwf + # 50 m x 50 m pand wholly inside cell 0 of the 2x2 grid of 100 m x 100 m cells. + path = _write_panden_shp(tmp_path, [box(10, 110, 60, 160)], ["ICAS-noord"]) + + riv_from_oppervlakte_pwn(ds, gwf, data_path_panden=path) + + np.testing.assert_allclose(ds["panden_coverage"].values, [2500.0 / 10000.0, 0.0, 0.0, 0.0]) diff --git a/tests/test_polder.py b/tests/test_polder.py index 9373eef..03e1bf3 100644 --- a/tests/test_polder.py +++ b/tests/test_polder.py @@ -17,6 +17,7 @@ import geopandas as gpd import nlmod import numpy as np +import pytest import xarray as xr from shapely.geometry import box @@ -187,3 +188,20 @@ def test_drn_fallback_excludes_lake_cells(disv_grid, monkeypatch): assert CELL_REAL_STAGE not in excl_cells assert CELL_FALLBACK in excl_cells # a fallback (maaiveld) cell is unaffected assert excl_cells == base_cells - {CELL_REAL_STAGE} + + +def test_drn_fractional_exclude_scales_land_share(disv_grid, monkeypatch): + """A fractional exclude scales the conductance by ``1 - fraction`` and drops fully covered cells.""" + ds, gwf, _geoms = disv_grid() + frac = xr.zeros_like(ds["top"]).astype(float) + frac.loc[{"icell2d": CELL_REAL_STAGE}] = 0.4 # 40% open water at the real-stage cell + frac.loc[{"icell2d": CELL_FALLBACK}] = 1.0 # fully open water -> reach dropped + + drn = _run_drn(monkeypatch, ds, gwf, exclude=frac) + conds = {rec["cellid"][-1]: float(rec["cond"]) for rec in drn.stress_period_data.data[0]} + + # cbot=1.0, so the unscaled conductance equals the cell area. + assert conds[CELL_REAL_STAGE] == pytest.approx(0.6 * float(ds["area"].sel(icell2d=CELL_REAL_STAGE))) + assert CELL_FALLBACK not in conds + untouched = next(c for c in conds if c not in {CELL_REAL_STAGE, CELL_FALLBACK}) + assert conds[untouched] == pytest.approx(float(ds["area"].sel(icell2d=untouched)))