Allow polygon shape to consist of disjoint parts

src/ert/config/_shapes.py

@@ -1,7 +1,7 @@
 from __future__ import annotations
 
 from functools import cached_property
-from typing import Annotated, ClassVar, Literal, Self
+from typing import Annotated, ClassVar, Literal, Self, cast
 
 import shapely
 import xtgeo
@@ -41,16 +41,16 @@ def __eq__(self, other: object) -> bool:
 
 
 class PolygonShapeConfig(ShapeConfig):
-    """Configuration for a polygonal shape.
+    """Configuration for a (multi)polygonal shape.
 
     Attributes:
-        vertices: List of (east, north) tuples defining the polygon vertices. Vertices
-        are expected to be normalized in shapely's sense (first vertex is the lowest,
-        and vertices are ordered clockwise).
+        wkt: Well-Known Text representation of the multipolygon. Vertices are expected
+        to be normalized in shapely's sense (first vertex is the lowest, and vertices
+        are ordered clockwise).
     """
 
     type: Literal["polygon"] = "polygon"
-    vertices: list[tuple[float, float]]
+    wkt: str  # Well-Known Text representation of the multipolygon
 
     TOLERANCE: ClassVar[float] = 0.1
 
@@ -62,37 +62,48 @@ def from_file(cls, filepath: str) -> Self:
             filepath: Path to a file containing polygon definition. Supported formats
             are the ones supported by xtgeo.polygons_from_file
             https://xtgeo.readthedocs.io/en/latest/api-points-polygons.html#xtgeo.polygons_from_file.
-            Expected to contain exactly one polygon with no holes.
+            Expected to contain one or more polygons with no holes. Multiple polygons
+            (disjoint or overlapping) are preserved; overlapping polygons are merged.
         """
-        xtgeo_polygon = xtgeo.polygons_from_file(filepath)
-
-        if len(set(xtgeo_polygon.dataframe["POLY_ID"])) != 1:
-            raise ValueError(
-                "Multiple polygons found in the file. "
-                "Behavior is defined only for one polygon."
-            )
+        xtgeo_polygons = xtgeo.polygons_from_file(filepath)
+        line_strings = xtgeo_polygons.get_shapely_objects()
 
         try:
-            shapely_polygon = (
-                shapely.Polygon(xtgeo_polygon.get_xyz_arrays())
-                .simplify(tolerance=cls.TOLERANCE)
-                .normalize()
-            )
+            separate_polygons = [shapely.Polygon(poly.coords) for poly in line_strings]
+            polygon_union = shapely.union_all(separate_polygons)
         except Exception as e:
             raise ValueError(f"Failed to create polygon from file {filepath}") from e
 
-        vertices = shapely.get_coordinates(shapely_polygon).tolist()
+        if isinstance(polygon_union, shapely.MultiPolygon):
+            multipolygon = polygon_union
+        elif isinstance(polygon_union, shapely.Polygon):
+            multipolygon = shapely.MultiPolygon([polygon_union])
+        else:
+            raise ValueError(
+                f"Shapes in the file '{filepath}' could not be converted to polygons. "
+                f"Unexpected geometry type {type(polygon_union).__name__}"
+            )
+
+        cleaned_polygons = [
+            cast(
+                shapely.Polygon,
+                geom.simplify(tolerance=cls.TOLERANCE).normalize(),
+            )
+            for geom in multipolygon.geoms
+        ]
+        multipolygon = shapely.MultiPolygon(cleaned_polygons)
 
-        return cls(vertices=vertices)
+        return cls(wkt=multipolygon.wkt)
 
     @cached_property
-    def _polygon(self) -> shapely.Polygon:
-        poly = shapely.Polygon(self.vertices)
-        shapely.prepare(poly)
-        return poly
+    def _polygon(self) -> shapely.MultiPolygon:
+        multipolygon = shapely.from_wkt(self.wkt)
+        assert isinstance(multipolygon, shapely.MultiPolygon)
+        shapely.prepare(multipolygon)
+        return multipolygon
 
     def contains(self, east: float, north: float) -> bool:
-        """Check if a point is inside the polygon.
+        """Check if a point is inside any of the internal polygons.
 
         Args:
             east: UTM_X-coordinate of the point

tests/ert/unit_tests/config/test_observation_declaration.py

@@ -1248,14 +1248,8 @@ def test_that_seismic_observation_reads_boundary_file(file_context_token):
     ]
     boundary = shape_registry.get(0)
     assert isinstance(boundary, PolygonShapeConfig)
-    expected = [
-        (0.0, 0.0),
-        (0.0, 1.0),
-        (1.0, 1.0),
-        (1.0, 0.0),
-        (0.0, 0.0),
-    ]
-    assert boundary.vertices == expected
+    expected = "MULTIPOLYGON Z (((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0)))"
+    assert boundary.wkt == expected
 
 
 @pytest.mark.usefixtures("use_tmpdir")

tests/ert/unit_tests/config/test_observation_quality_control.py

@@ -174,15 +174,9 @@ def test_that_observation_without_zones_are_not_disabled_by_zone_check():
 
 
 def test_that_observations_within_boundary_stay_while_outside_are_removed():
-    boundary = PolygonShapeConfig(
-        vertices=[
-            (0.0, 0.0),
-            (0.0, 1.0),
-            (1.0, 1.0),
-            (1.0, 0.0),
-            (0.0, 0.0),
-        ]
-    )
+    polygon1 = "((0 0 0, 0 1 0, 1 1 0, 1 0 0, 0 0 0))"
+    polygon2 = "((3 3 0, 3 4 0, 4 4 0, 4 3 0, 3 3 0))"
+    boundary = PolygonShapeConfig(wkt=f"MULTIPOLYGON Z ({polygon1}, {polygon2})")
     shape_registry = ShapeRegistry()
     boundary_id = shape_registry.register(boundary)
 
@@ -207,8 +201,13 @@ def test_that_observations_within_boundary_stay_while_outside_are_removed():
                 east=2.5,
                 north=2.5,
             ),
+            create_seismic_observation(
+                east=3.5,
+                north=3.5,
+                boundary_id=boundary_id,
+            ),
         ]
     )
 
     qc = qc_seismic_observations(observations, shape_registry)
-    assert qc["east"].to_list() == [0.5, 2.5]
+    assert qc["east"].to_list() == [0.5, 2.5, 3.5]