From 50ca61219ddd24e5c92ef1a5bdd646d07acb4630 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 25 Mar 2026 12:17:30 +0100 Subject: [PATCH 1/4] Make polygonize work for polygons with holes --- docs/api/changelog.rst | 1 + imod/tests/test_spatial.py | 27 +++++++++++++++++++++++++++ imod/util/spatial.py | 8 +++++++- 3 files changed, 35 insertions(+), 1 deletion(-) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 703a719b7..1e06933d3 100644 --- a/docs/api/changelog.rst +++ b/docs/api/changelog.rst @@ -55,6 +55,7 @@ Fixed layer, when the stage and bottom_elevation were exactly equal to the bottom of a layer in the model discretization, which would cause these cells to be dropped when distributing conductances later. +- Fixed :func:`imod.prepare.spatial.polygonize` for polygons with holes. Changed ~~~~~~~ diff --git a/imod/tests/test_spatial.py b/imod/tests/test_spatial.py index e6e7cc06f..c2a2ea36c 100644 --- a/imod/tests/test_spatial.py +++ b/imod/tests/test_spatial.py @@ -7,6 +7,7 @@ import shapely.geometry as sg import xarray as xr from pytest_cases import parametrize_with_cases +from shapely.geometry import Polygon import imod from imod.testing import assert_frame_equal @@ -118,6 +119,32 @@ def test_polygonize(): assert sorted(gdf["value"]) == [1.0, 2.0, 3.0] +def fits_exact_hole(parent: Polygon, child: Polygon, tol=0.0) -> bool: + hole_polys = Polygon(parent.interiors[0]) + exterior = Polygon(child.exterior) + return exterior.equals(hole_polys) + + +def test_polygonize_with_holes(): + nrow, ncol = 5, 5 + dx, dy = 1.0, -1.0 + xmin, xmax = 0.0, 5.0 + ymin, ymax = 0.0, 5.0 + coords = imod.util.spatial._xycoords((xmin, xmax, ymin, ymax), (dx, dy)) + kwargs = {"name": "test", "coords": coords, "dims": ("y", "x")} + data = np.ones((nrow, ncol), dtype=np.float32) + data[1:-1, 1:-1] = 2.0 + data[2, 2] = 3.0 + da = xr.DataArray(data, **kwargs) + gdf = imod.prepare.polygonize(da) + gdf = gdf.sort_values("value").reset_index(drop=True) + assert len(gdf) == 3 + assert sorted(gdf["value"]) == [1.0, 2.0, 3.0] + p1, p2, p3 = gdf.geometry + ok = fits_exact_hole(p1, p2) and fits_exact_hole(p2, p3) + assert ok + + def test_handle_dtype(): assert imod.prepare.spatial._handle_dtype(np.uint8, None) == (1, 0) assert imod.prepare.spatial._handle_dtype(np.uint16, None) == (2, 0) diff --git a/imod/util/spatial.py b/imod/util/spatial.py index 6f755a502..fd5d23482 100644 --- a/imod/util/spatial.py +++ b/imod/util/spatial.py @@ -730,7 +730,13 @@ def _polygonize(da: xr.DataArray) -> "gpd.GeoDataFrame": geometries = [] colvalues = [] for geom, colval in shapes: - geometries.append(shapely.geometry.Polygon(geom["coordinates"][0])) + coordinates = geom["coordinates"] + exterior = coordinates[0] + if len(coordinates) > 1: + interiors = coordinates[1:] + else: + interiors = None + geometries.append(shapely.geometry.Polygon(exterior, interiors)) colvalues.append(colval) gdf = gpd.GeoDataFrame({"value": colvalues, "geometry": geometries}) From 0adb0399bf1a39d458a7ceb9c2e6ff3f4b741a9d Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 25 Mar 2026 13:49:57 +0100 Subject: [PATCH 2/4] Add description to test --- imod/tests/test_spatial.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/imod/tests/test_spatial.py b/imod/tests/test_spatial.py index c2a2ea36c..8e43d8c9d 100644 --- a/imod/tests/test_spatial.py +++ b/imod/tests/test_spatial.py @@ -126,6 +126,21 @@ def fits_exact_hole(parent: Polygon, child: Polygon, tol=0.0) -> bool: def test_polygonize_with_holes(): + """ + Polygonize should be able to handle polygons with holes. In this test, we + create a 5x5 grid of values as follows: + + 1 1 1 1 1 + 1 2 2 2 1 + 1 2 3 2 1 + 1 2 2 2 1 + 1 1 1 1 1 + + The polygon for the value 1.0 should be a large square that contains the + entire grid, but has a hole in the middle that corresponds to the polygon + for the value 2.0. The polygon for the value 2.0 should have a hole that + corresponds exactly to the polygon for the value 3.0. + """ nrow, ncol = 5, 5 dx, dy = 1.0, -1.0 xmin, xmax = 0.0, 5.0 From 866cc577fef326963abee8937ecfb80c59139856 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 25 Mar 2026 13:50:34 +0100 Subject: [PATCH 3/4] format --- imod/tests/test_spatial.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/tests/test_spatial.py b/imod/tests/test_spatial.py index 8e43d8c9d..44c8c6d6d 100644 --- a/imod/tests/test_spatial.py +++ b/imod/tests/test_spatial.py @@ -135,7 +135,7 @@ def test_polygonize_with_holes(): 1 2 3 2 1 1 2 2 2 1 1 1 1 1 1 - + The polygon for the value 1.0 should be a large square that contains the entire grid, but has a hole in the middle that corresponds to the polygon for the value 2.0. The polygon for the value 2.0 should have a hole that From 52c330003976c7670cc3d3c8b01a167561d13657 Mon Sep 17 00:00:00 2001 From: JoerivanEngelen Date: Wed, 25 Mar 2026 13:52:38 +0100 Subject: [PATCH 4/4] Add comments --- imod/tests/test_spatial.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/imod/tests/test_spatial.py b/imod/tests/test_spatial.py index 44c8c6d6d..7c04f2bf7 100644 --- a/imod/tests/test_spatial.py +++ b/imod/tests/test_spatial.py @@ -141,6 +141,7 @@ def test_polygonize_with_holes(): for the value 2.0. The polygon for the value 2.0 should have a hole that corresponds exactly to the polygon for the value 3.0. """ + # Arrange nrow, ncol = 5, 5 dx, dy = 1.0, -1.0 xmin, xmax = 0.0, 5.0 @@ -151,8 +152,10 @@ def test_polygonize_with_holes(): data[1:-1, 1:-1] = 2.0 data[2, 2] = 3.0 da = xr.DataArray(data, **kwargs) + # Act gdf = imod.prepare.polygonize(da) gdf = gdf.sort_values("value").reset_index(drop=True) + # Assert assert len(gdf) == 3 assert sorted(gdf["value"]) == [1.0, 2.0, 3.0] p1, p2, p3 = gdf.geometry