From 14a93c703b1f7fc3d4c1b1a6a4a0be7d488a1cee 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 | 95 ++++++++++++- src/xtgeo/surface/regular_surface.py | 40 +++--- .../test_regular_surface_in_other.py | 130 ++++++++++++++++-- 4 files changed, 275 insertions(+), 27 deletions(-) diff --git a/docs/usextgeoroxar.rst b/docs/usextgeoroxar.rst index ca8a88a3a..e1402de2b 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(): + """Main work, looping wells and make APS relevant logs""" + + 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..d3dc4c2e4 100644 --- a/src/xtgeo/surface/_regsurf_oper.py +++ b/src/xtgeo/surface/_regsurf_oper.py @@ -5,6 +5,7 @@ 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 +453,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 +471,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 +520,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 +541,92 @@ 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, opname="add", inside=True): + """Operations restricted to polygons, using matplotlib (much faster)""" + + proxy = _proxy_map_polygons(self, poly, inside=inside) + print("\n", proxy.values) + 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": + print("R1\n", result.values) + result.values = np.ma.where( + proxy.values == 1, result.values + value, result.values + ) + print("R2\n", 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..09c9b12f2 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=1): """A generic function for map operations inside or outside polygon(s). Args: @@ -2051,59 +2051,67 @@ 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): """Add a value (scalar or other map) inside polygons.""" - self.operation_polygons(poly, value, opname="add", inside=True) + self.operation_polygons(poly, value, opname="add", inside=True, _version=2) def add_outside(self, poly, value): """Add a value (scalar or other map) outside polygons.""" - self.operation_polygons(poly, value, opname="add", inside=False) + self.operation_polygons(poly, value, opname="add", inside=False, _version=2) def sub_inside(self, poly, value): """Subtract a value (scalar or other map) inside polygons.""" - self.operation_polygons(poly, value, opname="sub", inside=True) + self.operation_polygons(poly, value, opname="sub", inside=True, _version=2) def sub_outside(self, poly, value): """Subtract a value (scalar or other map) outside polygons.""" - self.operation_polygons(poly, value, opname="sub", inside=False) + self.operation_polygons(poly, value, opname="sub", inside=False, _version=2) def mul_inside(self, poly, value): """Multiply a value (scalar or other map) inside polygons.""" - self.operation_polygons(poly, value, opname="mul", inside=True) + self.operation_polygons(poly, value, opname="mul", inside=True, _version=2) def mul_outside(self, poly, value): """Multiply a value (scalar or other map) outside polygons.""" - self.operation_polygons(poly, value, opname="mul", inside=False) + self.operation_polygons(poly, value, opname="mul", inside=False, _version=2) def div_inside(self, poly, value): """Divide a value (scalar or other map) inside polygons.""" - self.operation_polygons(poly, value, opname="div", inside=True) + self.operation_polygons(poly, value, opname="div", inside=True, _version=2) def div_outside(self, poly, value): """Divide a value (scalar or other map) outside polygons.""" - self.operation_polygons(poly, value, opname="div", inside=False) + self.operation_polygons(poly, value, opname="div", inside=False, _version=2) def set_inside(self, poly, value): """Set a value (scalar or other map) inside polygons.""" - self.operation_polygons(poly, value, opname="set", inside=True) + self.operation_polygons(poly, value, opname="set", inside=True, _version=2) def set_outside(self, poly, value): """Set a value (scalar or other map) outside polygons.""" - self.operation_polygons(poly, value, opname="set", inside=False) + self.operation_polygons(poly, value, opname="set", inside=False, _version=2) def eli_inside(self, poly): """Eliminate current map values inside polygons.""" - self.operation_polygons(poly, 0, opname="eli", inside=True) + self.operation_polygons(poly, 0, opname="eli", inside=True, _version=2) def eli_outside(self, poly): """Eliminate current map values outside polygons.""" - self.operation_polygons(poly, 0, opname="eli", inside=False) + self.operation_polygons(poly, 0, opname="eli", inside=False, _version=2) # ================================================================================== # Operation with secondary map 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):