From f81d034ce2d0b4e52bb1c662d698f76dadadc05c Mon Sep 17 00:00:00 2001 From: verkaik Date: Wed, 11 Mar 2026 16:27:30 +0100 Subject: [PATCH 1/2] Bug fix writing clipped idf_svat. --- docs/api/changelog.rst | 2 ++ imod/msw/idf_mapping.py | 20 ---------------- imod/msw/pkgbase.py | 24 +++++++++++++++++++ imod/tests/test_msw/test_idf_mapping.py | 8 +++---- imod/tests/test_msw/test_output_control.py | 28 ++++++++++++++++++++++ 5 files changed, 58 insertions(+), 24 deletions(-) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 2fc2228bf..6a8dd8256 100644 --- a/docs/api/changelog.rst +++ b/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 diff --git a/imod/msw/idf_mapping.py b/imod/msw/idf_mapping.py index a156ae0da..933df0ef4 100644 --- a/imod/msw/idf_mapping.py +++ b/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) diff --git a/imod/msw/pkgbase.py b/imod/msw/pkgbase.py index 442644fc2..a43dcacba 100644 --- a/imod/msw/pkgbase.py +++ b/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"] diff --git a/imod/tests/test_msw/test_idf_mapping.py b/imod/tests/test_msw/test_idf_mapping.py index db95e8fd1..bcd7b90e7 100644 --- a/imod/tests/test_msw/test_idf_mapping.py +++ b/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) diff --git a/imod/tests/test_msw/test_output_control.py b/imod/tests/test_msw/test_output_control.py index 8cca4ce6c..4641f4712 100644 --- a/imod/tests/test_msw/test_output_control.py +++ b/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 From 93a3ba9990e93e343698b1939b20833ecbc3a284 Mon Sep 17 00:00:00 2001 From: verkaik Date: Tue, 17 Mar 2026 08:30:12 +0100 Subject: [PATCH 2/2] Moved modification of the render function according to PR. --- imod/msw/idf_mapping.py | 41 +++++++++++++++++++++++++++++++++++++++-- imod/msw/pkgbase.py | 24 ------------------------ 2 files changed, 39 insertions(+), 26 deletions(-) diff --git a/imod/msw/idf_mapping.py b/imod/msw/idf_mapping.py index 933df0ef4..a8e46c26f 100644 --- a/imod/msw/idf_mapping.py +++ b/imod/msw/idf_mapping.py @@ -1,15 +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 @@ -132,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) diff --git a/imod/msw/pkgbase.py b/imod/msw/pkgbase.py index a43dcacba..442644fc2 100644 --- a/imod/msw/pkgbase.py +++ b/imod/msw/pkgbase.py @@ -139,30 +139,6 @@ 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"]