diff --git a/src/xtgeo/surface/_regsurf_boundary.py b/src/xtgeo/surface/_regsurf_boundary.py index 3e1bd0a7b..05dced1bd 100644 --- a/src/xtgeo/surface/_regsurf_boundary.py +++ b/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) diff --git a/src/xtgeo/xyz/_polygons_oper.py b/src/xtgeo/xyz/_polygons_oper.py index e7fafc7e0..d710ccbd7 100644 --- a/src/xtgeo/xyz/_polygons_oper.py +++ b/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: diff --git a/tests/test_surface/test_regular_surface.py b/tests/test_surface/test_regular_surface.py index fe1c0f34e..cbec51803 100644 --- a/tests/test_surface/test_regular_surface.py +++ b/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."""