Issue #1772 Bug fix writing clipped idf_svat.

docs/api/changelog.rst

@@ -45,6 +45,8 @@ Fixed
   produced incorrect results when point water heads were below elevation levels
   for unstructured grids.
 - Support pandas 3.0.
+- :class:`imod.msw.IdfMapping` when model clipping is applied, the global
+  row/column indices are converted to local indices, as writen in idf_svat.inp.
 
 
 Changed

imod/msw/idf_mapping.py

@@ -1,16 +1,18 @@
 from dataclasses import asdict
-from typing import Optional
+from typing import Optional, TextIO
 
 import numpy as np
 import xarray as xr
 
 from imod.common.interfaces.iregridpackage import IRegridPackage
 from imod.common.utilities.dataclass_type import DataclassType
 from imod.common.utilities.regrid import _regrid_array
+from imod.mf6.dis import StructuredDiscretization
+from imod.mf6.mf6_wel_adapter import Mf6Wel
 from imod.msw.fixed_format import VariableMetaData
 from imod.msw.pkgbase import MetaSwapPackage
 from imod.msw.regrid.regrid_schemes import IdfMappingRegridMethod
-from imod.typing import GridDataArray
+from imod.typing import GridDataArray, IntArray
 from imod.util.regrid import RegridderWeightsCache
 from imod.util.spatial import spatial_reference
 
@@ -47,25 +49,6 @@ def __init__(self, area, nodata):
         self.dataset["area"] = area
         self.dataset["nodata"] = nodata
 
-        nrow = self.dataset.coords["y"].size
-        ncol = self.dataset.coords["x"].size
-
-        y_index = xr.DataArray(
-            np.arange(1, nrow + 1), coords={"y": self.dataset.coords["y"]}, dims=("y",)
-        )
-        x_index = xr.DataArray(
-            np.arange(1, ncol + 1), coords={"x": self.dataset.coords["x"]}, dims=("x",)
-        )
-        rows, columns = xr.broadcast(y_index, x_index)
-
-        self.dataset["rows"] = rows
-        self.dataset["columns"] = columns
-
-        y_grid, x_grid = xr.broadcast(self.dataset["y"], self.dataset["x"])
-
-        self.dataset["x_grid"] = x_grid
-        self.dataset["y_grid"] = y_grid
-
     def _get_output_settings(self):
         grid = self.dataset["area"]
         dx, xmin, _, dy, ymin, _ = spatial_reference(grid)
@@ -152,3 +135,37 @@ def regrid_like(
         )
 
         return type(self)(regridded_area, nodata)
+
+    def _render(
+        self,
+        file: TextIO,
+        index: IntArray,
+        svat: xr.DataArray,
+        mf6_dis: StructuredDiscretization,
+        mf6_well: Mf6Wel,
+    ) -> None:
+
+        nrow = self.dataset.coords["y"].size
+        ncol = self.dataset.coords["x"].size
+
+        y_index = xr.DataArray(
+            np.arange(1, nrow + 1),
+            coords={"y": self.dataset.coords["y"]},
+            dims=("y",),
+        )
+        x_index = xr.DataArray(
+            np.arange(1, ncol + 1),
+            coords={"x": self.dataset.coords["x"]},
+            dims=("x",),
+        )
+        rows, columns = xr.broadcast(y_index, x_index)
+
+        self.dataset["rows"] = rows
+        self.dataset["columns"] = columns
+
+        y_grid, x_grid = xr.broadcast(self.dataset["y"], self.dataset["x"])
+
+        self.dataset["x_grid"] = x_grid
+        self.dataset["y_grid"] = y_grid
+
+        super()._render(file, index, svat, mf6_dis, mf6_well)

imod/tests/test_msw/test_idf_mapping.py

@@ -16,8 +16,8 @@ def test_idf_mapping_regrid():
 
     idf_mapping_regridded = idf_mapping.regrid_like(new_grid, regrid_context)
 
-    assert len(idf_mapping_regridded.dataset["rows"]) == 5
-    assert len(idf_mapping_regridded.dataset["columns"]) == 5
+    assert len(idf_mapping_regridded.dataset["y"]) == 5
+    assert len(idf_mapping_regridded.dataset["x"]) == 5
     assert np.isnan(idf_mapping_regridded.dataset["nodata"].values[()])
 
 
@@ -27,8 +27,8 @@ def test_idf_mapping_clip():
     idf_mapping = IdfMapping(area, np.nan)
 
     idf_mapping_selected = idf_mapping.clip_box(
-        x_min=1.0, x_max=2.5, y_min=1.0, y_max=2.5
+        x_min=2.0, x_max=3.5, y_min=1.0, y_max=2.5
     )
 
-    expected_area = area.sel(x=slice(1.0, 2.5), y=slice(2.5, 1.0))
+    expected_area = area.sel(x=slice(2.0, 3.5), y=slice(2.5, 1.0))
     xr.testing.assert_allclose(idf_mapping_selected.dataset["area"], expected_area)

imod/tests/test_msw/test_output_control.py

@@ -114,6 +114,34 @@ def test_idf_oc_write_simple_model(fixed_format_parser):
     assert_almost_equal(results["x_grid"], np.array([2.0, 2.0, 2.0, 2.0, 4.0]))
 
 
+def test_idf_oc_write_simple_model_clip(fixed_format_parser):
+    area, index, svat = grid()
+    nodata = -9999.0
+    x_min, x_max, y_min, y_max = 1.5, 4.5, 1.5, 3.5
+
+    idf_output_control = IdfMapping(area, nodata)
+    idf_output_control_clip = idf_output_control.clip_box(
+        x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max
+    )
+
+    svat_clip = svat.sel(x=slice(x_min, x_max), y=slice(y_min, y_max))
+    index_clip = (svat_clip != 0).values.ravel()
+
+    with tempfile.TemporaryDirectory() as output_dir:
+        output_dir = Path(output_dir)
+        idf_output_control_clip.write(output_dir, index_clip, svat_clip, None, None)
+
+        results = fixed_format_parser(
+            output_dir / IdfMapping._file_name, IdfMapping._metadata_dict
+        )
+
+    assert_equal(results["svat"], np.array([2, 4, 5]))
+    assert_equal(results["rows"], np.array([2, 1, 1]))
+    assert_equal(results["columns"], np.array([1, 1, 3]))
+    assert_almost_equal(results["y_grid"], np.array([3.0, 2.0, 2.0]))
+    assert_almost_equal(results["x_grid"], np.array([2.0, 2.0, 4.0]))
+
+
 def test_idf_oc_mapping_order(fixed_format_parser):
     area, index, svat = grid()
     nodata = -9999.0