Add drain_from_df helper

[bdestombe]

Summary

• add nlmod.gwf.drain.drain_from_df for DRN packages from vector/tabular drain input
• compute clipped line/polygon conductance and preserve separate drain activation thresholds
• expose deterministic MVR provider mapping and helper perioddata construction
• document overlap with nlmod.gwf.surface_water.gdf_to_seasonal_pkg
• extend nlmod.layers.get_layer_of_z with sparse/vector cellid support and nearest-active IDOMAIN-aware lookup
• use vectorized get_layer_of_z in drain_from_df for batched explicit 3D DRN remapping

Relationship to surface_water

There is deliberate overlap with nlmod.gwf.surface_water.gdf_to_seasonal_pkg for polygon-to-DRN conversion. Use gdf_to_seasonal_pkg for surface-water polygons with winter/summer stages and seasonal conductance timeseries. Use drain_from_df for fixed drain features (line, polygon, point, or direct cellid), and when deterministic MVR provider IDs are needed for routing drain outflow.

drain_from_df reuses surface_water.build_spd for vector and 2D-cellid layer placement, so those inputs are checked against get_idomain(ds): fully inactive columns are omitted and drains are placed only in active cells (idomain > 0). FloPy/MODFLOW 6 receive explicit 3D cellids and do not relocate DRN boundaries. Therefore explicit 3D cellids targeting inactive (idomain == 0) or vertical pass-through (idomain < 0) cells are remapped by calling nlmod.layers.get_layer_of_z(..., nearest_active=True) for all affected direct rows in one vectorized batch.

get_layer_of_z now remains backward compatible for dense full-grid calls, and additionally supports scalar/vector sparse cellid + z calls. This makes the same IDOMAIN-aware layer lookup reusable for other boundary packages such as wells.

Validation

• uv run -q pytest tests\test_009_layers.py::test_get_layer_of_z tests\test_009_layers.py::test_get_layer_of_z_above_model tests\test_009_layers.py::test_get_layer_of_z_below_model tests\test_009_layers.py::test_get_layer_of_z_nearest_active_cellids tests\test_009_layers.py::test_get_layer_of_z_nearest_active_vector_cellids tests\test_027_drain.py tests\test_013_surface_water.py::test_gdf_to_seasonal_pkg -v
• uv run -q ruff format nlmod\dims\layers.py nlmod\gwf\drain.py tests\test_009_layers.py tests\test_027_drain.py
• uv run -q ruff check --fix nlmod\gwf\drain.py tests\test_027_drain.py
• uv run -q ruff check --fix --extend-ignore D103,D205,D400,B028,B905,E501,PT018,PT028 nlmod\dims\layers.py tests\test_009_layers.py
• uv run -q python -m py_compile nlmod\dims\layers.py nlmod\gwf\drain.py tests\test_009_layers.py tests\test_027_drain.py

[bdestombe]

Review notes from exercising this branch against the NHFLO 09pwnmodel2 local-drains use case. Three items, verified against the branch head.

1. An inactive explicit 3D cellid with a NaN elevation aborts the whole call

_build_spd_with_provider_mapping sends inactive 3D cellids into get_layer_of_z(z=stage, nearest_active=True, ...) to remap them (nlmod/gwf/drain.py:386). When that stage/elevation is NaN, get_layer_of_z raises z must be finite to determine the layer, which is re-raised as Cannot remap DRN cellids; ... and kills the entire drain_from_df call. But an active 3D cellid with a NaN elevation is quietly omitted (see test_drain_from_df_omits_3d_cellids_without_elevation). So the same NaN row is omitted or fatal depending only on whether its cell happens to be active.

Reproduction on the current branch (fails):

def test_inactive_3d_cellid_with_nan_elevation_is_omitted_not_aborted():
    ds = test_010_wells.get_model_ds()
    ds["active_domain"] = ds["botm"].notnull()
    ds["active_domain"].data[0, 0, 0] = False           # make (0,0,0) inactive
    assert nlmod.dims.layers.get_idomain(ds).data[:, 0, 0].tolist() == [0, 1, 1]
    _, gwf = test_010_wells.get_sim_and_gwf(ds)
    drains = pd.DataFrame({"cellid": [(0, 0, 0)], "elevation": [np.nan], "cond": [5.0]})
    drn, mapping = nlmod.gwf.drain.drain_from_df(
        drains, gwf, ds, silent=True, return_provider_mapping=True)
    assert drn is None and mapping.empty

ValueError: Cannot remap DRN cellids; z must be finite to determine the layer

Suggestion: drop NaN-elevation rows before the remap batch (or skip them in it), so they are omitted consistently regardless of idomain.

2. get_layer_of_z silently ignores nearest_active / preferred_layer / idomain in dense (cellid=None) calls

The new IDOMAIN-aware logic lives entirely inside the if cellid is not None: block. The dense full-grid path (cellid=None) falls through to the original z-only lookup and never consults idomain, nearest_active, or preferred_layer. A caller writing get_layer_of_z(ds, z, nearest_active=True) (no cellid) gets a silent plain lookup with no IDOMAIN awareness. Suggest either honoring these in the dense path too, or raising when they are combined with cellid=None so the no-op can't pass unnoticed.

3. New nlmod.gwf.drain module is missing from the API docs

docs/modules.rst lists nlmod.gwf.surface_water, ...horizontal_flow_barrier, ...wells, etc., but has no .. automodule:: nlmod.gwf.drain, so drain_from_df won't appear in the rendered API reference. Add an entry alongside the other nlmod.gwf.* modules.

[bdestombe]

Reviewed this PR while assessing how the NHFLO 09pwnmodel2 model could adopt drain_from_df (its manual local-drains block maps onto it 1:1). Overall the design is sound and the adoption fit is good — a few defects and generalization gaps below. I ran the new code in a scratch harness against installed nlmod 0.11.3dev / flopy 3.11.0.dev0 and exercised the edge cases empirically.

Defects

1. NaN-elevation 3D-cellid rows targeting inactive cells crash the whole call (nlmod/gwf/drain.py, remap path ~L389). A drain row with an explicit 3D cellid and NaN elevation is (correctly) omitted when the cell is active — matching test_drain_from_df_omits_3d_cellids_without_elevation and build_spd's NaN-stage drop — but the same row targeting an idomain<=0 cell reaches the remap and raises ValueError("... z must be finite to determine the layer"), aborting the entire drain_from_df call and discarding all other valid rows. Fix: exclude pd.isna(stage) rows from remap_positions (e.g. remap_positions = positions[~active & ~celldata["stage"].isna().to_numpy()[positions]]) so NaN-elevation rows are uniformly omitted regardless of the target cell's idomain.

2. nearest_active can silently return layer 0 on a NaN botm (nlmod/dims/layers.py, get_layer_of_z nearest_active branch ~L1704). A NaN botm produces NaN distances (and poisons the derived top of the layer below via layer_tops = vstack((top0, layer_botms[:-1]))); NaN comparisons make the candidate mask all-False for that column, and np.argmax(all-False) returns 0 — so the function silently returns layer 0 even when z is unambiguously inside a deeper active layer. Fix: after distances[~active] = np.inf, also set non-finite distances to inf (distances[~np.isfinite(distances)] = np.inf) and raise if a column's minimum distance is inf, mirroring the existing "no active layers" error.

3. Sparse-cellid path diverges from the dense path for NaN z, and above_model/below_model are silently inert under nearest_active (nlmod/dims/layers.py ~L1666). The dense full-grid path returns below_model (-999) for NaN z, but the new sparse path raises ValueError("z must be finite..."); and in nearest_active mode the above_model/below_model parameters have no effect (out-of-range z is clipped to the nearest active layer) with no note of either divergence. Fix: document both behaviors, or align the sparse non-nearest path with the dense semantics.

4. A GeoDataFrame with a pre-existing cellid column skips grid clipping and computes conductance from unclipped geometry (nlmod/gwf/drain.py, _geodataframe_celldata ~L326). Intersection with the grid only happens when cellid is absent; a GeoDataFrame carrying a cellid column with unclipped (whole-feature) geometries gets cond = full_length_or_area * rate per row, silently double-counting conductance for every feature spanning multiple cells — contradicting the docstring's unconditional "Line and polygon geometries are intersected with the model grid." Fix: always re-intersect GeoDataFrame input, or document that a supplied cellid implies pre-clipped per-cell geometries and warn when a duplicated source index suggests otherwise.

5. Per-row build_spd recomputes full-grid idomain for every 2D record (nlmod/gwf/drain.py, _build_spd_with_provider_mapping ~L427). Each 2D-cellid row is wrapped in a one-row DataFrame and passed to surface_water.build_spd, which re-runs get_idomain(ds) (full-grid thickness) plus ds.top/.botm/.kh extraction every call — O(n_rows × grid_size). Measured ~45x slowdown at 200 rows on a toy grid, scaling with grid size. build_spd already accepts multi-row celldata and returns records in input order, so the 2D block can be built in one call with provider ids reconstructed from a per-row record count. This is also a precondition for later migrating polder.py (which covers the whole model domain).

6. mvr_provider_id validity is tied to the persistent period-0 DRN block, but the docstring implies unconditional validity (nlmod/gwf/drain.py ~L437). MF6 defines the MVR provider id as the boundary number as listed per stress period. The stored mapping is correct only while the period-0 block persists to all periods; if a user later sets drn.stress_period_data for period > 0 (e.g. seasonal conductance), the ids silently become stale for those periods. Fix: state in the drain_from_df / mvr_perioddata_from_provider_mapping docstrings that the mapping holds for all periods only while the period-0 DRN block persists.

Generalization gaps (to fully replace the manual NHFLO blocks)

• Expose pkg + separate stage/rbot + aux pass-through so drain_from_df can also emit RIV/GHB. build_spd already supports pkg in {RIV, DRN, GHB}, a distinct rbot column, and an aux column; surfacing them would let this replace nhflotools.panden.riv_from_oppervlakte_pwn (which needs RIV with rbot = stage - 2.0, a CONCENTRATION aux, and ssm_sources registration). Currently _base_celldata forces rbot = stage = elevation and the function is DRN-only.
• Add a layer_method="first_active" option (backed by the existing get_first_active_layer_from_idomain). nhflotools.polder.drn_from_waterboard_data places every polder/maaiveld drain in the first active layer; migrating it to drain_from_df with the current lay_of_rbot would silently move drains whose stage lies below the top active layer's bottom into deeper layers. A first_active option makes that migration behavior-preserving.

Adoption notes (for the NHFLO side — no change needed here)

Verified non-blockers for adopting drain_from_df in 09pwnmodel2: dropping ds["thickness"] before the call is fine (build_spd/get_idomain recompute from top/botm); vector-input layer placement is unchanged (both routes use build_spd lay_of_rbot — the new get_layer_of_z(nearest_active=True) path only fires for explicit 3D cellids). The intended behavior change to accept is per-feature activation thresholds instead of the current per-cell collapse to min elevation + summed conductance (a physics improvement — the PR's own test shows drainage decreases where features of different elevation share a cell). Combined with mover_destinations + mvr_perioddata_from_provider_mapping, this is exactly what would wire the Bergen pumping-station drains into a basin lake via MVR.