From 2369745837fecdf23268ca78183fb7ff20cc3c2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jan=20C=2E=20Riven=C3=A6s?= Date: Tue, 15 Aug 2023 07:58:25 +0200 Subject: [PATCH] ENH: improve points in polygons methods Improve both speed and logics. Speed by using MPath from matplotlib, and logics to adress the strange behaviour in previous versions, when e.g. ovelaping polygons areas gets cumulated values. --- src/xtgeo/xyz/_xyz.py | 397 ++++++++++++++++++++++++++--- src/xtgeo/xyz/_xyz_oper.py | 89 ++++++- tests/test_xyz/test_points_poly.py | 339 ++++++++++++++++++------ 3 files changed, 712 insertions(+), 113 deletions(-) diff --git a/src/xtgeo/xyz/_xyz.py b/src/xtgeo/xyz/_xyz.py index 8183f3a6a..7b9e33293 100644 --- a/src/xtgeo/xyz/_xyz.py +++ b/src/xtgeo/xyz/_xyz.py @@ -1,6 +1,8 @@ # -*- coding: utf-8 -*- """XTGeo XYZ module (abstract base class)""" from abc import ABC, abstractmethod +from typing import List, TypeVar, Union +from warnings import warn import numpy as np import pandas as pd @@ -198,7 +200,7 @@ def protected_columns(self): Columns not deleted by :meth:`delete_columns`, for instance the coordinate columns. """ - return [self.xname, self.yname, self.zname, self.pname] + return [self.xname, self.yname, self.zname] def geometry_columns(self): """ @@ -271,7 +273,36 @@ def get_boundary(self): return (xmin, xmax, ymin, ymax, zmin, zmax) - def operation_polygons(self, poly, value, opname="add", inside=True): + def mark_in_polygons( + self, + poly: Union["Polygons", List["Polygons"]], + name: str = "pstatus", + inside_value: int = 1, + outside_value: int = 0, + ): + """Add a column that assign values if points are inside or outside polygons. + + This is a generic function that adds a column in the points dataframe with + a flag for values being inside or outside polygons in a Polygons instance. + + Args: + poly: One single xtgeo Polgons instance, or a list of Polygons instances. + name: Name of column that flags inside or outside status + inside_value: Flag value for being inside polygons + outside_value: Flag value for being outside polygons + + ..versionadded:: 3.2 + """ + _xyz_oper.mark_in_polygons_mpl(self, poly, name, inside_value, outside_value) + + def operation_polygons( + self, + poly: Union["Polygons", List["Polygons"]], + value: float, + opname: str = "add", + inside: bool = True, + version: int = 1, + ): """A generic function for operations restricted to inside or outside polygon(s). The operations are performed on the Z values, while the 'inside' or 'outside' @@ -291,103 +322,393 @@ def operation_polygons(self, poly, value, opname="add", inside=True): * ``eli``: eliminate; here value is not applied Args: - poly (Polygons): A XTGeo Polygons instance - value(float): Value to add, subtract etc - opname (str): Name of operation... 'add', 'sub', etc - inside (bool): If True do operation inside polygons; else outside. Note + poly: A single Polygons instance or a list of Polygons instances. + The list option is only allowed when version = 2 + value: Value to add, subtract etc + opname: Name of operation... 'add', 'sub', etc + inside: If True do operation inside polygons; else outside. Note that boundary is treated as 'inside' + version: The algorithm version, see notes below. Although version 1 + is default, version 2 is recommended as it is much faster and works + intuitively when have multiple polygons and/or using the + `is_inside=False` (i.e. outside) Note: - This function works only intuitively when using one single polygon - in the ``poly`` instance. When having several polygons the + ``version=1``: This function works only intuitively when using one single + polygon in the ``poly`` instance. When having several polygons the operation is done sequentially per polygon which may lead to surprising results. For instance, using "add inside" into two overlapping polygons, the addition will be doubled in the - overlapping part. Similarly using e.g. "eli, outside" will completely + 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 + much faster and intuitively for both inside and outside, overlapping and + multiple polygons within a Polygons instance. + + .. versionchanged:: 3.2 Add ``version`` option which defaults to 1. + Also allow that ``poly`` option can be a list of Polygons + when version is 2. + """ - _xyz_oper.operation_polygons(self, poly, value, opname=opname, inside=inside) + if version == 2: + _xyz_oper.operation_polygons_v2( + self, poly, value, opname=opname, inside=inside + ) + else: + _xyz_oper.operation_polygons_v1( + self, poly, value, opname=opname, inside=inside + ) + if version == 0: + # using version 0 INTERNALLY to mark that "add_inside" etc has been + # applied, and now raise a generic deprecation warning: + itxt = "inside" if inside else "outside" + warn( + f"You are using the method '{opname}_{itxt}()'; this will " + "be deprecated in future versions. Consider using " + f"'{opname}_{itxt}_polygons()' instead which is both faster " + "and works intuitively when several and/or overlapping polygons", + DeprecationWarning, + ) def add_inside(self, poly, value): - """Add a value (scalar) to points inside polygons. + """Add a value (scalar) to points inside polygons (old behaviour). + + Args: + poly: A xtgeo Polygons instance + value: Value to add to Z values inside polygons. - See notes under :meth:`operation_polygons`. + See notes under :meth:`operation_polygons()` and consider instead + :meth:`add_inside polygons()`. """ - self.operation_polygons(poly, value, opname="add", inside=True) + self.operation_polygons(poly, value, opname="add", inside=True, version=0) + + def add_inside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Add a value (scalar) to points inside polygons (new behaviour). + + This is an improved implementation than :meth:`add_inside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to add to Z values inside polygons. + + .. versionadded:: 3.2 + """ + self.operation_polygons(poly, value, opname="add", inside=True, version=2) def add_outside(self, poly, value): - """Add a value (scalar) to points outside polygons. + """Add a value (scalar) to points outside polygons (old behaviour). + + Args: + poly: A xtgeo Polygons instance + value: Value to add to Z values outside polygons. + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`add_outside polygons()`. + """ + self.operation_polygons(poly, value, opname="add", inside=False, version=0) + + def add_outside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Add a value (scalar) to points outside polygons (new behaviour). + + This is an improved implementation than :meth:`add_outside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to add to Z values outside polygons. - See notes under :meth:`operation_polygons`. + .. versionadded:: 3.2 """ - 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) to points inside polygons. - See notes under :meth:`operation_polygons`. + Args: + poly: A xtgeo Polygons instance + value: Value to subtract to Z values inside polygons. + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`sub_inside polygons()`. + """ + self.operation_polygons(poly, value, opname="sub", inside=True, version=1) + + def sub_inside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Subtract a value (scalar) for points inside polygons (new behaviour). + + This is an improved implementation than :meth:`sub_inside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to subtract to Z values inside polygons. + + .. versionadded:: 3.2 """ - 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) to points outside polygons. - See notes under :meth:`operation_polygons`. + Args: + poly: A xtgeo Polygons instance + value: Value to subtract to Z values outside polygons. + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`sub_outside polygons()`. """ - self.operation_polygons(poly, value, opname="sub", inside=False) + self.operation_polygons(poly, value, opname="sub", inside=False, version=0) + + def sub_outside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Subtract a value (scalar) for points outside polygons (new behaviour). + + This is an improved implementation than :meth:`sub_outside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to subtract to Z values outside polygons. + + .. versionadded:: 3.2 + """ + self.operation_polygons(poly, value, opname="sub", inside=False, version=2) def mul_inside(self, poly, value): """Multiply a value (scalar) to points inside polygons. - See notes under :meth:`operation_polygons`. + Args: + poly: A xtgeo Polygons instance + value: Value to multiply to Z values inside polygons. + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`mul_inside polygons()`. + """ + self.operation_polygons(poly, value, opname="mul", inside=True, version=0) + + def mul_inside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Multiply a value (scalar) for points inside polygons (new behaviour). + + This is an improved implementation than :meth:`mul_inside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to multiply to Z values inside polygons. + + .. versionadded:: 3.2 """ - 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) to points outside polygons. - See notes under :meth:`operation_polygons`. + Args: + poly: A xtgeo Polygons instance + value: Value to multiply to Z values outside polygons. + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`mul_outside polygons()`. + """ + self.operation_polygons(poly, value, opname="mul", inside=False, version=0) + + def mul_outside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Multiply a value (scalar) for points outside polygons (new behaviour). + + This is an improved implementation than :meth:`mul_outside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to multiply to Z values outside polygons. + + .. versionadded:: 3.2 """ - 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) to points inside polygons. - See notes under :meth:`operation_polygons`. + Args: + poly: A xtgeo Polygons instance + value: Value to divide Z values inside polygons. + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`div_inside polygons()`. """ - self.operation_polygons(poly, value, opname="div", inside=True) + self.operation_polygons(poly, value, opname="div", inside=True, version=0) + + def div_inside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Divide a value (scalar) for points inside polygons (new behaviour). + + This is an improved implementation than :meth:`div_inside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to divide to Z values inside polygons. + + .. versionadded:: 3.2 + """ + self.operation_polygons(poly, value, opname="div", inside=True, version=2) def div_outside(self, poly, value): """Divide a value (scalar) outside polygons (value 0.0 will give result 0). - See notes under :meth:`operation_polygons`. + Args: + poly: A xtgeo Polygons instance + value: Value to divide Z values outside polygons. + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`div_outside polygons()`. """ - self.operation_polygons(poly, value, opname="div", inside=False) + self.operation_polygons(poly, value, opname="div", inside=False, version=0) + + def div_outside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Divide a value (scalar) for points outside polygons (new behaviour). + + Note if input value is 0.0 (division on zero), the result will be 0.0. + + This is an improved implementation than :meth:`div_outside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to divide to Z values outside polygons. + + .. versionadded:: 3.2 + """ + self.operation_polygons(poly, value, opname="div", inside=False, version=2) def set_inside(self, poly, value): """Set a value (scalar) to points inside polygons. - See notes under :meth:`operation_polygons`. + Args: + poly: A xtgeo Polygons instance + value: Value to set Z values inside polygons. + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`set_inside polygons()`. + """ + self.operation_polygons(poly, value, opname="set", inside=True, version=0) + + def set_inside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Set a value (scalar) for points inside polygons (new behaviour). + + This is an improved implementation than :meth:`set_inside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to set as Z values inside polygons. + + .. versionadded:: 3.2 """ - 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) to points outside polygons. - See notes under :meth:`operation_polygons`. + Args: + poly: A xtgeo Polygons instance + value: Value to set Z values outside polygons. + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`set_outside polygons()`. """ - self.operation_polygons(poly, value, opname="set", inside=False) + self.operation_polygons(poly, value, opname="set", inside=False, version=0) + + def set_outside_polygons( + self, poly: Union["Polygons", List["Polygons"]], value: float + ): + """Set a value (scalar) for points outside polygons (new behaviour). + + This is an improved implementation than :meth:`set_outside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + value: Value to set as Z values inside polygons. + + .. versionadded:: 3.2 + """ + self.operation_polygons(poly, value, opname="set", inside=False, version=2) def eli_inside(self, poly): - """Eliminate current points inside polygons. + """Eliminate current points inside polygons (old implentation). - See notes under :meth:`operation_polygons`. + Args: + poly: A xtgeo Polygons instance + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`eli_inside polygons()`. + """ + self.operation_polygons(poly, 0, opname="eli", inside=True, version=0) + + def eli_inside_polygons(self, poly: Union["Polygons", List["Polygons"]]): + """Remove points inside polygons. + + This is an improved implementation than :meth:`eli_inside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances + + .. versionadded:: 3.2 """ - 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 points outside polygons. + """Eliminate current points outside polygons (old implentation). + + Args: + poly: A xtgeo Polygons instance + + See notes under :meth:`operation_polygons()` and consider instead + :meth:`eli_outside polygons()`. + """ + self.operation_polygons(poly, 0, opname="eli", inside=False, version=0) + + def eli_outside_polygons(self, poly: Union["Polygons", List["Polygons"]]): + """Remove points outside polygons. + + This is an improved implementation than :meth:`eli_outside()`, and is now the + recommended method, as it is both faster and works similar for all single + and overlapping sub-polygons within one or more Polygons instances. + + Args: + poly: A xtgeo Polygons instance, or a list of xtgeo Polygons instances - See notes under :meth:`operation_polygons`. + .. versionadded:: 3.2 """ - self.operation_polygons(poly, 0, opname="eli", inside=False) + self.operation_polygons(poly, 0, opname="eli", inside=False, version=2) diff --git a/src/xtgeo/xyz/_xyz_oper.py b/src/xtgeo/xyz/_xyz_oper.py index 10b51ab4d..076e53fdc 100644 --- a/src/xtgeo/xyz/_xyz_oper.py +++ b/src/xtgeo/xyz/_xyz_oper.py @@ -5,6 +5,7 @@ import numpy as np import pandas as pd import shapely.geometry as sg +from matplotlib.path import Path as MPath from scipy.interpolate import UnivariateSpline, interp1d import xtgeo @@ -19,7 +20,37 @@ # pylint: disable=protected-access -def operation_polygons(self, poly, value, opname="add", inside=True, where=True): +def mark_in_polygons_mpl(self, poly, name, inside_value, outside_value): + """Make column to mark if XYZ df is inside/outside polygons, using matplotlib.""" + points = np.array( + [self.dataframe[self.xname].values, self.dataframe[self.yname].values] + ).T + + if name in self.protected_columns(): + raise ValueError("The proposed name: {name}, is protected and cannot be used") + + # 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") + + self.dataframe[name] = outside_value + + 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) + self.dataframe.loc[is_inside, name] = inside_value + + +def operation_polygons_v1(self, poly, value, opname="add", inside=True, where=True): """ Operations re restricted to closed polygons, for points or polyline points. @@ -30,8 +61,10 @@ def operation_polygons(self, poly, value, opname="add", inside=True, where=True) will be the intersection. The "where" filter is reserved for future use. + + This is now the legacy version, which will be deprecated later. """ - logger.warning("Where is not imeplented: %s", where) + logger.warning("Where is not implemented: %s", where) oper = {"set": 1, "add": 2, "sub": 3, "mul": 4, "div": 5, "eli": 11} @@ -41,7 +74,7 @@ def operation_polygons(self, poly, value, opname="add", inside=True, where=True) logger.info("Operations of points inside polygon(s)...") if not isinstance(poly, xtgeo.xyz.Polygons): - raise ValueError("The poly input is not a Polygons instance") + raise ValueError("The poly input is not a single Polygons instance") idgroups = poly.dataframe.groupby(poly.pname) @@ -80,6 +113,56 @@ def operation_polygons(self, poly, value, opname="add", inside=True, where=True) logger.info("Operations of points inside polygon(s)... done") +def operation_polygons_v2(self, poly, value, opname="add", inside=True, where=True): + """ + Operations re restricted to closed polygons, for points or polyline points. + + This (from v1) is NOT implemented here: "If value is not float but 'poly', then the + avg of each polygon Z value will be used instead." + + NOTE! Both 'inside' and "outside" several polygons will become a union. Hence v2 + it will TREAT OUTSIDE POINTS DIFFERENT than v1! + + The "where" filter is reserved for future use. + + The version v2 uses mark_in_polygons_mpl(), and should be much faster. + """ + logger.info("Operations of points inside polygon(s), version 2") + + allowed_opname = ("set", "add", "sub", "mul", "div", "eli") + + logger.warning("Where is not implemented: %s", where) + if not isinstance(poly, xtgeo.xyz.Polygons): + raise ValueError("The poly input is not a Polygons instance") + + tmp = self.copy() + tmpdf = tmp.dataframe + ivalue = 1 if inside else 0 + ovalue = 0 if inside else 1 + mark_in_polygons_mpl(tmp, poly, "_TMP", inside_value=ivalue, outside_value=ovalue) + + if opname == "add": + self.dataframe[self.zname][tmpdf._TMP == 1] += value + elif opname == "sub": + self.dataframe[self.zname][tmpdf._TMP == 1] -= value + elif opname == "mul": + self.dataframe[self.zname][tmpdf._TMP == 1] *= value + elif opname == "div": + if value != 0.0: + self.dataframe[self.zname][tmpdf._TMP == 1] /= value + else: + self.dataframe[self.zname][tmpdf._TMP == 1] = 0.0 + elif opname == "set": + self.dataframe[self.zname][tmpdf._TMP == 1] = value + elif opname == "eli": + self.dataframe = self.dataframe[tmpdf._TMP == 0] + self.dataframe.reset_index(inplace=True, drop=True) + else: + raise KeyError(f"The opname={opname} is not one of {allowed_opname}") + + logger.info("Operations of points inside polygon(s)... done") + + def rescale_polygons(self, distance=10, addlen=False, kind="simple", mode2d=False): """Rescale (resample) a polygons segment Default settings will make it backwards compatible with 2.0 diff --git a/tests/test_xyz/test_points_poly.py b/tests/test_xyz/test_points_poly.py index a02694a7a..4cf396126 100644 --- a/tests/test_xyz/test_points_poly.py +++ b/tests/test_xyz/test_points_poly.py @@ -1,6 +1,7 @@ import pathlib import sys +import numpy as np import pytest import xtgeo @@ -50,44 +51,101 @@ ] -def test_points_in_polygon(testpath): - """Import XYZ points and do operations if inside or outside""" - +@pytest.fixture(name="reekset") +def fixture_reekset(testpath): + """Read a test set from disk.""" poi = xtgeo.points_from_file(testpath / POINTSET2) pol = xtgeo.polygons_from_file(testpath / POLSET2) + return poi, pol + + +def test_mark_points_polygon(reekset): + """Import XYZ points and get marks of points inside or outside""" + poi, pol = reekset + poi2 = poi.copy() + # remove points in polygons + poi2.mark_in_polygons(pol, name="_ptest", inside_value=1, outside_value=0) + assert poi2.dataframe["_ptest"].values.sum() == 11 + + # a.a, but allow multiple polygons (reuse same here) + poi3 = poi.copy() + poi3.mark_in_polygons([pol, pol], name="_ptest", inside_value=1, outside_value=0) + assert poi3.dataframe["_ptest"].values.sum() == 11 + + +def test_mark_points_polygon_protected(reekset): + """Raise a ValuError if a protected name is attempted""" + poi, pol = reekset + poi2 = poi.copy() + with pytest.raises(ValueError, match="is protected and cannot be used"): + poi2.mark_in_polygons(pol, name=poi.xname, inside_value=1, outside_value=0) + + +def test_points_in_polygon(reekset): + """Import XYZ points and do operations if inside or outside""" + + _poi, pol = reekset + poi = _poi.copy() assert poi.nrow == 30 # remove points in polygon - poi.operation_polygons(pol, 0, opname="eli") + poi.operation_polygons(pol, 0, opname="eli", version=2) assert poi.nrow == 19 - poi = xtgeo.points_from_file(testpath / POINTSET2) # remove points outside polygon - poi.operation_polygons(pol, 0, opname="eli", inside=False) - assert poi.nrow == 1 + poi.operation_polygons(pol, 0, opname="eli", inside=False, version=1) + assert poi.nrow == 0 + + # new behaviour with version = 2: + poi = _poi.copy() # refresh + poi.operation_polygons(pol, 0, opname="eli", inside=False, version=2) + print(poi.dataframe) + assert poi.nrow == 11 + + +def test_add_inside_old_new_behaviour(reekset): + """Add inside using old and new class behaviour""" + _poi, pol = reekset + poi1 = _poi.copy() + poi1.dataframe.Z_TVDSS = 0.0 + poi2 = _poi.copy() + poi2.dataframe.Z_TVDSS = 0.0 + + poi1.add_inside(pol, 1) + print(poi1.dataframe) + zvec = poi1.dataframe["Z_TVDSS"].values + assert 2.0 in zvec.tolist() # will be doubling where polygons overlap! + zvec = zvec[zvec < 1] + assert zvec.size == 19 + + poi2.add_inside_polygons(pol, 1) + zvec = poi2.dataframe["Z_TVDSS"].values + assert 2.0 not in zvec.tolist() # here NO doubling where polygons overlap! + zvec = zvec[zvec < 1] + print(poi2.dataframe) + assert zvec.size == 19 def test_check_single_point_in_polygon(): - pol = Polygons(SMALL_POLY_INNER) - poi = Points([(4.0, 4.0, 4.0)]) - poi.operation_polygons(pol, value=1, opname="eli", inside=False) - assert len(poi.dataframe) == 1 - - poi.operation_polygons(pol, value=1, opname="eli", inside=True) - assert len(poi.dataframe) == 0 + for ver in (1, 2): + pol = Polygons(SMALL_POLY_INNER) + poi = Points([(4.0, 4.0, 4.0)]) + poi.operation_polygons(pol, value=1, opname="eli", inside=False, version=ver) + assert len(poi.dataframe) == 1 def test_check_multi_point_single_polygon(): - pol = Polygons(SMALL_POLY_INNER) - poi = Points([(6.0, 4.0, 4.0), (4.0, 4.0, 4.0), (4.0, 6.0, 4.0)]) - assert len(poi.dataframe) == 3 + for ver in (1, 2): + pol = Polygons(SMALL_POLY_INNER) + poi = Points([(6.0, 4.0, 4.0), (4.0, 4.0, 4.0), (4.0, 6.0, 4.0)]) + assert len(poi.dataframe) == 3 - poi.operation_polygons(pol, value=1, opname="eli", inside=False) - assert len(poi.dataframe) == 1 + poi.operation_polygons(pol, value=1, opname="eli", inside=False, version=ver) + assert len(poi.dataframe) == 1 - poi.operation_polygons(pol, value=1, opname="eli", inside=True) - assert len(poi.dataframe) == 0 + poi.operation_polygons(pol, value=1, opname="eli", inside=True, version=ver) + assert len(poi.dataframe) == 0 def test_check_multi_point_single_polygon_zdir(): @@ -103,15 +161,20 @@ def test_check_multi_point_single_polygon_zdir(): def test_check_multi_point_multi_polyon_inside_op(): pol = Polygons(SMALL_POLY_INNER + LARGE_POLY_SHIFTED) # Two points in small cube, one in large cube, one outside - poi = Points([(4.0, 4.0, 0.0), (4.5, 4.0, 0.0), (7.0, 7.0, 0.0), (20.0, 5.0, 0.0)]) - assert len(poi.dataframe) == 4 + for ver in (1, 2): + poi = Points( + [(4.0, 4.0, 0.0), (4.5, 4.0, 0.0), (7.0, 7.0, 0.0), (20.0, 5.0, 0.0)] + ) + assert len(poi.dataframe) == 4 - poi.operation_polygons(pol, value=1, opname="eli", inside=True) - assert len(poi.dataframe) == 1 + poi.operation_polygons(pol, value=1, opname="eli", inside=True, version=ver) + assert len(poi.dataframe) == 1 - poi = Points([(4.0, 4.0, 0.0), (4.5, 4.0, 0.0), (7.0, 7.0, 0.0), (20.0, 5.0, 0.0)]) - poi.operation_polygons(pol, value=1, opname="eli", inside=False) - assert len(poi.dataframe) == 0 + poi = Points( + [(4.0, 4.0, 0.0), (4.5, 4.0, 0.0), (7.0, 7.0, 0.0), (20.0, 5.0, 0.0)] + ) + poi.operation_polygons(pol, value=1, opname="eli", inside=False, version=ver) + assert len(poi.dataframe) == 0 if ver == 1 else 3 # dependent on version! def test_check_multi_point_multi_polyon_outside_op(): @@ -120,72 +183,94 @@ def test_check_multi_point_multi_polyon_outside_op(): poi = Points([(4.0, 4.0, 0.0), (4.5, 4.0, 0.0), (7.0, 7.0, 0.0), (20.0, 5.0, 0.0)]) assert len(poi.dataframe) == 4 - # Note the operation will loop over the polygons, and hence remove the points - # in the small polygon when considering the large polygon, and vice versa - poi.operation_polygons(pol, value=1, opname="eli", inside=False) + # Note the operation will loop over the polygons if version 1, and hence remove the + # points in the small polygon when considering the large polygon, and vice versa + poi.operation_polygons(pol, value=1, opname="eli", inside=False, version=1) assert len(poi.dataframe) == 0 + # Fixed in version 2 + poi = Points([(4.0, 4.0, 0.0), (4.5, 4.0, 0.0), (7.0, 7.0, 0.0), (20.0, 5.0, 0.0)]) + poi.operation_polygons(pol, value=1, opname="eli", inside=False, version=2) + assert len(poi.dataframe) == 3 + def test_check_single_polygon_in_single_polygon(): - inner_pol = Polygons(SMALL_POLY_INNER) - outer_pol = Polygons(SMALL_POLY_OUTER) + for ver in (1, 2): + inner_pol = Polygons(SMALL_POLY_INNER) + outer_pol = Polygons(SMALL_POLY_OUTER) - # Do not delete inner_pol when specified to delete polygons inside inner polygon - inner_pol.operation_polygons(outer_pol, value=1, opname="eli", inside=False) - assert len(inner_pol.dataframe) == 5 + # Do not delete inner_pol when specified to delete polygons outside outer polygon + inner_pol.operation_polygons( + outer_pol, value=1, opname="eli", inside=False, version=ver + ) + assert len(inner_pol.dataframe) == 5 - # Do not delete outer_pol when specified to delete polygons outside outer polygon - outer_pol.operation_polygons(inner_pol, value=1, opname="eli", inside=True) - assert len(outer_pol.dataframe) == 5 + # Do not delete outer_pol when specified to delete polygons inside inner polygon + outer_pol.operation_polygons( + inner_pol, value=1, opname="eli", inside=True, version=ver + ) + assert len(outer_pol.dataframe) == 5 - inner_pol.operation_polygons(outer_pol, value=1, opname="eli", inside=True) - assert len(inner_pol.dataframe) == 0 + inner_pol.operation_polygons( + outer_pol, value=1, opname="eli", inside=True, version=ver + ) + assert len(inner_pol.dataframe) == 0 - inner_pol = Polygons(SMALL_POLY_INNER) + inner_pol = Polygons(SMALL_POLY_INNER) - outer_pol.operation_polygons(inner_pol, value=1, opname="eli", inside=False) - assert len(outer_pol.dataframe) == 0 + outer_pol.operation_polygons( + inner_pol, value=1, opname="eli", inside=False, version=ver + ) + assert len(outer_pol.dataframe) == 0 def test_check_multi_polygon_in_single_polygon(): - inner_pols = Polygons(SMALL_POLY_INNER + SMALL_POLY_OVERLAP_INNER) - outer_pol = Polygons(SMALL_POLY_OUTER) + for ver in (1, 2): + inner_pols = Polygons(SMALL_POLY_INNER + SMALL_POLY_OVERLAP_INNER) + outer_pol = Polygons(SMALL_POLY_OUTER) - inner_pols.operation_polygons(outer_pol, value=1, opname="eli", inside=True) - assert len(inner_pols.dataframe) == 0 + inner_pols.operation_polygons( + outer_pol, value=1, opname="eli", inside=True, version=ver + ) + assert len(inner_pols.dataframe) == 0 def test_operation_inclusive_polygon(): - pol = Polygons(SMALL_POLY_INNER) - # We place a point on the edge of a polygon - poi = Points([(4.0, 4.0, 0.0)]) - poi.operation_polygons(pol, value=1, opname="eli", inside=True) - assert len(poi.dataframe) == 0 + for ver in (1, 2): + pol = Polygons(SMALL_POLY_INNER) + # We place a point on the edge of a polygon + poi = Points([(4.0, 4.0, 0.0)]) + poi.operation_polygons(pol, value=1, opname="eli", inside=True, version=ver) + assert len(poi.dataframe) == 0 - # We place a point on a corner of a polygon - poi = Points([(3.0, 3.0, 0.0)]) - poi.operation_polygons(pol, value=1, opname="eli", inside=True) - assert len(poi.dataframe) == 0 + # We place a point on a corner of a polygon + poi = Points([(3.0, 3.0, 0.0)]) + poi.operation_polygons(pol, value=1, opname="eli", inside=True) + assert len(poi.dataframe) == 0 def test_polygons_overlap(): - pol = Polygons(SMALL_POLY_INNER + SMALL_POLY_OVERLAP_INNER) - # The Four points are placed: within first poly, within the overlap, within the - # second poly, outside both poly - poi = Points([(3.5, 3.5, 0.0), (4.5, 4.5, 0.0), (5.5, 5.5, 0.0), (6.5, 6.5, 0.0)]) - poi.operation_polygons(pol, value=1, opname="eli", inside=True) - assert len(poi.dataframe) == 1 + for ver in (1, 2): + pol = Polygons(SMALL_POLY_INNER + SMALL_POLY_OVERLAP_INNER) + # The Four points are placed: within first poly, within the overlap, within the + # second poly, outside both poly + poi = Points( + [(3.5, 3.5, 0.0), (4.5, 4.5, 0.0), (5.5, 5.5, 0.0), (6.5, 6.5, 0.0)] + ) + poi.operation_polygons(pol, value=1, opname="eli", inside=True, version=ver) + assert len(poi.dataframe) == 1 @pytest.mark.parametrize( "oper, expected", [("add", 12), ("sub", 8), ("div", 5), ("mul", 20), ("set", 2)] ) def test_oper_single_point_in_polygon(oper, expected): - pol = Polygons(SMALL_POLY_INNER) - poi = Points([(4.0, 4.0, 10.0)]) - # operators work on z-values - poi.operation_polygons(pol, value=2, opname=oper, inside=True) - assert poi.dataframe[poi.zname].values[0] == expected + for ver in (1, 2): + pol = Polygons(SMALL_POLY_INNER) + poi = Points([(4.0, 4.0, 10.0)]) + # operators work on z-values + poi.operation_polygons(pol, value=2, opname=oper, inside=True, version=ver) + assert poi.dataframe[poi.zname].values[0] == expected @pytest.mark.parametrize( @@ -198,7 +283,7 @@ def test_oper_single_point_in_polygon(oper, expected): ("set", [2, 2, 2, 10]), ], ) -def test_oper_points_outside_overlapping_polygon(oper, expected): +def test_oper_points_outside_overlapping_polygon_v1(oper, expected): pol = Polygons(SMALL_POLY_INNER + SMALL_POLY_OVERLAP_INNER) # The Four points are placed: within first poly, within the overlap, within the # second poly, outside both poly @@ -212,7 +297,36 @@ def test_oper_points_outside_overlapping_polygon(oper, expected): ) # Operations are performed n times, where n reflects the number of polygons a # point is within - poi.operation_polygons(pol, value=2, opname=oper, inside=True) + poi.operation_polygons(pol, value=2, opname=oper, inside=True, version=1) + assert list(poi.dataframe[poi.zname].values) == expected + + +@pytest.mark.parametrize( + "oper, expected", + [ + ("add", [12, 12, 12, 10]), + ("sub", [8, 8, 8, 10]), + ("div", [5, 5, 5, 10]), + ("mul", [20, 20, 20, 10]), + ("set", [2, 2, 2, 10]), + ], +) +def test_oper_points_outside_overlapping_polygon_v2(oper, expected): + # different expected values her for version 2, compared with version 1 + pol = Polygons(SMALL_POLY_INNER + SMALL_POLY_OVERLAP_INNER) + # The Four points are placed: within first poly, within the overlap, within the + # second poly, outside both poly + poi = Points( + [ + (3.5, 3.5, 10.0), + (4.5, 4.5, 10.0), + (5.5, 5.5, 10.0), + (6.5, 6.5, 10.0), + ] + ) + # Operations are performed n times, where n reflects the number of polygons a + # point is within + poi.operation_polygons(pol, value=2, opname=oper, inside=True, version=2) assert list(poi.dataframe[poi.zname].values) == expected @@ -226,7 +340,7 @@ def test_oper_points_outside_overlapping_polygon(oper, expected): ("set", [2, 10, 2, 2]), ], ) -def test_oper_points_inside_overlapping_polygon(oper, expected): +def test_oper_points_inside_overlapping_polygon_v1(oper, expected): pol = Polygons(SMALL_POLY_INNER + SMALL_POLY_OVERLAP_INNER) # The Four points are placed: within first poly, within the overlap, within the # second poly, outside both poly @@ -235,10 +349,91 @@ def test_oper_points_inside_overlapping_polygon(oper, expected): ) # Operations are performed n times, where n reflects the number of polygons a # point is outside - poi.operation_polygons(pol, value=2, opname=oper, inside=False) + poi.operation_polygons(pol, value=2, opname=oper, inside=False, version=1) assert list(poi.dataframe[poi.zname].values) == expected +@pytest.mark.parametrize( + "oper, expected", + [ + ("add", [10, 10, 10, 12]), + ("sub", [10, 10, 10, 8]), + ("div", [10, 10, 10, 5]), + ("mul", [10, 10, 10, 20]), + ("set", [10, 10, 10, 2]), + ], +) +def test_oper_points_inside_overlapping_polygon_v2(oper, expected): + pol = Polygons(SMALL_POLY_INNER + SMALL_POLY_OVERLAP_INNER) + # The Four points are placed: within first poly, within the overlap, within the + # second poly, outside both poly + poi = Points( + [(3.5, 3.5, 10.0), (4.5, 4.5, 10.0), (5.5, 5.5, 10.0), (6.5, 6.5, 10.0)] + ) + # Operations are performed n times, where n reflects the number of polygons a + # point is outside + poi.operation_polygons(pol, value=2, opname=oper, inside=False, version=2) + assert list(poi.dataframe[poi.zname].values) == expected + + +def test_add_inside_polygons_etc(): + pol = Polygons(SMALL_POLY_INNER + SMALL_POLY_OVERLAP_INNER) + # The Four points are placed: within first poly, within the overlap, within the + # second poly, outside both poly + _poi = Points( + [(3.5, 3.5, 10.0), (4.5, 4.5, 10.0), (5.5, 5.5, 10.0), (6.5, 6.5, 10.0)] + ) + poi = _poi.copy() + poi.add_inside_polygons(pol, 10.0) + assert list(poi.dataframe[poi.zname].values) == [20.0, 20.0, 20.0, 10.0] + + poi = _poi.copy() + poi.add_outside_polygons(pol, 10.0) + assert list(poi.dataframe[poi.zname].values) == [10.0, 10.0, 10.0, 20.0] + + poi = _poi.copy() + poi.sub_inside_polygons(pol, 10.0) + assert list(poi.dataframe[poi.zname].values) == [0.0, 0.0, 0.0, 10.0] + + poi = _poi.copy() + poi.sub_outside_polygons(pol, 10.0) + assert list(poi.dataframe[poi.zname].values) == [10.0, 10.0, 10.0, 0.0] + + poi = _poi.copy() + poi.mul_inside_polygons(pol, 10.0) + assert list(poi.dataframe[poi.zname].values) == [100.0, 100.0, 100.0, 10.0] + + poi = _poi.copy() + poi.mul_outside_polygons(pol, 10.0) + assert list(poi.dataframe[poi.zname].values) == [10.0, 10.0, 10.0, 100.0] + + poi = _poi.copy() + poi.div_inside_polygons(pol, 10.0) + assert list(poi.dataframe[poi.zname].values) == [1.0, 1.0, 1.0, 10.0] + + poi = _poi.copy() + poi.div_outside_polygons(pol, 10.0) + assert list(poi.dataframe[poi.zname].values) == [10.0, 10.0, 10.0, 1.0] + + # zero division + poi = _poi.copy() + poi.div_outside_polygons(pol, 0.0) + assert list(poi.dataframe[poi.zname].values) == [10.0, 10.0, 10.0, 0.0] + + # close to zero division + poi = _poi.copy() + poi.div_outside_polygons(pol, 1e-10) + assert list(poi.dataframe[poi.zname].values) == [10.0, 10.0, 10.0, 1e11] + + poi = _poi.copy() + poi.eli_inside_polygons(pol) + assert list(poi.dataframe[poi.zname].values) == [10.0] + + poi = _poi.copy() + poi.eli_outside_polygons(pol) + assert list(poi.dataframe[poi.zname].values) == [10.0, 10.0, 10.0] + + @pytest.mark.skipif(sys.version_info < (3, 8), reason="Different order in python 3.7") def test_boundary_from_points_simple(): """Test deriving a boundary around points (classmethod)."""