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

@@ -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()

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

src/xtgeo/surface/regular_surface.py

@@ -2043,18 +2043,26 @@ 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:
             poly (Polygons): A XTGeo Polygons instance
             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):

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