GmMultiPolyIntersector for python

CMakeLists.txt

@@ -161,6 +161,7 @@ set(xmsgrid_py
   xmsgrid/python/xmsgrid_py.cpp
   # Geometry
   xmsgrid/python/geometry/geometry_py.cpp
+  xmsgrid/python/geometry/GmMultiPolyIntersector_py.cpp
   xmsgrid/python/geometry/GmTriSearch_py.cpp
   # Triangulate
   xmsgrid/python/triangulate/triangulate_py.cpp

_package/tests/unit_tests/geometry_tests/multi_poly_intersector_pyt.py

@@ -0,0 +1,141 @@
+"""Tests for the MultiPolyIntersector class."""
+
+__copyright__ = "(C) Copyright Aquaveo 2023"
+__license__ = "See accompanying file LICENSE or https://aquaveo.com/bsd/license.txt"
+
+# 1. Standard python modules
+import math
+import unittest
+
+# 2. Third party modules
+
+# 3. Aquaveo modules
+
+# 4. Local modules
+from xms import grid
+
+
+def _zeroed_z(points) -> list[list[float]]:
+    """Returns the points with 0.0 for z values."""
+    return [[point[0], point[1], 0.0] for point in points]
+
+
+class TestMultiPolyIntersector(unittest.TestCase):
+    """Tests for the MultiPolyIntersector class."""
+
+    def setUp(self):
+        """Runs before each test case."""
+        pass
+
+    def _run_test(self, pt1, pt2, poly_points, polys, poly_ids, t_vals, pts, starting_id=1, query='covered_by'):
+        """Runs the test."""
+        mpi = grid.geometry.MultiPolyIntersector(poly_points, polys, starting_id, query)
+
+        # Traverse the line segment
+        out_poly_ids, out_t_vals, out_pts = mpi.traverse_line_segment(pt1, pt2)
+
+        # Check poly ids
+        assert len(out_poly_ids) == len(poly_ids)
+        assert_str = f'Expected {poly_ids}. Got {out_poly_ids}'
+        assert out_poly_ids == poly_ids, assert_str
+
+        # Check t vals
+        assert len(out_t_vals) == len(t_vals)
+        for out_t_val, t_val in zip(out_t_vals, t_vals):
+            assert_str = f'Expected {t_val}. Got {out_t_val}'
+            assert math.isclose(out_t_val, t_val), assert_str
+
+        # Check points
+        assert len(out_pts) == len(pts)
+        for idx, (out_point, point) in enumerate(zip(_zeroed_z(out_pts), _zeroed_z(pts))):
+            for i in range(3):
+                assert_str = f'Point {idx}[{i}], : Expected {point[i]}. Got {out_point[i]}'
+                assert math.isclose(out_point[i], point[i]), assert_str
+
+    def test_traverse_line_segment_1_out_out(self):
+        r"""A test.
+
+                      (10,10)
+           3-------------2
+           |             |
+         0------------------1
+           |             |
+           |      1      |
+           |             |
+           0-------------1
+        (0,0)
+        """
+        # Use lists to prove that we can
+        pts = [[0, 0, 0], [10, 0, 0], [10, 10, 0], [0, 10, 0]]
+        polys = [[0, 1, 2, 3]]
+        expected_ids = (1, -1)
+        expected_t_vals = (0.08333333333333333, 0.91666666666666667)
+        expected_pts = ((0.0, 5.0, 0.0), (10.0, 5.0, 0.0))
+        self._run_test([-1, 5], [11, 5], pts, polys, expected_ids, expected_t_vals, expected_pts)
+        expected_ids = (0, -1)
+        self._run_test([-1, 5], [11, 5], pts, polys, expected_ids, expected_t_vals, expected_pts, 0, 'intersects')
+
+    def test_traverse_line_segment_1_out_in(self):
+        """A test.
+
+                      (10,10)
+           3-------------2
+           |             |
+         0----------1    |
+           |             |
+           |      1      |
+           |             |
+           0-------------1
+        (0,0)
+        """
+        # Use tuples to prove that we can
+        pts = ((0, 0, 0), (10, 0, 0), (10, 10, 0), (0, 10, 0))
+        polys = [(0, 1, 2, 3)]
+        expected_ids = (1, -1)
+        expected_t_vals = (0.11111111111111111, 1.0)
+        expected_pts = ((0.0, 5.0, 0.0), (8.0, 5.0, 0.0))
+        self._run_test((-1, 5), (8, 5), pts, polys, expected_ids, expected_t_vals, expected_pts)
+        expected_ids = (2, -1)
+        self._run_test((-1, 5), (8, 5), pts, polys, expected_ids, expected_t_vals, expected_pts, 2, 'intersects')
+
+    def test_polygon_from_point(self):
+        r"""A test.
+
+         (0,20)
+           3-------------2
+           |             |
+           |             |
+           |      *1     |      *
+           |             |
+           |             |
+           0------*-----*1------*------6
+                         |             |
+                         |             |
+                         |      *2     |
+                  *      |             |
+                         |             |
+                         4-------------5
+        (0,0)                        (20,0)
+        """
+        pts = [(0, 10, 0), (10, 10, 0), (10, 20, 0), (0, 20, 0), (10, 0, 0), (20, 0, 0), (20, 10, 0)]
+        polys = [(0, 1, 2, 3), (4, 5, 6, 1)]
+        mpi = grid.geometry.MultiPolyIntersector(pts, polys)
+        assert mpi.polygon_from_point((5, 5, 0)) == -1
+        assert mpi.polygon_from_point((5, 10, 0)) == 1
+        assert mpi.polygon_from_point((5, 15, 0)) == 1
+        assert mpi.polygon_from_point((10, 10, 0)) == 1
+        assert mpi.polygon_from_point((15, 15, 0)) == -1
+        assert mpi.polygon_from_point((15, 10, 0)) == 2
+        assert mpi.polygon_from_point((15, 5, 0)) == 2
+        mpi = grid.geometry.MultiPolyIntersector(pts, polys, query='intersects')
+        assert mpi.polygon_from_point((5, 5, 0)) == -1
+        assert mpi.polygon_from_point((5, 10, 0)) == 1
+        assert mpi.polygon_from_point((5, 15, 0)) == 1
+        assert mpi.polygon_from_point((10, 10, 0)) == 1
+        assert mpi.polygon_from_point((15, 15, 0)) == -1
+        assert mpi.polygon_from_point((15, 10, 0)) == 2
+        assert mpi.polygon_from_point((15, 5, 0)) == 2
+
+
+if __name__ == '__main__':
+    unittest.main()

_package/xms/grid/__init__.py

@@ -3,4 +3,4 @@
 from . import triangulate  # NOQA: F401
 from . import ugrid  # NOQA: F401
 
-__version__ = '7.7.6'
+__version__ = '7.8.0'

_package/xms/grid/geometry/__init__.py

@@ -1,3 +1,4 @@
 """Initialize the module."""
 from . import geometry  # NOQA: F401
+from .multi_poly_intersector import MultiPolyIntersector  # NOQA: F401
 from .tri_search import TriSearch  # NOQA: F401

_package/xms/grid/geometry/multi_poly_intersector.py

@@ -0,0 +1,106 @@
+"""Pure Python wrapper for GmMultiPolyIntersector class."""
+
+__copyright__ = "(C) Copyright Aquaveo 2023"
+__license__ = "See accompanying file LICENSE or https://aquaveo.com/bsd/license.txt"
+
+# 1. Standard python modules
+
+# 2. Third party modules
+
+# 3. Aquaveo modules
+
+# 4. Local modules
+from .._xmsgrid.geometry import GmMultiPolyIntersector
+
+
+class MultiPolyIntersector(object):
+    """Intersects a line segment with any number of polygons in 2D and returns the polygons in order with t values."""
+
+    def __init__(self, points, polys, starting_id=1, query='covered_by', **kwargs) -> None:
+        """Constructor.
+
+        Args:
+            points (list): The points that make up the polygon.
+            polys (list): 0-based indexes into a_points array to define polygons. The first point is NOT repeated as
+             the last point.
+            starting_id (int): If the polygon IDs should start at something other than 1, specify the starting value.
+            query (str): The query to use ('covered_by', 'intersects')
+            **kwargs (dict): Generic keyword arguments
+        """
+        if 'instance' not in kwargs:
+            if not points:
+                raise ValueError('points is a required arguments.')
+            if not polys:
+                raise ValueError('polys is a required argument.')
+            if query not in {'covered_by', 'intersects'}:
+                raise ValueError('query must be either "covered_by" or "intersects".')
+            self._instance = GmMultiPolyIntersector(points, polys, starting_id, query)
+        else:
+            if not isinstance(kwargs['instance'], GmMultiPolyIntersector):
+                raise ValueError('"instance" must be of type _xmsgrid.geometry.GmMultiPolyIntersector')
+            self._instance = kwargs['instance']
+
+    def __eq__(self, other) -> bool:
+        """Equality operator.
+
+        Args:
+            other (MultiPolyIntersector): MultiPolyIntersector to compare
+
+        Returns:
+            bool: True if MultiPolyIntersector are equal
+        """
+        other_instance = getattr(other, '_instance', None)
+        if not other_instance or not isinstance(other_instance, GmMultiPolyIntersector):
+            print("not instance or no value")
+            return False
+        return other_instance == self._instance
+
+    def __ne__(self, other) -> bool:
+        """Equality operator.
+
+        Args:
+            other (MultiPolyIntersector): MultiPolyIntersector to compare
+
+        Returns:
+            bool: True if MultiPolyIntersector are not equal
+        """
+        result = self.__eq__(other)
+        return not result
+
+    # Why define these?
+    # def __repr__(self) -> str:
+    #     """Returns a string representation of the MultiPolyIntersector."""
+    #     return "<xms.grid.geometry.MultiPolyIntersector>"
+    #
+    # def __str__(self) -> str:
+    #     """Returns a string representation of the MultiPolyIntersector."""
+    #     return "<xms.grid.geometry.MultiPolyIntersector>"
+
+    def traverse_line_segment(self, pt1, pt2) -> tuple[tuple[int], tuple[float], tuple[tuple[float, float, float]]]:
+        """Intersect segment with polys and return intersected polys, t-values, and points.
+
+        Args:
+            pt1 (iterable): 1st point defining a line segment.
+            pt2 (iterable): 2nd point defining a line segment.
+
+        Returns:
+            tuple containing:
+            - Ids of polygons intersected by line segment. Can be zero or 1 based depending on starting_id.
+            - Values from 0.0 to 1.0 representing where on the line segment the intersection with the polygon occurs.
+            If there are any t values there are always at least 2 and all represent where the line enters the polygon,
+            except the last which represents where it exited. There would therefore be one more t value than poly id
+            but we make the sizes equal by always making the last poly id -1.
+            - Intersection points.
+        """
+        return self._instance.TraverseLineSegment(pt1, pt2)
+
+    def polygon_from_point(self, point) -> int:
+        """Returns the ID of the polygon containing the point.
+
+        Args:
+            point (iterable): The point.
+
+        Returns:
+            The polygon id.
+        """
+        return self._instance.PolygonFromPoint(point)

pydocs/source/index.rst

@@ -36,6 +36,7 @@ Index
 
 * :doc:`modules/ugrid/UGrid`
 * :doc:`modules/ugrid/ugrid_utils`
+* :doc:`modules/geometry/MultiPolyIntersector`
 * :doc:`modules/geometry/TriSearch`
 * :doc:`modules/triangulate/Tin`
 
@@ -46,6 +47,7 @@ Index
 
    modules/ugrid/UGrid.rst
    modules/ugrid/ugrid_utils.rst
+   modules/geometry/MultiPolyIntersector.rst
    modules/geometry/TriSearch.rst
    modules/triangulate/Tin.rst
 

pydocs/source/modules/geometry/MultiPolyIntersector.rst

@@ -0,0 +1,10 @@
+********************
+MultiPolyIntersector
+********************
+
+Intersects a line segment with any number of polygons in 2D and returns the polygons in order with t values.
+
+.. autoclass:: xms.grid.geometry.MultiPolyIntersector
+   :members:
+
+   .. automethod:: __init__

xmsgrid/geometry/GmBoostTypes.h

@@ -4,7 +4,7 @@
 /// \brief boost::geometry types
 /// \ingroup geometry
 /// \copyright (C) Copyright Aquaveo 2018. Distributed under FreeBSD License
-/// (See accompanying file LICENSE or https://aqaveo.com/bsd/license.txt)
+/// (See accompanying file LICENSE or https://aquaveo.com/bsd/license.txt)
 //------------------------------------------------------------------------------
 
 //----- Included files ---------------------------------------------------------

xmsgrid/geometry/GmExtents.cpp

@@ -2,7 +2,7 @@
 /// \file
 /// \ingroup geometry
 /// \copyright (C) Copyright Aquaveo 2018. Distributed under FreeBSD License
-/// (See accompanying file LICENSE or https://aqaveo.com/bsd/license.txt)
+/// (See accompanying file LICENSE or https://aquaveo.com/bsd/license.txt)
 //------------------------------------------------------------------------------
 
 //----- Included files ---------------------------------------------------------

xmsgrid/geometry/GmExtents.h

@@ -3,7 +3,7 @@
 /// \file
 /// \ingroup geometry
 /// \copyright (C) Copyright Aquaveo 2018. Distributed under FreeBSD License
-/// (See accompanying file LICENSE or https://aqaveo.com/bsd/license.txt)
+/// (See accompanying file LICENSE or https://aquaveo.com/bsd/license.txt)
 //------------------------------------------------------------------------------
 
 //----- Included files ---------------------------------------------------------

xmsgrid/geometry/GmExtents.t.h

@@ -1,7 +1,7 @@
 //------------------------------------------------------------------------------
 /// \file
 /// \copyright (C) Copyright Aquaveo 2018. Distributed under FreeBSD License
-/// (See accompanying file LICENSE or https://aqaveo.com/bsd/license.txt)
+/// (See accompanying file LICENSE or https://aquaveo.com/bsd/license.txt)
 //------------------------------------------------------------------------------
 #pragma once
 

xmsgrid/geometry/GmMultiPolyIntersectionSorter.h

@@ -3,7 +3,7 @@
 /// \file
 /// \ingroup geometry
 /// \copyright (C) Copyright Aquaveo 2018. Distributed under FreeBSD License
-/// (See accompanying file LICENSE or https://aqaveo.com/bsd/license.txt)
+/// (See accompanying file LICENSE or https://aquaveo.com/bsd/license.txt)
 //------------------------------------------------------------------------------
 
 //----- Included files ---------------------------------------------------------

xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.cpp

@@ -2,7 +2,7 @@
 /// \file
 /// \ingroup geometry
 /// \copyright (C) Copyright Aquaveo 2018. Distributed under FreeBSD License
-/// (See accompanying file LICENSE or https://aqaveo.com/bsd/license.txt)
+/// (See accompanying file LICENSE or https://aquaveo.com/bsd/license.txt)
 //------------------------------------------------------------------------------
 
 //----- Included files ---------------------------------------------------------

xmsgrid/geometry/GmMultiPolyIntersectionSorterTerse.h

@@ -3,7 +3,7 @@
 /// \file
 /// \ingroup geometry
 /// \copyright (C) Copyright Aquaveo 2018. Distributed under FreeBSD License
-/// (See accompanying file LICENSE or https://aqaveo.com/bsd/license.txt)
+/// (See accompanying file LICENSE or https://aquaveo.com/bsd/license.txt)
 //------------------------------------------------------------------------------
 
 //----- Included files ---------------------------------------------------------

xmsgrid/geometry/GmMultiPolyIntersector.cpp

@@ -2338,7 +2338,7 @@ void GmMultiPolyIntersector2IntermediateTests::testPointOnPolygonVertex()
     {118.66666666666669, 93.33333333333333, 0.0},
     {100.0, 100.0, 0.0}
   };
-  
+
   VecInt2d expectedPolyIds = {
     {167, -1},
     {166, -1},
@@ -2412,7 +2412,7 @@ void GmMultiPolyIntersector2IntermediateTests::testPointOnPolygonVertex()
     {138.66666666666669, 93.33333333333333, 0.0},
     {120.0, 100.0, 0.0}
   };
-  
+
   expectedPolyIds = {
     {},
     {13, -1},
@@ -2497,7 +2497,7 @@ void GmMultiPolyIntersector2IntermediateTests::testPointsNearEdgePoints()
   VecDbl tValues;
   VecPt3d points;
   mpi->TraverseLineSegment(segmentPoints[0].x, segmentPoints[0].y, segmentPoints[1].x, segmentPoints[1].y, polyIds, tValues, points);
-  
+
   VecInt expectedPolyIds = {1836, -1};
   VecPt3d expectedPoints = {
     {958452.39285714, 498553.35714286, 0.0},
@@ -2525,4 +2525,4 @@ void GmMultiPolyIntersector2IntermediateTests::testPointsNearEdgePoints()
 
 //} // namespace xms
 
-#endif
+#endif

xmsgrid/geometry/GmMultiPolyIntersector.h

@@ -77,4 +77,4 @@ class GmMultiPolyIntersector
 
 }; // class GmMultiPolyIntersector
 
-} // namespace xms
+} // namespace xms