Aggregate lake pieces per cell in lake_from_gdf, like RIV/DRN celldata aggregation - #582
Aggregate lake pieces per cell in lake_from_gdf, like RIV/DRN celldata aggregation#582bdestombe wants to merge 1 commit into
Conversation
… strt A lake commonly intersects a grid cell in multiple polygon pieces (e.g. straight from nlmod.grid.gdf_to_grid). lake_from_gdf turned each row into its own VERTICAL connection over the full cell area with bedleak = 1/clake, multiply-counting the lake-aquifer exchange. Pieces of one lake within a cell now collapse to a single connection with clake = cell_area / sum(piece_area / clake), the same area-weighted aggregation nlmod.gwf.surface_water.aggregate applies to RIV and DRN celldata. Numeric per-lake settings are now compared with np.allclose(rtol=1e-8) so aggregated inputs carrying float-level noise (e.g. area-weighted stages) are not rejected by the single-value check. Both regression tests fail on unchanged dev: three connections with bedleak 0.1 instead of two with area-weighted bedleak, and an AssertionError on a 1e-12 strt difference.
OnnoEbbens
left a comment
OnnoEbbens
left a comment
There was a problem hiding this comment.
Thanks Bas, hadn't realised this was happening for lakes. I have just some small textual comments.
Another thing is that I had a bit of a hard time understanding the PR text. It only clicked after I looked into the code. It looks as if it is AI generated. Not necessarily something I oppose but maybe the AI can make it more concise next time :)
| the first piece per cell. | ||
| """ | ||
| if "geometry" not in lake_gdf.columns or lake_gdf.geometry.isna().any(): | ||
| raise ValueError( |
There was a problem hiding this comment.
I had a bit of trouble reading the error message. Is this a good alternative:
"found multiple elements of one lake in a single cell; provide polygon "
"geometries of the lake to aggregate these elements"
| lakeno : with the number of the lake | ||
| strt : with the starting head of the lake | ||
| clake : with the bed resistance of the lake | ||
| A lake may have multiple rows (polygon pieces) per cellid, e.g. straight |
There was a problem hiding this comment.
| A lake may have multiple rows (polygon pieces) per cellid, e.g. straight | |
| A single lake may have multiple elements (polygon pieces) in one cell, e.g. straight |
OnnoEbbens
left a comment
OnnoEbbens
left a comment
There was a problem hiding this comment.
Thanks Bas, hadn't realised this was happening for lakes. I have just some small textual comments.
Another thing is that I had a bit of a hard time understanding the PR text. It only clicked after I looked into the code. It looks as if it is AI generated. Not necessarily something I oppose but maybe the AI can make it more concise next time :)
There was a problem hiding this comment.
This is a good improvement. Existing models will change, as the bedleak (conductance) at the edge of lakes is now scaled by the part of the cell that is covered by the lake. This is a good thing though, and will maybe help with model convergence as well.
When Onno's comments are addressed, we can merge this.
3 participants
What
lake_from_gdfnow accepts per-piece input — multiple rows per cellid for one lake, as produced directly bynlmod.grid.gdf_to_grid— and combines the pieces into one lake-GWF connection per cell with an area-weighted bed resistance. It also tolerates float-level noise in numeric per-lake settings.Why
Two defects when the input is not pre-aggregated by the caller:
lake_from_gdfturns every row into its own connection withbedleak = 1/clake. A cell holding two pieces of one lake (6000 m² atclake=10and 2000 m² atclake=20in a 10000 m² cell) gets two full-cell-area connections totalling0.15 * cell_areaof conductance, where the pieces physically justify(6000/10 + 2000/20) / 10000 = 0.07 * cell_area._get_and_check_single_valuecomparesstrt(and outlet settings) with==. Area-weighted per-cell stages carry float-level noise (a 1e-12 difference), which raisedAssertionError: A single lake should have a single strtfor no physical reason.RIV and DRN already have exactly this aggregation step in nlmod:
nlmod.gwf.surface_water.aggregatecollapses per-piece celldata to one reach per cell with area-weighted parameters beforebuild_spd. The LAK path lacked the equivalent, forcing every caller to hand-roll it. (Encountered wiring the NHFLO Bergen ponds — NHFLO/models#127, NHFLO/tools#58 — where the workaround currently lives downstream.)How
_aggregate_connections_per_cell(lake_gdf, ds), invoked per lake when the cellid index has duplicates: combinedclake = cell_area / sum(piece_area / clake), so the VERTICAL connection conductancecell_area / clakeequals the summed piece conductance; other columns take the first piece's value (they are per-lake settings, still validated by_get_and_check_single_value). Raises a clearValueErrorwhen pieces come without geometry, since the areas are then unknown._get_and_check_single_valueusesnp.allclose(rtol=1e-8, atol=0)for numeric columns; exact comparison is kept for strings (boundnames, timeseries names,lakeout).nlakeconnconsequently counts cells, not pieces.Behaviour for already-aggregated input (one row per cell per lake) is unchanged.
Tests
Two regression tests in
test_013_surface_water.py, both shown to fail on unchangeddevbefore the fix:test_lake_from_gdf_aggregates_pieces_per_cell— fails on dev with 3 connections atbedleak=0.1/0.05instead of 2 connections at the area-weighted0.07/0.1; also pinsnlakeconn == 2.test_lake_from_gdf_accepts_floating_point_strt_noise— fails on dev with theAssertionErroron a1e-12strt difference.test_gdf_lake(existing) still passes;ruff check/ruff formatadd no findings relative to dev.