Bug fix writing clipped idf_svat.

Author: verkaik

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

@@ -2,7 +2,6 @@
 from typing import Optional
 
 import numpy as np
-import xarray as xr
 
 from imod.common.interfaces.iregridpackage import IRegridPackage
 from imod.common.utilities.dataclass_type import DataclassType
@@ -47,25 +46,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)

imod/msw/pkgbase.py

@@ -139,6 +139,30 @@ def _render(
         Collect to be written data in a DataFrame and call
         ``self.write_dataframe_fixed_width``
         """
+        if self.__class__.__name__ == "IdfMapping":
+            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
+
         data_dict: DataDictType = {"svat": svat.values.ravel()[index]}
 
         subunit = svat.coords["subunit"]

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