From 099bd2e3eac0608dee9f121c57ee75f7bcf16e46 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fr=C3=A9d=C3=A9ric=20BRIOL?= Date: Sun, 19 Nov 2023 10:58:34 +0100 Subject: [PATCH] =?UTF-8?q?=E2=9C=A8=20Add=20intersection=20method=20to=20?= =?UTF-8?q?LineString=20class?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/pyinterp/core/geodetic.pyi | 7 ++++ .../include/pyinterp/geodetic/line_string.hpp | 7 ++++ .../core/lib/geodetic/line_string.cpp | 17 ++++++++++ src/pyinterp/core/module/geodetic.cpp | 33 ++++++++++++++++--- src/pyinterp/tests/core/test_geodetic.py | 26 +++++++++++++-- 5 files changed, 83 insertions(+), 7 deletions(-) diff --git a/src/pyinterp/core/geodetic.pyi b/src/pyinterp/core/geodetic.pyi index 994240f8..d8402509 100644 --- a/src/pyinterp/core/geodetic.pyi +++ b/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: diff --git a/src/pyinterp/core/include/pyinterp/geodetic/line_string.hpp b/src/pyinterp/core/include/pyinterp/geodetic/line_string.hpp index 39aa497b..df76f429 100644 --- a/src/pyinterp/core/include/pyinterp/geodetic/line_string.hpp +++ b/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 { public: @@ -66,6 +69,10 @@ class LineString : public boost::geometry::model::linestring { const std::optional& wgs) const -> LineString; + [[nodiscard]] auto intersection(const Polygon& rhs, + const std::optional& wgs) const + -> std::vector; + /// Get the length of the linestring /// /// @return the length of the line string diff --git a/src/pyinterp/core/lib/geodetic/line_string.cpp b/src/pyinterp/core/lib/geodetic/line_string.cpp index 957f7bae..1603c30e 100644 --- a/src/pyinterp/core/lib/geodetic/line_string.cpp +++ b/src/pyinterp/core/lib/geodetic/line_string.cpp @@ -7,6 +7,8 @@ #include #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& wgs) const + -> std::vector { + std::vector output; + if (wgs) { + boost::geometry::intersection( + *this, rhs, output, + boost::geometry::strategy::intersection::geographic_segments<>( + static_cast>(*wgs))); + } else { + boost::geometry::intersection(*this, rhs, output); + } + return output; +} + auto LineString::getstate() const -> pybind11::tuple { auto lon = pybind11::array_t(pybind11::array::ShapeContainer{size()}); auto lat = pybind11::array_t(pybind11::array::ShapeContainer{size()}); diff --git a/src/pyinterp/core/module/geodetic.cpp b/src/pyinterp/core/module/geodetic.cpp index 554d41c0..5e068ca1 100644 --- a/src/pyinterp/core/module/geodetic.cpp +++ b/src/pyinterp/core/module/geodetic.cpp @@ -1157,7 +1157,7 @@ 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. @@ -1165,9 +1165,13 @@ Test if this linestring intersects with another linestring. True if the linestring intersects this instance. )__doc__", py::call_guard()) - .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 &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::call_guard()) + .def( + "intersection", + [](const geodetic::LineString &self, const geodetic::Polygon &rhs, + const std::optional &wgs) + -> std::vector { + 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()) .def( "simplify", [](const geodetic::LineString &self, const double tolerance, diff --git a/src/pyinterp/tests/core/test_geodetic.py b/src/pyinterp/tests/core/test_geodetic.py index f8f8ed3d..3c5b2b4d 100644 --- a/src/pyinterp/tests/core/test_geodetic.py +++ b/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)