diff --git a/docs/usextgeoroxar.rst b/docs/usextgeoroxar.rst index ca8a88a3a..141051752 100644 --- a/docs/usextgeoroxar.rst +++ b/docs/usextgeoroxar.rst @@ -521,4 +521,39 @@ be input to Equinor's APS module. Line point data --------------- -Examples to come... +Add to or remove points inside or outside polygons +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In the following example, remove or add to points being inside or outside polygons on clipboard. + +.. code-block:: python + + import xtgeo + + PRJ = project + + POLYGONS = ["mypolygons", "myfolder"] # mypolygons in folder myfolder on clipboard + POINTSET1 = ["points1", "myfolder"] + POINTSET2 = ["points2", "myfolder"] + + POINTSET1_UPDATED = ["points1_edit", "myfolder"] + POINTSET2_UPDATED = ["points2_edit", "myfolder"] + + def main(): + """Operations on points inside or outside polygons.""" + + poly = xtgeo.polygons_from_roxar(PRJ, *POLYGONS, stype="clipboard") + po1 = xtgeo.points_from_roxar(PRJ, *POINTSET1, stype="clipboard") + po2 = xtgeo.points_from_roxar(PRJ, *POINTSET2, stype="clipboard") + + po1.eli_inside_polygons(poly) + po1.to_roxar(PRJ, *POINTSET1_UPDATED, stype="clipboard") # store + + # now add 100 inside polugons for POINTSET2, and then remove all points outside + po2.add_inside_polygons(poly, 100) + po2.eli_outside_polygons(poly) + po2.to_roxar(PRJ, *POINTSET2_UPDATED, stype="clipboard") # store + + + if __name__ == "__main__": + main() diff --git a/src/xtgeo/surface/_regsurf_oper.py b/src/xtgeo/surface/_regsurf_oper.py index 1bf8c6b81..a68f23fcd 100644 --- a/src/xtgeo/surface/_regsurf_oper.py +++ b/src/xtgeo/surface/_regsurf_oper.py @@ -2,9 +2,11 @@ """Various operations""" import numbers +from typing import Any, Union import numpy as np import numpy.ma as ma +from matplotlib.path import Path as MPath import xtgeo import xtgeo.cxtgeo._cxtgeo as _cxtgeo # type: ignore @@ -452,6 +454,7 @@ def _get_randomline_fence(self, fencespec, hincrement, atleast, nextend): def operation_polygons(self, poly, value, opname="add", inside=True): """Operations restricted to polygons""" + # keep this for a while (e.g. mid 2024), and then replace it with _v2 below. if not isinstance(poly, Polygons): raise ValueError("The poly input is not a Polygons instance") if opname not in VALID_OPER_POLYS: @@ -469,7 +472,7 @@ def operation_polygons(self, poly, value, opname="add", inside=True): if isinstance(value, type(self)): if not self.compare_topology(value): - raise ValueError("Input is RegularSurface, but not same map " "topology") + raise ValueError("Input is RegularSurface, but not same map topology") value = value.values.copy() else: # turn scalar value into numpy array @@ -518,7 +521,7 @@ def operation_polygons(self, poly, value, opname="add", inside=True): # as result in these cases if 0.0 in value: xtg.warn( - "Dividing a surface with value or surface with zero " + "Dividing a surface with value=0.0 or surface with zero " "elements; may get unexpected results, try to " "achieve zero values as result!" ) @@ -539,3 +542,95 @@ def operation_polygons(self, poly, value, opname="add", inside=True): self.values[proxyv == proxytarget] = tmp[proxyv == proxytarget] del tmp + + +def _proxy_map_polygons(surf, poly, inside=True): + """Return a proxy map where on one to do operations, as 0 and 1.""" + inside_value = 1 if inside else 0 + outside_value = 0 if inside else 1 + + proxy = surf.copy() + + # allow a single Polygons instance or a list of Polygons instances + if isinstance(poly, xtgeo.Polygons): + usepolys = [poly] + elif isinstance(poly, list) and all( + isinstance(pol, xtgeo.Polygons) for pol in poly + ): + usepolys = poly + else: + raise ValueError("The poly values is not a Polygons or a list of Polygons") + + proxy.values = outside_value + xvals, yvals = proxy.get_xy_values(asmasked=False) + points = np.array([xvals.ravel(), yvals.ravel()]).T + + for pol in usepolys: + idgroups = pol.dataframe.groupby(pol.pname) + for _, grp in idgroups: + singlepoly = np.array([grp[pol.xname].values, grp[pol.yname].values]).T + poly_path = MPath(singlepoly) + is_inside = poly_path.contains_points(points) + is_inside = is_inside.reshape(proxy.ncol, proxy.nrow) + proxy.values = np.where(is_inside, inside_value, proxy.values) + + return proxy + + +def operation_polygons_v2( + self, poly, value: Union[float, Any], opname="add", inside=True +): + """Operations restricted to polygons, using matplotlib (much faster). + + The 'value' can be a number or another regular surface (with same design) + + """ + + proxy = _proxy_map_polygons(self, poly, inside=inside) + result = self.copy() + + if isinstance(value, type(result)): + if not result.compare_topology(value): + raise ValueError("Input is RegularSurface, but not same map topology") + value = value.values.copy() + else: + # turn scalar value into numpy array + value = result.values.copy() * 0 + value + + if opname == "add": + result.values = np.ma.where( + proxy.values == 1, result.values + value, result.values + ) + elif opname == "sub": + result.values = np.ma.where( + proxy.values == 1, result.values - value, result.values + ) + elif opname == "mul": + result.values = np.ma.where( + proxy.values == 1, result.values * value, result.values + ) + elif opname == "div": + # Dividing a map of zero is always a hazzle; try to obtain 0.0 + # as result in these cases + if 0.0 in value: + xtg.warn( + "Dividing a surface with value = 0.0 or surface with zero " + "elements; may get unexpected results, will try to " + "achieve zero values as result!" + ) + + result.values = np.ma.where(value == 0.0, 0.0, result.values) + proxy.values = np.ma.where(value == 0.0, 0, proxy.values) + result.values = np.ma.where( + proxy.values == 1, result.values / value, result.values + ) + + elif opname == "set": + result.values = np.ma.where(proxy.values == 1, value, result.values) + elif opname == "eli": + result.values = np.ma.where(proxy.values == 1, xtgeo.UNDEF, result.values) + result.values = ma.masked_greater(result.values, xtgeo.UNDEF_LIMIT) + else: + raise KeyError(f"The opname={opname} is not one of {VALID_OPER_POLYS}") + + self.values = result.values diff --git a/src/xtgeo/surface/regular_surface.py b/src/xtgeo/surface/regular_surface.py index f1b5238c1..01dfc680b 100644 --- a/src/xtgeo/surface/regular_surface.py +++ b/src/xtgeo/surface/regular_surface.py @@ -2043,7 +2043,7 @@ def operation(self, opname, value): # Operations restricted to inside/outside polygons # ================================================================================== - def operation_polygons(self, poly, value, opname="add", inside=True): + def operation_polygons(self, poly, value, opname="add", inside=True, _version=2): """A generic function for map operations inside or outside polygon(s). Args: @@ -2051,10 +2051,18 @@ def operation_polygons(self, poly, value, opname="add", inside=True): value(float or RegularSurface): Value to add, subtract etc opname (str): Name of operation... 'add', 'sub', etc inside (bool): If True do operation inside polygons; else outside. + _version (int): Algorithm version, 2 will be much faster when many points + on polygons (this key will be removed in later versions and shall not + be applied) """ - _regsurf_oper.operation_polygons( - self, poly, value, opname=opname, inside=inside - ) + if _version == 2: + _regsurf_oper.operation_polygons_v2( + self, poly, value, opname=opname, inside=inside + ) + else: + _regsurf_oper.operation_polygons( + self, poly, value, opname=opname, inside=inside + ) # shortforms def add_inside(self, poly, value): diff --git a/src/xtgeo/xyz/_xyz.py b/src/xtgeo/xyz/_xyz.py index 0158b6a13..c69d0a783 100644 --- a/src/xtgeo/xyz/_xyz.py +++ b/src/xtgeo/xyz/_xyz.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- """XTGeo XYZ module (abstract base class)""" from abc import ABC, abstractmethod -from typing import List, Union +from typing import List, TypeVar, Union from warnings import warn import numpy as np @@ -14,6 +14,8 @@ xtg = XTGeoDialog() logger = xtg.functionlogger(__name__) +Polygons = TypeVar("Polygons") + class XYZ(ABC): """Abstract base class for XYZ objects, i.e. Points and Polygons in XTGeo. @@ -275,7 +277,7 @@ def get_boundary(self): def mark_in_polygons( self, - poly: Union["Polygons", List["Polygons"]], # noqa: F821 + poly: Union[Polygons, List[Polygons]], # noqa: F821 name: str = "pstatus", inside_value: int = 1, outside_value: int = 0, @@ -297,7 +299,7 @@ def mark_in_polygons( def operation_polygons( self, - poly: Union["Polygons", List["Polygons"]], # noqa: F821 + poly: Union[Polygons, List[Polygons]], # noqa: F821 value: float, opname: str = "add", inside: bool = True, @@ -342,7 +344,7 @@ def operation_polygons( overlapping part. Similarly, using e.g. "eli, outside" will completely remove all points of two non-overlapping polygons are given as input. - ``version=2``: This is a new and recommended implemenation. It works + ``version=2``: This is a new and recommended implementation. It works much faster and intuitively for both inside and outside, overlapping and multiple polygons within a Polygons instance. @@ -384,7 +386,7 @@ def add_inside(self, poly, value): self.operation_polygons(poly, value, opname="add", inside=True, version=0) def add_inside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Add a value (scalar) to points inside polygons (new behaviour). @@ -413,7 +415,7 @@ def add_outside(self, poly, value): self.operation_polygons(poly, value, opname="add", inside=False, version=0) def add_outside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Add a value (scalar) to points outside polygons (new behaviour). @@ -442,7 +444,7 @@ def sub_inside(self, poly, value): self.operation_polygons(poly, value, opname="sub", inside=True, version=1) def sub_inside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Subtract a value (scalar) for points inside polygons (new behaviour). @@ -471,7 +473,7 @@ def sub_outside(self, poly, value): self.operation_polygons(poly, value, opname="sub", inside=False, version=0) def sub_outside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Subtract a value (scalar) for points outside polygons (new behaviour). @@ -500,7 +502,7 @@ def mul_inside(self, poly, value): self.operation_polygons(poly, value, opname="mul", inside=True, version=0) def mul_inside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Multiply a value (scalar) for points inside polygons (new behaviour). @@ -529,7 +531,7 @@ def mul_outside(self, poly, value): self.operation_polygons(poly, value, opname="mul", inside=False, version=0) def mul_outside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Multiply a value (scalar) for points outside polygons (new behaviour). @@ -558,7 +560,7 @@ def div_inside(self, poly, value): self.operation_polygons(poly, value, opname="div", inside=True, version=0) def div_inside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Divide a value (scalar) for points inside polygons (new behaviour). @@ -587,7 +589,7 @@ def div_outside(self, poly, value): self.operation_polygons(poly, value, opname="div", inside=False, version=0) def div_outside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Divide a value (scalar) for points outside polygons (new behaviour). @@ -618,7 +620,7 @@ def set_inside(self, poly, value): self.operation_polygons(poly, value, opname="set", inside=True, version=0) def set_inside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Set a value (scalar) for points inside polygons (new behaviour). @@ -647,7 +649,7 @@ def set_outside(self, poly, value): self.operation_polygons(poly, value, opname="set", inside=False, version=0) def set_outside_polygons( - self, poly: Union["Polygons", List["Polygons"]], value: float # noqa: F821 + self, poly: Union[Polygons, List[Polygons]], value: float # noqa: F821 ): """Set a value (scalar) for points outside polygons (new behaviour). @@ -674,9 +676,7 @@ def eli_inside(self, poly): """ self.operation_polygons(poly, 0, opname="eli", inside=True, version=0) - def eli_inside_polygons( - self, poly: Union["Polygons", List["Polygons"]] # noqa: F821 - ): + def eli_inside_polygons(self, poly: Union[Polygons, List[Polygons]]): # noqa: F821 """Remove points inside polygons. This is an improved implementation than :meth:`eli_inside()`, and is now the @@ -701,9 +701,7 @@ def eli_outside(self, poly): """ self.operation_polygons(poly, 0, opname="eli", inside=False, version=0) - def eli_outside_polygons( - self, poly: Union["Polygons", List["Polygons"]] # noqa: F821 - ): + def eli_outside_polygons(self, poly: Union[Polygons, List[Polygons]]): # noqa: F821 """Remove points outside polygons. This is an improved implementation than :meth:`eli_outside()`, and is now the diff --git a/tests/test_surface/test_regular_surface_in_other.py b/tests/test_surface/test_regular_surface_in_other.py index cd128543f..df3c9bce0 100644 --- a/tests/test_surface/test_regular_surface_in_other.py +++ b/tests/test_surface/test_regular_surface_in_other.py @@ -2,7 +2,10 @@ from os.path import join +import numpy as np +import numpy.ma as ma import pytest + import xtgeo from xtgeo.common import XTGeoDialog @@ -21,18 +24,129 @@ SURF1 = TPATH / "surfaces/reek/1/topreek_rota.gri" POLY1 = TPATH / "polygons/reek/1/closedpoly1.pol" - -def test_operations_inside_outside_polygon_generic(): - """Do perations in relation to polygons in a generic way""" +# a surface with nodes in X position 0.0 3.0 6.0 and Y 0.0 3.0 6.0 +# and 12 values +SF1 = xtgeo.RegularSurface( + xori=0, + yori=0, + ncol=3, + nrow=3, + xinc=3.0, + yinc=3.0, + values=np.array([10, 10, 10, 20, 0, 10, 10, 10, 1]), +) + +SF2 = xtgeo.RegularSurface( + xori=0, + yori=0, + ncol=3, + nrow=3, + xinc=3.0, + yinc=3.0, + values=np.array([100, 200, 300, 400, 0, 500, 600, 700, 800]), +) + +SMALL_POLY = [ + (2.9, 2.9, 0.0, 0), + (4.9, 2.9, 0.0, 0), + (6.1, 6.1, 0.0, 0), + (2.9, 4.0, 0.0, 0), + (2.9, 2.9, 0.0, 0), +] + +SMALL_POLY_OVERLAP = [ + (3.0, 0.0, 0.0, 2), + (6.1, -1.0, 0.0, 2), + (6.1, 3.2, 0.0, 2), + (2.0, 4.0, 0.0, 2), + (3.0, 0.0, 0.0, 2), +] + + +@pytest.fixture(name="reekset") +def fixture_reekset(): + """Read a test set from disk.""" + srf = xtgeo.surface_from_file(SURF1) + pol = xtgeo.polygons_from_file(POLY1) + return srf, pol + + +@pytest.mark.parametrize( + "oper, value, inside, expected", + [ + # O O O I I O I I I + # orig [10, 10, 10, 20, 0, 10, 10, 10, 1] + ("add", 1, True, ma.array([10, 10, 10, 21, 1, 10, 11, 11, 2])), + ("add", 1, False, ma.array([11, 11, 11, 20, 0, 11, 10, 10, 1])), + ("add", SF2, True, ma.array([10, 10, 10, 420, 0, 10, 610, 710, 801])), + ("add", SF2, False, ma.array([110, 210, 310, 20, 0, 510, 10, 10, 1])), + ("sub", 1, True, ma.array([10, 10, 10, 19, -1, 10, 9, 9, 0])), + ("sub", 1, False, ma.array([9, 9, 9, 20, 0, 9, 10, 10, 1])), + ("sub", SF2, True, ma.array([10, 10, 10, -380, 0, 10, -590, -690, -799])), + ("sub", SF2, False, ma.array([-90, -190, -290, 20, 0, -490, 10, 10, 1])), + ("mul", 2, True, ma.array([10, 10, 10, 40, 0, 10, 20, 20, 2])), + ("mul", 2, False, ma.array([20, 20, 20, 20, 0, 20, 10, 10, 1])), + ("div", 2, True, ma.array([10, 10, 10, 10, 0, 10, 5, 5, 0.5])), + ("div", 2, False, ma.array([5, 5, 5, 20, 0, 5, 10, 10, 1])), + ( + "eli", + 0, + True, + ma.array([False, False, False, True, True, False, True, True, True]), + ), + ( + "eli", + 0, + False, + ma.array([True, True, True, False, False, True, False, False, False]), + ), + ], +) +def test_operations_inside_outside_polygon_minimal(oper, value, inside, expected): + """Do operations in relation to polygons, minimal example""" + + poly = xtgeo.Polygons(SMALL_POLY + SMALL_POLY_OVERLAP) + + for ver in (1, 2): + surf = SF1.copy() + surf.operation_polygons(poly, value, inside=inside, opname=oper, _version=ver) + if oper == "eli": + np.testing.assert_array_equal(surf.values.mask.ravel(), expected) + else: + np.testing.assert_array_equal(surf.values.ravel(), expected) + + +@pytest.mark.parametrize( + "oper, inside, expected", + [ + ("add", True, 31.2033), + ("add", False, 70.7966), + ("sub", True, -29.2033), + ("sub", False, -68.7966), + ("mul", True, 30.9013), + ("mul", False, 70.0988), + ("div", True, 0.7010), + ("div", False, 0.3090), + ("set", True, 30.9013), + ("set", False, 70.0987), + ("eli", True, 1.0), + ("eli", False, 1.0), + ], +) +def test_operations_inside_outside_polygon_generic(reekset, oper, inside, expected): + """Do operations in relation to polygons in a generic way""" logger.info("Simple case...") - surf = xtgeo.surface_from_file(SURF1) - assert surf.values.mean() == pytest.approx(1698.65, abs=0.01) - poly = xtgeo.polygons_from_file(POLY1) + _surf, poly = reekset + _surf.values = 1.0 - surf.operation_polygons(poly, 100.0) # add 100 inside poly - assert surf.values.mean() == pytest.approx(1728.85, abs=0.01) + for version in (1, 2): + surf = _surf.copy() + surf.operation_polygons( + poly, 100.0, inside=inside, opname=oper, _version=version + ) + assert surf.values.mean() == pytest.approx(expected, abs=0.001) def test_operations_inside_outside_polygon_shortforms(tmpdir):