[bodemlagen_pwn_2024] Missing information

[bdestombe]

• Check Vincent's okay before uploading bodemlagen report. Agreed.
• Provide a mask for which the bodemlagen here provide better information than REGISII. We need to know where to use bodemlagen_pwn_2024 and where REGIS. With the current provided data we cannot make the distinction whether a layer is knowingly missing or that it falls outside the region of interest. Bergen and Koster mask marks where the layer is not present within the study area.
• Check bodemlagen_pwn_regis_koppeltabel.
• Merging of DS21 and DS21_0. No. In the Eem formatie are two clay layers present. S21_0 will be renamed to S20, is located in the top of the Eem formatie. Then 20m sand before the top of the S21. Koppeltabel aanpassen: EEK1 (of BEk1/2?)=S20 and EEK2=S21 -> Vincent doet de check. W21 opsplitsen in W20 en W21 en koppeltabel daarop aanpassen.
• Only seem to have looked at aquitards. What needs to be done for the aquifers? Vincent looks at botms af Bergen aquifers.
• Check basis aquitards for Bergen region is translated to top aquitards? Vincent's product only contains top of aquitards.
• bodemlagen_pwn_nhdz/v1.0.0/Bodemparams/Cwaarden_aquitards/DWAT_Boringen_Selectie.shp contains elevations and resistivities not included in the interpolation points.
• .

[bdestombe]

@vincentpost Let's discuss the above points in person in the near future.

[bdestombe]

Only use drilling data from areas in which the layer is present.

I am not sure if I agree with this approach. Drilling data in areas where the layer is not present, I expect the values there to be zero. It helps the thickness of the clay layers to go to zero near the edges.

Do we need to add extra points with thickness zero to enforce layer thickness to go to zero near its edges?

[bdestombe]

@vincentpost in NHFLO we only use vertex grids. Whenever you write code for NHFLO could you adapt it to vertex grids?

The following code is my translation of your code adapted to vertex grids:

"""Module containing functions to retrieve PWN bodemlagen."""

import logging
from pathlib import Path

import geopandas as gpd
import numpy as np
import pykrige.ok
import xarray as xr
from flopy.utils.gridintersect import GridIntersect
from nlmod.dims.grid import modelgrid_from_ds
from shapely.ops import unary_union

logger = logging.getLogger(__name__)


def get_pwn_aquitard_data(ds_regis: xr.Dataset, data_dir: Path, ix: GridIntersect, transition_length: float) -> dict:
    """
    Interpolate the thickness of the aquitard layers and the top of the aquitard layers using Kriging.

    The thickness of the aquitard layers is interpolated using the points in the file
    `dikte_aquitard/D{layer_name}/D{layer_name}_interpolation_points.geojson`.
    The top of the aquitard layers is interpolated using the points in the file
    `top_aquitard/T{layer_name}/T{layer_name}_interpolation_points.geojson`.
    The mask of the aquitard layers is defined in the file
    `dikte_aquitard/D{layer_name}/D{layer_name}_mask_combined.geojson`.

    Parameters
    ----------
    ds_regis : xr.Dataset
        The REGIS modellayer that contains the vertex grid.
    data_dir : Path
        The directory containing the data. Contains folders `dikte_aquitard` and `top_aquitard`.
    ix : flopy.utils.GridIntersect
        The index of the model grid.
    transition_length : float
        The length of the transition zone in meters.

    Returns
    -------
    dict
        A dictionary containing the interpolated values of the aquitard layers.
    """
    verbose = logger.level <= logging.DEBUG

    if ix is None:
        modelgrid = modelgrid_from_ds(ds_regis)
        ix = GridIntersect(modelgrid, method="vertex")

    ncell = ds_regis.sizes["icell2d"]
    layer_names = ["S11", "S12", "S13", "S21", "S22", "S31", "S32"]
    data = {}

    for name in layer_names:
        # Compute where the layer is _not_ present
        logger.info("Interpolating aquitard layer %s data and its transition zone", name)
        fp_mask = data_dir / "dikte_aquitard" / f"D{name}" / f"D{name}_mask_combined.geojson"
        gdf_mask = gpd.read_file(fp_mask)

        multipolygon = unary_union(gdf_mask.geometry)
        ids = ix.intersect(multipolygon, contains_centroid=False, min_area_fraction=0.5).cellids.astype(int)
        data[f"{name}_mask"] = np.zeros(ncell, dtype=bool)
        data[f"{name}_mask"][ids] = True

        # Compute where the layer transitions to REGIS
        multipolygon_transition = multipolygon.buffer(transition_length).difference(multipolygon)
        ids_trans = ix.intersect(
            multipolygon_transition, contains_centroid=False, min_area_fraction=0.5
        ).cellids.astype(int)
        data[f"{name}_transition"] = np.zeros(ncell, dtype=bool)
        data[f"{name}_transition"][ids_trans] = True

        # Interpolate thickness points using Kriging
        fp_pts = data_dir / "dikte_aquitard" / f"D{name}" / f"D{name}_interpolation_points.geojson"
        gdf_pts = gpd.read_file(fp_pts)
        ok = pykrige.ok.OrdinaryKriging(
            gdf_pts.geometry.x.values,
            gdf_pts.geometry.y.values,
            gdf_pts.value.values,
            variogram_model="linear",
            verbose=verbose,
            enable_plotting=False,
        )
        xq = ds_regis.x.values[~data[f"{name}_mask"]]
        yq = ds_regis.y.values[~data[f"{name}_mask"]]
        kriging_result = ok.execute("points", xq, yq)
        data[f"D{name}_value"] = np.zeros(ncell)
        data[f"D{name}_value"][~data[f"{name}_mask"]] = kriging_result[0]
        data[f"D{name}_value_unc"] = np.zeros(ncell)
        data[f"D{name}_value_unc"][~data[f"{name}_mask"]] = kriging_result[1]

        # Interpolate top aquitard points using Kriging
        fp_pts = data_dir / "top_aquitard" / f"T{name}" / f"T{name}_interpolation_points.geojson"
        gdf_pts = gpd.read_file(fp_pts)
        ok = pykrige.ok.OrdinaryKriging(
            gdf_pts.geometry.x.values,
            gdf_pts.geometry.y.values,
            gdf_pts.value.values,
            variogram_model="linear",
            verbose=verbose,
            enable_plotting=False,
        )
        kriging_result = ok.execute("points", xq, yq)
        data[f"T{name}_value"] = np.zeros(ncell)
        data[f"T{name}_value"][~data[f"{name}_mask"]] = kriging_result[0]
        data[f"T{name}_value_unc"] = np.zeros(ncell)
        data[f"T{name}_value_unc"][~data[f"{name}_mask"]] = kriging_result[1]

    return data