diff --git a/src/xtgeo/xyz/_xyz.py b/src/xtgeo/xyz/_xyz.py index 8183f3a6a..f40bc4fce 100644 --- a/src/xtgeo/xyz/_xyz.py +++ b/src/xtgeo/xyz/_xyz.py @@ -1,10 +1,12 @@ # -*- coding: utf-8 -*- """XTGeo XYZ module (abstract base class)""" +import warnings from abc import ABC, abstractmethod import numpy as np import pandas as pd +import xtgeo from xtgeo.common import XTGDescription, XTGeoDialog from . import _xyz_oper @@ -30,6 +32,14 @@ class XYZ(ABC): :class:`Polygons` classes! """ + _NEW_POLYGON_BEHAVIOR = False + + @classmethod + def new_polygon_behaviour(cls, value: bool = False): + """Set default behaviour for inside/outside polygons, set_inside(), etc""" + cls._NEW_POLYGON_BEHAVIOR = value + print(f"Behaviour set to {value}") + def __init__( self, xname: str = "X_UTME", @@ -271,7 +281,31 @@ 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: "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. + """ + _xyz_oper.mark_in_polygons_mpl(self, poly, name, inside_value, outside_value) + + def _getversion_polygon_behaviour(self): + """Private routine, to detect if old or new bheaviour is ran""" + if not self._NEW_POLYGON_BEHAVIOR: + warnings.warn( + "From xtgeo version 5, the behaviour will change!", DeprecationWarning + ) + return 1 + + return 2 + + def operation_polygons(self, poly, value, opname="add", inside=True, version=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' @@ -296,98 +330,112 @@ def operation_polygons(self, poly, value, opname="add", inside=True): opname (str): Name of operation... 'add', 'sub', etc inside (bool): If True do operation inside polygons; else outside. Note that boundary is treated as 'inside' + version (int): 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 remove all points of two non-overlapping polygons are given as input. + + ``version=2``: This function works intuitively! (from August 2023) """ - _xyz_oper.operation_polygons(self, poly, value, opname=opname, inside=inside) + if version == 1: + _xyz_oper.operation_polygons_v1( + self, poly, value, opname=opname, inside=inside + ) + else: + _xyz_oper.operation_polygons_v2( + self, poly, value, opname=opname, inside=inside + ) def add_inside(self, poly, value): """Add a value (scalar) to points inside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, value, opname="add", inside=True) + ver = self._getversion_polygon_behaviour() + self.operation_polygons(poly, value, opname="add", inside=True, version=ver) def add_outside(self, poly, value): """Add a value (scalar) to points outside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, value, opname="add", inside=False) + self.operation_polygons(poly, value, opname="add", inside=False, version=1) def sub_inside(self, poly, value): """Subtract a value (scalar) to points inside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, value, opname="sub", inside=True) + self.operation_polygons(poly, value, opname="sub", inside=True, version=1) def sub_outside(self, poly, value): """Subtract a value (scalar) to points outside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, value, opname="sub", inside=False) + self.operation_polygons(poly, value, opname="sub", inside=False, version=1) def mul_inside(self, poly, value): """Multiply a value (scalar) to points inside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, value, opname="mul", inside=True) + self.operation_polygons(poly, value, opname="mul", inside=True, version=1) def mul_outside(self, poly, value): """Multiply a value (scalar) to points outside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, value, opname="mul", inside=False) + self.operation_polygons(poly, value, opname="mul", inside=False, version=1) def div_inside(self, poly, value): """Divide a value (scalar) to points inside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, value, opname="div", inside=True) + self.operation_polygons(poly, value, opname="div", inside=True, version=1) 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`. """ - self.operation_polygons(poly, value, opname="div", inside=False) + self.operation_polygons(poly, value, opname="div", inside=False, version=1) def set_inside(self, poly, value): """Set a value (scalar) to points inside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, value, opname="set", inside=True) + self.operation_polygons(poly, value, opname="set", inside=True, version=1) def set_outside(self, poly, value): """Set a value (scalar) to points outside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, value, opname="set", inside=False) + self.operation_polygons(poly, value, opname="set", inside=False, version=1) def eli_inside(self, poly): """Eliminate current points inside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, 0, opname="eli", inside=True) + self.operation_polygons(poly, 0, opname="eli", inside=True, version=1) def eli_outside(self, poly): """Eliminate current points outside polygons. See notes under :meth:`operation_polygons`. """ - self.operation_polygons(poly, 0, opname="eli", inside=False) + self.operation_polygons(poly, 0, opname="eli", inside=False, version=1) diff --git a/src/xtgeo/xyz/_xyz_oper.py b/src/xtgeo/xyz/_xyz_oper.py index 10b51ab4d..50c5100e9 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,25 @@ # 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 + + self.dataframe[name] = outside_value + + idgroups = poly.dataframe.groupby(poly.pname) + for _, grp in idgroups: + # print(grp[poly.xname].values) + # print(grp[poly.yname].values.dtype) + singlepoly = np.array([grp[poly.xname].values, grp[poly.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 +49,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. """ - 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} @@ -80,6 +101,54 @@ 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. + + 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, approx 3x. + """ + logger.info("Operations of points inside polygon(s)...") + 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[self.zname][tmp._TMP == 1] *= value + elif opname == "div": + if value != 0: + self.dataframe[self.zname][tmpdf._TMP == 1] /= value + elif opname == "set": + self.dataframe[self.zname][tmpdf._TMP == 1] = value + elif opname == "eli": + self.dataframe[self.zname][tmpdf._TMP == 1] = np.nan + self.dataframe.dropna(how="any", subset=[self.zname], inplace=True) + 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..23a9ce24b 100644 --- a/tests/test_xyz/test_points_poly.py +++ b/tests/test_xyz/test_points_poly.py @@ -50,6 +50,17 @@ ] +def test_mark_points_polygon(testpath): + """Import XYZ points and get marks of points inside or outside""" + + poi = xtgeo.points_from_file(testpath / POINTSET2) + pol = xtgeo.polygons_from_file(testpath / POLSET2) + + # remove points in polygon + poi.mark_in_polygons(pol, name="_ptest", inside_value=1, outside_value=0) + print(poi.dataframe) + + def test_points_in_polygon(testpath): """Import XYZ points and do operations if inside or outside""" @@ -58,14 +69,24 @@ def test_points_in_polygon(testpath): assert poi.nrow == 30 # remove points in polygon - poi.operation_polygons(pol, 0, opname="eli") + poi.operation_polygons(pol, 0, opname="eli", version=2) + poi.to_file("_tmp1.pol") 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=2) + poi.to_file("_tmp2.pol") + # assert poi.nrow == 1 + + +def test_add_inside_old_new_behaviour(testpath): + """Add inside using old and new class behaviour""" + + xtgeo.Points.new_polygon_behaviour(False) + poi = xtgeo.points_from_file(testpath / POINTSET2) + pol = xtgeo.polygons_from_file(testpath / POLSET2) def test_check_single_point_in_polygon():