Process KD/C masks and improve overlap handling

kh: Rework kh processing to compute C and KD layers from NHDZ and Bergen sources. Added mask-based clipping of KD inputs and assembled KD formulas into summed parts before filling boundaries; write overlap and fill GeoJSONs for C layers. Adjusted Bergen constants and name mapping and changed several fill/overlap priorities (use 'last' and 'sum' where appropriate).

params_helper_functions: Add _compute_pairwise_overlaps and clip_polygons_to_mask helpers, strengthen geometry buffering/validation, and extend fill_boundary_with_polygons to return (result, overlaps_gdf, fill_gdf), detect/report overlaps more robustly, and support an overlap_priority 'sum' mode that sums values in overlapping regions. Overall improves handling of overlapping geometries and mask-based clipping for kh generation.

Author: bdestombe

src/nhflodata/data/mockup/bodemlagen_pwn_2024/v2.0.0/kh/kh.py

@@ -1,10 +1,8 @@
 import logging
 
 import geopandas as gpd
-import matplotlib.pyplot as plt
-import numpy as np
-from params_helper_functions import fill_boundary_with_polygons
-from shapely import MultiPolygon, Polygon, make_valid
+import pandas as pd
+from params_helper_functions import clip_polygons_to_mask, fill_boundary_with_polygons
 
 from nhflodata.get_paths import get_abs_data_path
 
@@ -16,17 +14,34 @@
 data_path_nhdz = get_abs_data_path("bodemlagen_pwn_nhdz", "1.0.0")
 
 layer_names = ["W11", "S11", "W12", "S12", "W13", "S13", "W21", "S21", "W22", "S22", "W31", "S31", "W32", "S32"]
+c_nhdz_fill_values = {  # From bodemlagen_pwn_nhdz/v1.0.0/Ontleding_model_v2.xlsx
+    "S11": 1.0,
+    "S12": 10.0,
+    "S13": 10.0,
+    "S21": 10.0,
+    "S22": 10.0,
+    "S31": 10.0,
+    "S32": 10.0,
+}
+
 kh_bergen_fill_values = {
-    "W11": 8.0,
-    "S11": np.nan,
+    "W11": 7.0,  # Triwaco Excel uses 7.0, but prior nhflo implementations used 8.0.
     "W12": 7.0,
-    "S12": np.nan,
     "W13": 12.0,
-    "S13": np.nan,
     "W21": 15.0,
-    "S21": np.nan,
 }
-
+c_bergen_fill_values = {
+    "S11": 1.0,
+    "S12": 100.0,
+    "S13": 100.0,
+    "S21": 1.0,
+}
+bergen_name_mapping = {
+    "11": "1A",
+    "12": "1B",
+    "13": "1C",
+    "21": "2",
+}
 
 # Get boundaries that we need to fill up with the data from Bergen and NHDZ
 gdf_bergen = gpd.read_file(data_path / "boundaries" / "triwaco_model_bergen.geojson")
@@ -38,31 +53,127 @@
     if name.startswith("W"):
         fp = data_path_nhdz / "Bodemparams" / "Kwaarden_aquifers" / f"K{name}.shp"
         k = gpd.read_file(fp).set_crs("EPSG:28992")
-        kh_nhdz[name], _ = fill_boundary_with_polygons(
+        kh_nhdz[name], _, _ = fill_boundary_with_polygons(
             boundary_gdf=gdf_nhdz,
             source_gdf=k,
             value_col="VALUE",
-            fill_value=kh_bergen_fill_values[name],  # Not used as fill method is 'mean'
+            # fill_value=kh_bergen_fill_values[name],  # Not used as fill method is 'mean'
             fill_method="mean",
-            overlap_priority="smallest",
+            overlap_priority="last",
             override_gdf=None,
         )
 
-# Bergen
-kh_bergen = {}
+kd_nhdz = {}  # horizontal transmissivity of the aquitards from nhdz
+folder = data_path_nhdz / "Bodemparams" / "KDwaarden_aquitards"
+s12kd = gpd.read_file(folder / "s12kd.shp").set_crs("EPSG:28992")
+s13kd = gpd.read_file(folder / "s13kd.shp").set_crs("EPSG:28992")
+s21kd = gpd.read_file(folder / "s21kd.shp").set_crs("EPSG:28992")
+s22kd = gpd.read_file(folder / "s22kd.shp").set_crs("EPSG:28992")
+folder = data_path_bergen / "Bodemparams" / "Masker_kdwaarden_aquitards"
+ms12kd = gpd.read_file(folder / "masker_aquitard12_kd.shp").set_crs("EPSG:28992")
+ms13kd = gpd.read_file(folder / "masker_aquitard13_kd.shp").set_crs("EPSG:28992")
+ms21kd = gpd.read_file(folder / "masker_aquitard21_kd.shp").set_crs("EPSG:28992")
+ms22kd = gpd.read_file(folder / "masker_aquitard22_kd.shp").set_crs("EPSG:28992")
+
+# Each formula is a list of (source_gdf, mask_gdf, mask_value, coefficient) terms
+#   kd12s = s12kd*(ms12kd==1) + 0.5*s12kd*(ms12kd==2) + 3*s12kd*(ms12kd==3)
+#   kd13s = (s13kd*(ms13kd==1) + 0.5*s12kd*(ms12kd==2)) * 1.04
+#   kd21s = (s21kd*(ms21kd==1) + s13kd*(ms13kd==2)) * 0.6
+#   kd22s = s22kd*(ms22kd==1) + s21kd*(ms21kd==2)
+#   kd31s = s22kd*(ms22kd==2)
+kd_formulas = {
+    "S12": [
+        (s12kd, ms12kd, 1, 1.0),
+        (s12kd, ms12kd, 2, 0.5),
+        (s12kd, ms12kd, 3, 3.0),
+    ],
+    "S13": [
+        (s13kd, ms13kd, 1, 1.04),
+        (s12kd, ms12kd, 2, 0.5 * 1.04),
+    ],
+    "S21": [
+        (s21kd, ms21kd, 1, 0.6),
+        (s13kd, ms13kd, 2, 0.6),
+    ],
+    "S22": [
+        (s22kd, ms22kd, 1, 1.0),
+        (s21kd, ms21kd, 2, 1.0),
+    ],
+    "S31": [
+        (s22kd, ms22kd, 2, 1.0),
+    ],
+}
+
+for name, terms in kd_formulas.items():
+    parts = [clip_polygons_to_mask(src, mask, mask_value=val, coefficient=coeff) for src, mask, val, coeff in terms]
+    source = gpd.GeoDataFrame(pd.concat(parts, ignore_index=True), crs="EPSG:28992")
+    kd_nhdz[name], _, _ = fill_boundary_with_polygons(
+        boundary_gdf=gdf_nhdz,
+        source_gdf=source,
+        value_col="VALUE",
+        fill_method="fill_value",
+        fill_value=0.01,  # From bodemlagen_pwn_nhdz/v1.0.0/Ontleding_model_v2.xlsx 
+        overlap_priority="sum",
+    )
+
+c_nhdz = {}
 for name in layer_names:
-    if name.startswith("W"):
-        kh_bergen[name], _ = fill_boundary_with_polygons(
-            boundary_gdf=gdf_bergen,
-            source_gdf=None,  # Are constants for Bergen, so use fill value only
-            value_col="VALUE",  # Not used as source_gdf is None
-            fill_value=kh_bergen_fill_values[name],
+    if name.startswith("S"):
+        fp = data_path_nhdz / "Bodemparams" / "Cwaarden_aquitards" / f"C{name}.shp"
+        c = gpd.read_file(fp).set_crs("EPSG:28992")
+        if name == "S22":
+            c["VALUE"] *= 1.02  # From bodemlagen_pwn_nhdz/v1.0.0/Ontleding_model_v2.xlsx
+
+        c_nhdz[name], overlaps_gdf, fill_gdf = fill_boundary_with_polygons(
+            boundary_gdf=gdf_nhdz,
+            source_gdf=c,
+            value_col="VALUE",
+            fill_value=c_nhdz_fill_values[name],
             fill_method="fill_value",
-            overlap_priority="smallest",  # Not used as source_gdf is None
+            overlap_priority="last",
             override_gdf=None,
         )
+        overlaps_gdf.to_file(data_path / "kh" / f"C{name}_NHDZ_overlaps.geojson", driver="GeoJSON")
+        fill_gdf.to_file(data_path / "kh" / f"C{name}_NHDZ_fill.geojson", driver="GeoJSON")
+
+# Bergen
+kh_bergen = {}
+for name in bergen_name_mapping:
+    kh_bergen[f"W{name}"], _, _ = fill_boundary_with_polygons(
+        boundary_gdf=gdf_bergen,
+        source_gdf=None,  # Are constants for Bergen, so use fill value only
+        value_col="VALUE",  # Not used as source_gdf is None
+        fill_value=kh_bergen_fill_values[f"W{name}"],
+        fill_method="fill_value",
+        overlap_priority="last",  # Not used as source_gdf is None
+        override_gdf=None,
+    )
+
+c_bergen = {}
+for name, bergen_name in bergen_name_mapping.items():
+    fp = data_path_bergen / "Bodemparams" / f"C{bergen_name}.shp"
+    c = gpd.read_file(fp).set_crs("EPSG:28992")
+    c_bergen[f"S{name}"], overlaps_gdf, fill_gdf = fill_boundary_with_polygons(
+        boundary_gdf=gdf_bergen,
+        source_gdf=c,
+        value_col="VALUE",  # Not used as source_gdf is None
+        fill_value=c_bergen_fill_values[f"S{name}"],
+        fill_method="fill_value",
+        overlap_priority="last",  # Not used as source_gdf is None
+        override_gdf=None,
+    )
+    overlaps_gdf.to_file(data_path / "kh" / f"C{bergen_name}_Bergen_overlaps.geojson", driver="GeoJSON")
+    fill_gdf.to_file(data_path / "kh" / f"C{bergen_name}_Bergen_fill.geojson", driver="GeoJSON")
 
 """
+# NHDZ Bergen
+S21 S2
+S13 S1C
+S12 S1B
+S11 S1A
+
+
+
     kh[0] = 8.0
     kh[1] = thickness[1] / clist[0] / f_anisotropy
     kh[2] = 7.0