ENH: Add grid boundary polygons (equinor#1169)

docs/tutorial/examples_rms.rst

@@ -277,6 +277,54 @@ Edit a 3D grid porosity inside polygons
    # Save in RMS as a new icon
    myprop.to_roxar(project, "Reek_sim", "NEWPORO_setinside")
 
+   
+Create region polygons from the grid
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. code-block:: python
+
+   import numpy as np
+   import xtgeo
+
+   GNAME = "Simgrid"
+   REGNAME = "Regions"
+   ZONENAME = "Zone"
+
+   CB_FOLDER = "Region_polygons"
+   
+   ZONE_FILTER = [2, 3]
+
+   # factor that controls the precision of the polygons
+   # higher value gives smoother more convex polygons.
+   ALPHA_FACTOR = 1 
+
+   def create_region_polygons():
+       """Create region polygons and store them on the clipboard"""
+       grid = xtgeo.grid_from_roxar(project, GNAME)
+       reg = xtgeo.gridproperty_from_roxar(project, GNAME, REGNAME)
+       zone = xtgeo.gridproperty_from_roxar(project, GNAME, ZONENAME)
+
+       for regnum, regname in reg.codes.items():
+           print(f"Creating boundary polygon for region {regname}")
+
+           # create a filter array to extract boundaries for the region
+           # zone was used as filter to minimice overlap of the final polygons
+           # Note: a layer filter could have been applied instead
+           filter_array = (reg.values==regnum) & (np.isin(zone.values, [ZONE_FILTER]))
+
+           pol = grid.get_boundary_polygons(ALPHA_FACTOR, filter_array=filter_array)
+
+           # in case of several polygons, keep only the largest (first)
+           pol.filter_byid([0])
+
+           # store polygon to the clipboard
+           pol.to_roxar(project, regname, CB_FOLDER, stype="clipboard")
+
+        print(f"Complete, region polygons are stored under clipboard folder {CB_FOLDER}")
+
+    if __name__ == "__main__":
+        create_region_polygons()
+
 .. _hybrid:
 
 Make a hybrid grid

src/xtgeo/grid3d/_grid_boundary.py

@@ -0,0 +1,76 @@
+from __future__ import annotations
+
+import math
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
+
+from xtgeo.xyz.points import Points
+from xtgeo.xyz.polygons import Polygons
+
+if TYPE_CHECKING:
+    from xtgeo.grid3d import Grid
+
+
+def create_boundary(
+    grid: Grid,
+    alpha_factor: float = 1.0,
+    convex: bool = False,
+    simplify: bool | dict[str, Any] = True,
+    filter_array: np.ndarray | None = None,
+) -> Polygons:
+    """Create boundary polygons for a grid."""
+
+    xval, yval, zval = (prop.values for prop in grid.get_xyz())
+
+    if filter_array is not None:
+        if filter_array.shape != grid.dimensions:
+            raise ValueError(
+                "The filter_array needs to have the same dimensions as the grid. "
+                f"Found: {filter_array.shape=} {grid.dimensions=}"
+            )
+        xval = np.ma.masked_where(~filter_array, xval)
+        yval = np.ma.masked_where(~filter_array, yval)
+        zval = np.ma.masked_where(~filter_array, zval)
+
+    # for performance create average points along layers
+    xval = np.ma.mean(xval, axis=2)
+    yval = np.ma.mean(yval, axis=2)
+    zval = np.ma.mean(zval, axis=2)
+
+    xyz_values = np.column_stack(
+        (
+            xval[~xval.mask].ravel(),
+            yval[~yval.mask].ravel(),
+            zval[~zval.mask].ravel(),
+        )
+    )
+
+    pol = Polygons.boundary_from_points(
+        points=Points(xyz_values),
+        alpha_factor=alpha_factor,
+        alpha=_estimate_alpha_for_grid(grid),
+        convex=convex,
+    )
+
+    if simplify:
+        if isinstance(simplify, bool):
+            pol.simplify(tolerance=0.1)
+        elif isinstance(simplify, dict) and "tolerance" in simplify:
+            pol.simplify(**simplify)
+        else:
+            raise ValueError("Invalid values for simplify keyword")
+
+    return pol
+
+
+def _estimate_alpha_for_grid(grid: Grid) -> float:
+    """
+    Estimate an alpha based on grid resolution.
+    Max dx and dy is used as basis for calculation to ensure that the alpha
+    computed is always high enough to prevent polygons appearing around areas
+    of the grid where cells have larger than average dx/dy increments.
+    """
+    dx, dy = grid.get_dx(), grid.get_dy()
+    xinc, yinc = dx.values.max(), dy.values.max()
+    return math.ceil(math.sqrt(xinc**2 + yinc**2) / 2)

src/xtgeo/grid3d/_grid_etc1.py

@@ -3,6 +3,7 @@
 from __future__ import annotations
 
 from copy import deepcopy
+from functools import lru_cache
 from math import atan2, degrees
 from typing import TYPE_CHECKING, Literal
 
@@ -134,6 +135,7 @@ def get_dz(
     )
 
 
+@lru_cache(maxsize=1)
 def get_dx(
     self: Grid, name: str = "dX", asmasked: bool = False, metric: METRIC = "horizontal"
 ) -> GridProperty:
@@ -166,6 +168,7 @@ def get_dx(
     )
 
 
+@lru_cache(maxsize=1)
 def get_dy(
     self: Grid, name: str = "dY", asmasked: bool = False, metric: METRIC = "horizontal"
 ) -> GridProperty:
@@ -384,6 +387,7 @@ def get_ijk_from_points(
     return list(mydataframe.itertuples(index=False, name=None))
 
 
+@lru_cache(maxsize=1)
 def get_xyz(
     self: Grid,
     names: tuple[str, str, str] = ("X_UTME", "Y_UTMN", "Z_TVDSS"),

src/xtgeo/grid3d/grid.py

@@ -22,6 +22,7 @@
 
 from . import (
     _grid3d_fence,
+    _grid_boundary,
     _grid_etc1,
     _grid_export,
     _grid_hybrid,
@@ -517,6 +518,10 @@ def __str__(self) -> str:
 
         return self.describe(flush=False) or ""
 
+    def __hash__(self):
+        """The __hash__ method."""
+        return hash(self.generate_hash())
+
     # ==================================================================================
     # Public Properties:
     # ==================================================================================
@@ -1567,6 +1572,62 @@ def get_zoneprop_from_subgrids(self) -> NoReturn:
         """Make a XTGeo GridProperty instance for a Zone property subgrid index."""
         raise NotImplementedError("Not yet; todo")
 
+    def get_boundary_polygons(
+        self: Grid,
+        alpha_factor: float = 1.0,
+        convex: bool = False,
+        simplify: bool | dict[str, Any] = True,
+        filter_array: np.ndarray | None = None,
+    ) -> Polygons:
+        """Extract boundary polygons from the grid cell centers.
+
+        A ``filter_array`` can be applied to extract boundaries around specific
+        parts of the grid e.g. a region or a zone.
+
+        The concavity and precision of the boundaries are controlled by the
+        ``alpha_factor``. A low ``alpha_factor`` makes more precise boundaries,
+        while a larger value makes more rough polygons.
+
+        Note that the ``alpha_factor`` is a multiplier (default value 1) on top
+        of an auto estimated value, derived from the maximum xinc and yinc from
+        the grid cells. Dependent on the regularity of the grid, tuning of the
+        alpha_factor (up/down) is sometimes necessary to get satisfactory results.
+
+        Args:
+            alpha_factor: An alpha multiplier, which controls the precision of the
+                boundaries. A higher number will produce smoother and less accurate
+                polygons. Not applied if convex is set to True.
+            convex: The default is False, which means that a "concave hull" algorithm
+                is used. If convex is True, the alpha factor is overridden to a large
+                number, producing a 'convex' shape boundary instead.
+            simplify: If True, a simplification is done in order to reduce the number
+                of points in the polygons, where tolerance is 0.1. Another
+                alternative to True is to input a Dict on the form
+                ``{"tolerance": 2.0, "preserve_topology": True}``, cf. the
+                :func:`Polygons.simplify()` method. For details on e.g. tolerance, see
+                Shapely's simplify() method.
+            filter_array: An numpy boolean array with equal shape as the grid dimension,
+                used to filter the grid cells and define where to extract boundaries.
+
+        Returns:
+            A XTGeo Polygons instance
+
+        Example::
+
+            grid = xtgeo.grid_from_roxar(project, "Simgrid")
+            # extract polygon for a specific region, here region 3
+            region = xtgeo.gridproperty_from_roxar(project, "Simgrid", "Regions")
+            filter_array = (region.values==3)
+            boundary = grid.get_boundary_polygons(filter_array=filter_array)
+
+        See also:
+            The :func:`Polygons.boundary_from_points()` class method.
+
+        """
+        return _grid_boundary.create_boundary(
+            self, alpha_factor, convex, simplify, filter_array
+        )
+
     def get_actnum_indices(
         self,
         order: Literal["C", "F", "A", "K"] = "C",
@@ -2047,7 +2108,7 @@ def get_xyz(
             )
             asmasked = self._evaluate_mask(mask)
 
-        return _grid_etc1.get_xyz(self, names=names, asmasked=asmasked)
+        return _grid_etc1.get_xyz(self, names=tuple(names), asmasked=asmasked)
 
     def get_xyz_cell_corners(
         self,

tests/test_grid3d/test_grid_operations.py

@@ -278,3 +278,36 @@ def test_reduce_to_one_layer(grd):
     grd.reduce_to_one_layer()
 
     assert grd.nlay == 1
+
+
+def test_get_boundary_polygons_grid_outline():
+    """Test getting a boundary for a grid."""
+
+    grid = xtgeo.grid_from_file(EMEGFILE)
+
+    boundary = grid.get_boundary_polygons(simplify=False)
+    df = boundary.get_dataframe(copy=False)
+
+    assert df["POLY_ID"].nunique() == 1
+    assert df[boundary.yname].min() == pytest.approx(5930003, abs=2)
+    assert df[boundary.yname].max() == pytest.approx(5937505, abs=2)
+    assert df[boundary.xname].min() == pytest.approx(459078, abs=2)
+    assert df[boundary.xname].max() == pytest.approx(466984, abs=2)
+
+
+def test_get_boundary_polygons_grid_region():
+    """Test getting a boundary for a region in a grid."""
+
+    grid = xtgeo.grid_from_file(EMEGFILE)
+    reg = xtgeo.gridproperty_from_file(EMERFILE, name="REGION")
+
+    boundary = grid.get_boundary_polygons(
+        simplify=False, filter_array=(reg.values == 1)
+    )
+    df = boundary.get_dataframe(copy=False)
+
+    assert df["POLY_ID"].nunique() == 1
+    assert df[boundary.yname].min() == pytest.approx(5932092, abs=2)
+    assert df[boundary.yname].max() == pytest.approx(5936223, abs=2)
+    assert df[boundary.xname].min() == pytest.approx(460344, abs=2)
+    assert df[boundary.xname].max() == pytest.approx(463765, abs=2)