From 8b7038ea2cdc92ff8b96c1714e69cd82e6b1da6a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20C=2E=20Riven=C3=A6s?= Date: Tue, 15 Aug 2023 08:05:15 +0200 Subject: [PATCH] ENH: improve speed for surf vs polygons operations Replace internal methods with matplotlib MPath for surface vs polygons operations. When many points in the polygons, the speed improvement is magnitude of orders! In contrast the points in polygons, the 'add_inside()' etc methods works as before, hence the ``_version`` key is marked intentionally as private --- docs/usextgeoroxar.rst | 37 ++++- src/xtgeo/surface/_regsurf_oper.py | 99 ++++++++++++- src/xtgeo/surface/regular_surface.py | 16 ++- src/xtgeo/xyz/_xyz.py | 38 +++-- .../test_regular_surface_in_other.py | 130 ++++++++++++++++-- 5 files changed, 285 insertions(+), 35 deletions(-) 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):