FIX: inrease robustness in boundary polygons (#1163)

src/xtgeo/surface/_regsurf_boundary.py

@@ -15,9 +15,6 @@ def create_boundary(self, alpha_factor, is_convex, simplify):
 
     pol = Polygons.boundary_from_points(points, alpha_factor, alpha, is_convex)
 
-    if pol is None:
-        raise RuntimeError("Could not create a valid Polygons instance for boundary.")
-
     if simplify:
         if isinstance(simplify, bool):
             pol.simplify(tolerance=0.1)

src/xtgeo/xyz/_polygons_oper.py

@@ -9,10 +9,13 @@
 Functions starting with '_' are local helper functions
 """
 
+from __future__ import annotations
+
 from math import ceil
 
 import numpy as np
 import pandas as pd
+import shapely
 import shapely.geometry as sg
 from scipy.spatial import Delaunay, cKDTree
 from shapely.ops import polygonize
@@ -121,12 +124,13 @@ def _create_boundary_polygon(
     yvec: np.ndarray,
     zvec: np.ndarray,
     alpha: float,
-):
+) -> list[list[float | int]] | None:
     xy = np.column_stack((xvec, yvec))
     coords = np.column_stack((xy, zvec))
 
     # return None if too few points
     if len(xy) < MINIMUM_NUMBER_POINTS:
+        logger.info("Too few points to derive boundary, need at least 4")
         return None
 
     edges = _alpha_shape(xy, alpha=alpha)
@@ -135,35 +139,33 @@ def _create_boundary_polygon(
             "Your alpha or alpha_factor value is too low, try increasing it!"
         )
 
-    polygons = _get_shapely_polygons_from_edges(coords, edges)
+    pol_coords = _get_polygon_coords_from_edges_shapely(coords, edges)
+    if not pol_coords:
+        raise BoundaryError("Could not create a valid Polygons instance for boundary.")
 
     # sort polygons by their size to get the largest polygons first
-    sorted(polygons, key=lambda pol: len(pol.exterior.coords), reverse=True)
+    pol_coords = sorted(pol_coords, key=len, reverse=True)
 
     data = []
-    for polid, pol in enumerate(polygons):
-        for coord in pol.exterior.coords:
+    for polid, polcoord in enumerate(pol_coords):
+        for coord in polcoord:
             data.append(list(coord) + [polid])
 
     return data
 
 
-def _alpha_shape(points, alpha):
+def _alpha_shape(points: np.ndarray, alpha: float) -> set[tuple[int, int]]:
     """Compute the alpha shape (concave hull) of a set of points.
 
     Args:
         points: np.array of shape (n,2) points.
         alpha: alpha value.
-        only_outer TODO?: boolean value to specify if we keep only the outer border
-            or also inner edges.
     Returns:
         Set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
             the indices in the points array.
     """
 
-    assert points.shape[0] >= MINIMUM_NUMBER_POINTS, "Need >= 4 pts to derive boundary"
-
-    def add_edge(edges, icv, jcv):
+    def add_edge(edges: set[tuple[int, int]], icv: int, jcv: int) -> None:
         """Add an edge between the i-th and j-th points, if not in the list already."""
         if (icv, jcv) in edges or (jcv, icv) in edges:
             # if both neighboring triangles are in shape, it is not a boundary edge
@@ -191,22 +193,30 @@ def add_edge(edges, icv, jcv):
         partial = partial[partial > 0.0]  # to avoid sqrt of negative number
         area = np.sqrt(partial)
 
-        # radius of circumcircle
-        circum_r = avv * bvv * cvv / (4.0 * area) if area.size > 0 else 0
+        if area.size > 0:
+            # radius of circumcircle
+            circum_r = avv * bvv * cvv / (4.0 * area)
 
-        # if radius less then alpha then add outer edge
-        if circum_r < alpha or alpha == 0:
-            add_edge(edges, ia, ib)
-            add_edge(edges, ib, ic)
-            add_edge(edges, ic, ia)
+            # if radius less then alpha then add outer edge
+            if circum_r < alpha or alpha == 0:
+                add_edge(edges, ia, ib)
+                add_edge(edges, ib, ic)
+                add_edge(edges, ic, ia)
 
     return edges
 
 
-def _get_shapely_polygons_from_edges(coords, edges):
+def _get_polygon_coords_from_edges_shapely(
+    coords: np.ndarray, edges: set[tuple[int, int]]
+) -> list[shapely.coords.CoordinateSequence]:
     """Split and connect edges into shapely polygons"""
+    logger.info("Getting polygon coordinates using shapely")
     lines = sg.MultiLineString([(coords[e[0]], coords[e[1]]) for e in edges])
-    return polygonize(lines)
+    # Need to use node workaround here to handle crossing lines
+    # see https://github.com/shapely/shapely/issues/1736
+    noded = shapely.node(shapely.GeometryCollection(lines))
+    polygons = polygonize(noded.geoms)
+    return [pol.exterior.coords for pol in polygons]
 
 
 def simplify_polygons(self, tolerance: float, preserve_topology: bool) -> bool:

tests/test_surface/test_regular_surface.py

@@ -1321,10 +1321,10 @@ def test_get_boundary_polygons_complex(show_plot):
     # for some reasons, macos tests gives slightly different result; that is why a large
     # tolerance is given
     assert boundary.get_dataframe()[boundary.xname].mean() == pytest.approx(
-        462392.0, abs=3.0
+        462230.0, abs=3.0
     )
     assert boundary.get_dataframe()[boundary.yname].mean() == pytest.approx(
-        5933152.0, abs=4.0
+        5933457.0, abs=4.0
     )
 
     # get the first (largest) polygon
@@ -1344,6 +1344,26 @@ def test_get_boundary_polygons_complex(show_plot):
         )
 
 
+def test_boundary_polygons_are_sorted():
+    """Test that boundary polygons are sorted from largest to smallest."""
+    xs = xtgeo.surface_from_file(TESTSET1)
+    xs.values = np.ma.masked_less(xs.values, 1700)
+    xs.values = np.ma.masked_greater(xs.values, 1800)
+
+    boundary = xs.get_boundary_polygons(simplify=False)
+
+    df = boundary.get_dataframe(copy=False)
+
+    # check that we have 7 unique boundaries for this surface
+    assert df["POLY_ID"].nunique() == 7
+
+    # check that the boundary are sorted from largest to smallest polygon
+    pol_lengths = [len(poldf) for _, poldf in df.groupby("POLY_ID")]
+    assert all(
+        pol_lengths[i] >= pol_lengths[i + 1] for i in range(len(pol_lengths) - 1)
+    )
+
+
 def test_regsurface_get_dataframe(default_surface):
     """Test getting a dataframe from a surface."""