Issue #965 npf from imod5 (#1010)

Fixes #965 

# Description
Adds a function to import an npf package from an imod5 project file.
Adds tests.

# Checklist


- [x] Links to correct issue
- [ ] Update changelog, if changes affect users
- [x] PR title starts with ``Issue #nr``, e.g. ``Issue #737``
- [x] Unit tests were added
- [ ] **If feature added**: Added/extended example

---------

Co-authored-by: Joeri van Engelen <[email protected]>

imod/mf6/npf.py

@@ -1,12 +1,19 @@
 import warnings
+from copy import deepcopy
 from typing import Optional, Tuple
 
 import numpy as np
+import xarray as xr
 
 from imod.logging import init_log_decorator
 from imod.mf6.interfaces.iregridpackage import IRegridPackage
 from imod.mf6.package import Package
-from imod.mf6.utilities.regrid import RegridderType
+from imod.mf6.utilities.imod5_converter import fill_missing_layers
+from imod.mf6.utilities.regrid import (
+    RegridderType,
+    RegridderWeightsCache,
+    _regrid_package_data,
+)
 from imod.mf6.validation import PKG_DIMS_SCHEMA
 from imod.schemata import (
     AllValueSchema,
@@ -17,6 +24,7 @@
     IndexesSchema,
 )
 from imod.typing import GridDataArray
+from imod.typing.grid import zeros_like
 
 
 def _dataarray_to_bool(griddataarray: GridDataArray) -> bool:
@@ -461,3 +469,74 @@ def _validate(self, schemata, **kwargs):
 
     def get_regrid_methods(self) -> Optional[dict[str, Tuple[RegridderType, str]]]:
         return self._regrid_method
+
+    @classmethod
+    def from_imod5_data(
+        cls,
+        imod5_data: dict[str, dict[str, GridDataArray]],
+        target_grid: GridDataArray,
+        regridder_types: Optional[dict[str, tuple[RegridderType, str]]] = None,
+    ) -> "NodePropertyFlow":
+        """
+        Construct an npf-package from iMOD5 data, loaded with the
+        :func:`imod.formats.prj.open_projectfile_data` function.
+
+        .. note::
+
+            The method expects the iMOD5 model to be fully 3D, not quasi-3D.
+
+        Parameters
+        ----------
+        imod5_data: dict
+            Dictionary with iMOD5 data. This can be constructed from the
+            :func:`imod.formats.prj.open_projectfile_data` method.
+        target_grid: GridDataArray
+            The grid that should be used for the new package. Does not
+            need to be identical to one of the input grids.
+        regridder_types: dict, optional
+            Optional dictionary with regridder types for a specific variable.
+            Use this to override default regridding methods.
+
+        Returns
+        -------
+        Modflow 6 npf package.
+
+        """
+
+        data = {
+            "k": imod5_data["khv"]["kh"],
+        }
+        has_vertical_anisotropy = (
+            "kva" in imod5_data.keys()
+            and "vertical_anisotropy" in imod5_data["kva"].keys()
+        )
+        has_horizontal_anisotropy = "ani" in imod5_data.keys()
+
+        if has_vertical_anisotropy:
+            data["k33"] = data["k"] * imod5_data["kva"]["vertical_anisotropy"]
+        if has_horizontal_anisotropy:
+            if not np.all(np.isnan(imod5_data["ani"]["factor"].values)):
+                factor = imod5_data["ani"]["factor"]
+                factor = fill_missing_layers(factor, target_grid, 1)
+                data["k22"] = data["k"] * factor
+            if not np.all(np.isnan(imod5_data["ani"]["angle"].values)):
+                angle1 = imod5_data["ani"]["angle"]
+                angle1 = 90.0 - angle1
+                angle1 = xr.where(angle1 < 0, 360.0 + angle1, angle1)
+                angle1 = fill_missing_layers(angle1, target_grid, 0)
+                data["angle1"] = angle1
+
+        icelltype = zeros_like(target_grid, dtype=int)
+
+        regridder_settings = deepcopy(cls._regrid_method)
+        if regridder_types is not None:
+            regridder_settings.update(regridder_types)
+
+        regrid_context = RegridderWeightsCache(data["k"], target_grid)
+
+        new_package_data = _regrid_package_data(
+            data, target_grid, regridder_settings, regrid_context, {}
+        )
+        new_package_data["icelltype"] = icelltype
+
+        return NodePropertyFlow(**new_package_data)

imod/mf6/utilities/imod5_converter.py

@@ -1,3 +1,5 @@
+from typing import Union
+
 import numpy as np
 import xarray as xr
 
@@ -26,3 +28,17 @@ def convert_ibound_to_idomain(
     )
     # Fill the remaining nans at tops and bottoms with 0
     return idomain_float.fillna(0).astype(int)
+
+
+def fill_missing_layers(
+    source: xr.DataArray, full: xr.DataArray, fillvalue: Union[float | int]
+) -> xr.DataArray:
+    """
+    This function takes a source grid in which the layer dimension is
+    incomplete. It creates a result-grid which has the same layers as the "full"
+    grid, which is assumed to have all layers. The result has the values in the
+    source for the layers that are in the source. For the other layers, the
+    fillvalue is assigned.
+    """
+    layer = full.coords["layer"]
+    return source.reindex(layer=layer, fill_value=fillvalue)

imod/mf6/utilities/regrid.py

@@ -198,7 +198,7 @@ def _regrid_package_data(
     target_grid: GridDataArray,
     regridder_settings: dict[str, tuple[RegridderType, str]],
     regrid_context: RegridderWeightsCache,
-    new_package_data: Optional[dict[str, GridDataArray]] = {},
+    new_package_data: dict[str, GridDataArray] = {},
 ) -> dict[str, GridDataArray]:
     """
     Regrid package data. Loops over regridder settings to regrid variables one

imod/tests/conftest.py

@@ -1,5 +1,6 @@
 import pytest
 
+from .fixtures.backward_compatibility_fixture import imod5_dataset
 from .fixtures.flow_basic_fixture import (
     basic_dis,
     basic_dis__topsystem,

imod/tests/fixtures/backward_compatibility_fixture.py

@@ -0,0 +1,20 @@
+import pytest
+import xarray as xr
+
+import imod
+
+
+@pytest.fixture(scope="module")
+def imod5_dataset():
+    tmp_path = imod.util.temporary_directory()
+    data = imod.data.imod5_projectfile_data(tmp_path)
+    _load_imod5_data_in_memory(data[0])
+    return data[0]
+
+
+def _load_imod5_data_in_memory(imod5_data):
+    """For debugging purposes, load everything in memory"""
+    for pkg in imod5_data.values():
+        for vardata in pkg.values():
+            if isinstance(vardata, xr.DataArray):
+                vardata.load()

imod/tests/test_mf6/test_mf6_dis.py

@@ -8,14 +8,9 @@
 import imod
 from imod.mf6.write_context import WriteContext
 from imod.schemata import ValidationError
-
-
-def _load_imod5_data_in_memory(imod5_data):
-    """For debugging purposes, load everything in memory"""
-    for pkg in imod5_data.values():
-        for vardata in pkg.values():
-            if isinstance(vardata, xr.DataArray):
-                vardata.load()
+from imod.tests.fixtures.backward_compatibility_fixture import (
+    _load_imod5_data_in_memory,
+)
 
 
 @pytest.fixture(scope="function")

imod/tests/test_mf6/test_mf6_npf.py

@@ -1,6 +1,7 @@
 import pathlib
 import re
 import textwrap
+from copy import deepcopy
 
 import numpy as np
 import pytest
@@ -208,3 +209,87 @@ def test_configure_xt3d(tmp_path):
     assert "xt3d" not in rendered
     assert "rhs" not in rendered
     assert not npf.get_xt3d_option()
+
+
+@pytest.mark.usefixtures("imod5_dataset")
+def test_npf_from_imod5_isotropic(imod5_dataset, tmp_path):
+    data = deepcopy(imod5_dataset)
+    # throw out kva (=vertical anisotropy array) and ani (=horizontal anisotropy array)
+    data.pop("kva")
+    data.pop("ani")
+
+    target_grid = data["khv"]["kh"]
+    npf = imod.mf6.NodePropertyFlow.from_imod5_data(data, target_grid)
+
+    # Test array values are the same  for k ( disregarding the locations where k == np.nan)
+    k_nan_removed = xr.where(np.isnan(npf.dataset["k"]), 0, npf.dataset["k"])
+    np.testing.assert_allclose(k_nan_removed, data["khv"]["kh"].values)
+
+    rendered_npf = npf.render(tmp_path, "npf", None, None)
+    assert "k22" not in rendered_npf
+    assert "k33" not in rendered_npf
+    assert "angle1" not in rendered_npf
+    assert "angle2" not in rendered_npf
+    assert "angle3" not in rendered_npf
+
+
+@pytest.mark.usefixtures("imod5_dataset")
+def test_npf_from_imod5_horizontal_anisotropy(imod5_dataset, tmp_path):
+    data = deepcopy(imod5_dataset)
+    # throw out kva (=vertical anisotropy array)
+    data.pop("kva")
+
+    target_grid = data["khv"]["kh"]
+    data["ani"]["angle"].values[:, :, :] = 135.0
+    data["ani"]["factor"].values[:, :, :] = 0.1
+    npf = imod.mf6.NodePropertyFlow.from_imod5_data(data, target_grid)
+
+    # Test array values  for k22 and angle1
+    for layer in npf.dataset["k"].coords["layer"].values:
+        ds_layer = npf.dataset.sel({"layer": layer})
+
+        ds_layer = ds_layer.fillna(0.0)
+
+        if layer in data["ani"]["factor"].coords["layer"].values:
+            np.testing.assert_allclose(
+                ds_layer["k"].values * 0.1, ds_layer["k22"].values, atol=1e-10
+            )
+            assert np.all(ds_layer["angle1"].values == 315.0)
+        else:
+            assert np.all(ds_layer["k"].values == ds_layer["k22"].values)
+            assert np.all(ds_layer["angle1"].values == 0.0)
+
+    rendered_npf = npf.render(tmp_path, "npf", None, None)
+    assert "k22" in rendered_npf
+    assert "k33" not in rendered_npf
+    assert "angle1" in rendered_npf
+    assert "angle2" not in rendered_npf
+    assert "angle3" not in rendered_npf
+
+
+@pytest.mark.usefixtures("imod5_dataset")
+def test_npf_from_imod5_vertical_anisotropy(imod5_dataset, tmp_path):
+    data = deepcopy(imod5_dataset)
+    # throw out ani (=horizontal anisotropy array)
+    data.pop("ani")
+
+    data["kva"]["vertical_anisotropy"].values[:] = 0.1
+    target_grid = data["khv"]["kh"]
+
+    npf = imod.mf6.NodePropertyFlow.from_imod5_data(data, target_grid)
+
+    # Test array values  for k33
+    for layer in npf.dataset["k"].coords["layer"].values:
+        k_layer = npf.dataset["k"].sel({"layer": layer})
+        k33_layer = npf.dataset["k33"].sel({"layer": layer})
+
+        k_layer = xr.where(np.isnan(k_layer), 0.0, k_layer)
+        k33_layer = xr.where(np.isnan(k33_layer), 0.0, k33_layer)
+        np.testing.assert_allclose(k_layer.values * 0.1, k33_layer.values, atol=1e-10)
+
+    rendered_npf = npf.render(tmp_path, "npf", None, None)
+    assert "k22" not in rendered_npf
+    assert "k33" in rendered_npf
+    assert "angle1" not in rendered_npf
+    assert "angle2" not in rendered_npf
+    assert "angle3" not in rendered_npf