✨ Add intersection method to LineString class

src/pyinterp/core/geodetic.pyi

@@ -208,11 +208,18 @@ class LineString:
     def from_geojson(array: List[List[float]]) -> LineString:
         ...
 
+    @overload
     def intersection(self,
                      rhs: LineString,
                      wgs: Optional[Spheroid] = None) -> LineString:
         ...
 
+    @overload
+    def intersection(self,
+                     rhs: Polygon,
+                     wgs: Optional[Spheroid] = None) -> List[LineString]:
+        ...
+
     def intersects(self,
                    rhs: LineString,
                    wgs: Optional[Spheroid] = None) -> bool:

src/pyinterp/core/include/pyinterp/geodetic/line_string.hpp

@@ -18,6 +18,9 @@
 
 namespace pyinterp::geodetic {
 
+/// Forward declaration
+class Polygon;
+
 /// A linestring (named so by OGC) is a collection of points.
 class LineString : public boost::geometry::model::linestring<Point> {
  public:
@@ -66,6 +69,10 @@ class LineString : public boost::geometry::model::linestring<Point> {
                                   const std::optional<Spheroid>& wgs) const
       -> LineString;
 
+  [[nodiscard]] auto intersection(const Polygon& rhs,
+                                  const std::optional<Spheroid>& wgs) const
+      -> std::vector<LineString>;
+
   /// Get the length of the linestring
   ///
   /// @return the length of the line string

src/pyinterp/core/lib/geodetic/line_string.cpp

@@ -7,6 +7,8 @@
 #include <boost/geometry/algorithms/simplify.hpp>
 
 #include "pyinterp/detail/broadcast.hpp"
+#include "pyinterp/geodetic/multipolygon.hpp"
+#include "pyinterp/geodetic/polygon.hpp"
 
 namespace pyinterp::geodetic {
 
@@ -54,6 +56,21 @@ auto LineString::intersection(const LineString& rhs,
   return output;
 }
 
+auto LineString::intersection(const Polygon& rhs,
+                              const std::optional<Spheroid>& wgs) const
+    -> std::vector<LineString> {
+  std::vector<LineString> output;
+  if (wgs) {
+    boost::geometry::intersection(
+        *this, rhs, output,
+        boost::geometry::strategy::intersection::geographic_segments<>(
+            static_cast<boost::geometry::srs::spheroid<double>>(*wgs)));
+  } else {
+    boost::geometry::intersection(*this, rhs, output);
+  }
+  return output;
+}
+
 auto LineString::getstate() const -> pybind11::tuple {
   auto lon = pybind11::array_t<double>(pybind11::array::ShapeContainer{size()});
   auto lat = pybind11::array_t<double>(pybind11::array::ShapeContainer{size()});

src/pyinterp/core/module/geodetic.cpp

@@ -1157,17 +1157,21 @@ Computes the curvilinear distance between the points of this instance.
 Test if this linestring intersects with another linestring.
 
 Args:
-    rhs: The linestring to test.
+    rhs: The linestring to use.
     wgs: If specified, searches for the intersection using geographic
     coordinates with the specified spheroid, otherwise searches for the
     intersection using spherical coordinates.
 Returns:
     True if the linestring intersects this instance.
 )__doc__",
            py::call_guard<py::gil_scoped_release>())
-      .def("intersection", &geodetic::LineString::intersection, py::arg("rhs"),
-           py::arg("wgs") = std::nullopt,
-           R"__doc__(
+      .def(
+          "intersection",
+          [](const geodetic::LineString &self, const geodetic::LineString &rhs,
+             const std::optional<geodetic::Spheroid> &wgs)
+              -> geodetic::LineString { return self.intersection(rhs, wgs); },
+          py::arg("rhs"), py::arg("wgs") = std::nullopt,
+          R"__doc__(
 Computes the intersection between this linestring and another linestring.
 
 Args:
@@ -1177,7 +1181,26 @@ Computes the intersection between this linestring and another linestring.
 Returns:
     The intersection between this linestring and the other linestring.
 )__doc__",
-           py::call_guard<py::gil_scoped_release>())
+          py::call_guard<py::gil_scoped_release>())
+      .def(
+          "intersection",
+          [](const geodetic::LineString &self, const geodetic::Polygon &rhs,
+             const std::optional<geodetic::Spheroid> &wgs)
+              -> std::vector<geodetic::LineString> {
+            return self.intersection(rhs, wgs);
+          },
+          py::arg("rhs"), py::arg("wgs") = std::nullopt,
+          R"__doc__(
+Computes the intersection between this linestring and a polygon.
+
+Args:
+    rhs: The polygon to use.
+    wgs: The World Geodetic System to use. Defaults to WGS84.
+
+Returns:
+    The intersection between this linestring and the polygon.
+)__doc__",
+          py::call_guard<py::gil_scoped_release>())
       .def(
           "simplify",
           [](const geodetic::LineString &self, const double tolerance,

src/pyinterp/tests/core/test_geodetic.py

@@ -2,6 +2,7 @@
 #
 # All rights reserved. Use of this source code is governed by a
 # BSD-style license that can be found in the LICENSE file.
+from typing import Tuple
 import json
 import math
 import pickle
@@ -701,8 +702,8 @@ def test_linestring():
     assert a == b
 
 
-def test_curvilinear_distance():
-    """Test the curvilinear distance calculation."""
+def get_coordinates() -> Tuple[np.ndarray, np.ndarray]:
+    """Get coordinates for a line string."""
     lon = np.array([
         -9.72270435, -9.58849868, -9.45429341, -9.32009001, -9.18588996,
         -9.05169472, -8.91750577, -8.78332458, -8.64915263, -8.51499137,
@@ -727,6 +728,12 @@ def test_curvilinear_distance():
         -77.6090989, -77.60641495, -77.60366627, -77.60085292, -77.59797493,
         -77.59503237, -77.59202528, -77.5889537, -77.5858177, -77.58261731
     ])
+    return lon, lat
+
+
+def test_curvilinear_distance():
+    """Test the curvilinear distance calculation."""
+    lon, lat = get_coordinates()
     ls = geodetic.LineString(lon, lat)
 
     distance = [0.0]
@@ -739,3 +746,18 @@ def test_curvilinear_distance():
 
     assert np.all((ls.curvilinear_distance(strategy='thomas') -
                    np.cumsum(np.array(distance))) == np.zeros(50))
+
+
+def test_linestring_intersection_with_polygon():
+    """Test the intersection of a line string with a polygon."""
+    lon, lat = get_coordinates()
+    ls = geodetic.LineString(lon, lat)
+    polygon = geodetic.Polygon.read_wkt(
+        'POLYGON((-60 -180,-60 180,60 180,60 -180,-60 -180))')
+    intersection = ls.intersection(polygon)
+    assert len(intersection) == 0
+    polygon = geodetic.Polygon.read_wkt(
+        'POLYGON((-90 -180,-90 180,90 180,90 -180,-90 -180))')
+    intersection = ls.intersection(polygon)
+    assert len(intersection) == 1
+    assert intersection[0] == geodetic.LineString(lon, lat)