feat: add for_each_point_pairwise_distance algorithm for pairwise distance calculations between geometries

Author: fbriol

cxx/include/pyinterp/geometry/geographic/algorithms/for_each_point_distance.hpp

@@ -35,6 +35,31 @@ template <typename SourceGeometry, typename TargetGeometry,
   return result;
 }
 
+/// @brief Calculate distance from each point to target using compile-time
+/// strategy for multi-point target
+/// @tparam Geometry Geometry type (MultiPoint, LineString, Ring)
+/// @tparam Method Strategy method to use
+/// @param[in] source Source geometry containing points
+/// @param[in] target Target geometry containing points (must have the same
+/// number of points as source)
+/// @param[in] spheroid Optional Spheroid for geodetic calculations
+/// @return Vector of distances for each point in meters
+template <typename Geometry, StrategyMethod Method>
+[[nodiscard]] inline auto for_each_point_pairwise_distance(
+    const Geometry& source, const Geometry& target,
+    const std::optional<Spheroid>& spheroid) -> Eigen::VectorXd {
+  if (source.size() != target.size()) {
+    throw std::invalid_argument(
+        "Source and target geometries must have the same number of points.");
+  }
+  Eigen::VectorXd result(source.size());
+  auto strategy = make_distance_strategy<Method>(make_spheroid(spheroid));
+  for (std::size_t i = 0; i < source.size(); ++i) {
+    result(i) = boost::geometry::distance(source[i], target[i], strategy);
+  }
+  return result;
+}
+
 /// @brief Calculate distance from each point to target using runtime strategy
 /// @tparam SourceGeometry Source geometry type (MultiPoint, LineString, Ring)
 /// @tparam TargetGeometry Target geometry type
@@ -66,4 +91,36 @@ template <typename SourceGeometry, typename TargetGeometry>
   std::unreachable();
 }
 
+/// @brief Calculate distance from each point to target using runtime strategy
+/// for multi-point target
+/// @tparam Geometry Geometry type (MultiPoint, LineString, Ring)
+/// @param[in] source Source geometry containing points
+/// @param[in] target Target geometry containing points (must have the same
+/// number of points as source)
+/// @param[in] spheroid Optional Spheroid for geodetic calculations
+/// @param[in] strategy Strategy method to use
+/// @return Vector of distances for each point in meters
+template <typename Geometry>
+[[nodiscard]] inline auto for_each_point_pairwise_distance(
+    const Geometry& source, const Geometry& target,
+    const std::optional<Spheroid>& spheroid, const StrategyMethod strategy)
+    -> Eigen::VectorXd {
+  using enum StrategyMethod;
+  switch (strategy) {
+    case kAndoyer:
+      return for_each_point_pairwise_distance<Geometry, kAndoyer>(
+          source, target, spheroid);
+    case kKarney:
+      return for_each_point_pairwise_distance<Geometry, kKarney>(source, target,
+                                                                 spheroid);
+    case kThomas:
+      return for_each_point_pairwise_distance<Geometry, kThomas>(source, target,
+                                                                 spheroid);
+    case kVincenty:
+      return for_each_point_pairwise_distance<Geometry, kVincenty>(
+          source, target, spheroid);
+  }
+  std::unreachable();
+}
+
 }  // namespace pyinterp::geometry::geographic

cxx/src/pybind/geometry/cartesian/algorithm/for_each_point_distance_cartesian.cpp

@@ -6,6 +6,11 @@
 #include <nanobind/eigen/dense.h>
 #include <nanobind/nanobind.h>
 
+#include <utility>
+
+#include "pyinterp/geometry/cartesian/linestring.hpp"
+#include "pyinterp/geometry/cartesian/multi_point.hpp"
+#include "pyinterp/geometry/cartesian/ring.hpp"
 #include "pyinterp/pybind/geometry/algorithm_binding_helpers.hpp"
 
 namespace nb = nanobind;
@@ -24,6 +29,19 @@ Calculate the distance from each point in a source geometry to a target geometry
     Array of distances in coordinate units.
 )doc";
 
+constexpr auto kForEachPointPairwiseDistanceDoc = R"doc(
+Calculate pairwise distances between points of two geometries.
+
+Args:
+    geometry1: Source geometry containing points (MultiPoint, LineString, or
+      Ring).
+    geometry2: Target geometry containing points (must have the same number of
+      points as geometry1).
+
+Returns:
+    Array of distances in coordinate units.
+)doc";
+
 // Calculate distance from each point in source to target geometry
 template <typename SourceGeometry, typename TargetGeometry>
 [[nodiscard]] inline auto for_each_point_distance(const SourceGeometry& source,
@@ -36,13 +54,34 @@ template <typename SourceGeometry, typename TargetGeometry>
   return result;
 }
 
+// Calculate pairwise distances for geometries of the same type.
+template <typename Geometry>
+[[nodiscard]] inline auto for_each_point_pairwise_distance(
+    const Geometry& source, const Geometry& target) -> Eigen::VectorXd {
+  if (source.size() != target.size()) {
+    throw std::invalid_argument(
+        "Source and target geometries must have the same number of points.");
+  }
+  Eigen::VectorXd result(source.size());
+  for (std::size_t i = 0; i < source.size(); ++i) {
+    result(i) = boost::geometry::distance(source[i], target[i]);
+  }
+  return result;
+}
+
 auto init_for_each_point_distance(nb::module_& m) -> void {
   auto distance_impl = [](const auto& source,
                           const auto& target) -> Eigen::VectorXd {
     nb::gil_scoped_release release;
     return for_each_point_distance(source, target);
   };
 
+  auto pairwise_distance_impl = [](const auto& source,
+                                   const auto& target) -> Eigen::VectorXd {
+    nb::gil_scoped_release release;
+    return for_each_point_pairwise_distance(source, target);
+  };
+
   // Bind for MultiPoint
   geometry::pybind::define_for_each_point_single_source<
       decltype(distance_impl), MultiPoint, CONTAINER_TYPES(cartesian)>(
@@ -57,6 +96,15 @@ auto init_for_each_point_distance(nb::module_& m) -> void {
   geometry::pybind::define_for_each_point_single_source<
       decltype(distance_impl), Ring, CONTAINER_TYPES(cartesian)>(
       m, "for_each_point_distance", kForEachPointDistanceDoc, distance_impl);
+
+  // Bind pairwise distances for same-type geometries
+  geometry::pybind::define_binary_predicate<
+      decltype(pairwise_distance_impl),
+      std::pair<cartesian::MultiPoint, cartesian::MultiPoint>,
+      std::pair<cartesian::Ring, cartesian::Ring>,
+      std::pair<cartesian::LineString, cartesian::LineString>>(
+      m, "for_each_point_pairwise_distance", kForEachPointPairwiseDistanceDoc,
+      std::move(pairwise_distance_impl));
 }
 
 }  // namespace pyinterp::geometry::cartesian::pybind

cxx/src/pybind/geometry/geographic/algorithm/for_each_point_distance_geographic.cpp

@@ -33,6 +33,22 @@ trade-offs.
     Array of distances in meters.
 )doc";
 
+constexpr auto kForEachPointPairwiseDistanceDoc = R"doc(
+Calculate pairwise distances between points of two geometries.
+
+Args:
+    geometry1: Source geometry containing points (MultiPoint, LineString, or
+        Ring).
+    geometry2: Target geometry containing points (must have the same number of
+        points as geometry1).
+    spheroid: Optional spheroid for geodetic calculations. If not provided, uses
+        WGS84 ellipsoid.
+    strategy: Calculation strategy.
+
+Returns:
+    Array of distances in meters.
+)doc";
+
 auto init_for_each_point_distance(nb::module_& m) -> void {
   auto distance_impl = [](const auto& source, const auto& target,
                           const std::optional<Spheroid>& spheroid,
@@ -41,6 +57,14 @@ auto init_for_each_point_distance(nb::module_& m) -> void {
     return for_each_point_distance(source, target, spheroid, strategy);
   };
 
+  auto pairwise_distance_impl =
+      [](const auto& source, const auto& target,
+         const std::optional<Spheroid>& spheroid,
+         const StrategyMethod strategy) -> Eigen::VectorXd {
+    nb::gil_scoped_release release;
+    return for_each_point_pairwise_distance(source, target, spheroid, strategy);
+  };
+
   // Bind for MultiPoint
   geometry::pybind::define_for_each_point_single_source_with_strategy<
       decltype(distance_impl), MultiPoint, Spheroid, StrategyMethod,
@@ -58,6 +82,15 @@ auto init_for_each_point_distance(nb::module_& m) -> void {
       decltype(distance_impl), Ring, Spheroid, StrategyMethod,
       CONTAINER_TYPES(geographic)>(m, "for_each_point_distance",
                                    kForEachPointDistanceDoc, distance_impl);
+
+  // Bind pairwise distances for same-type geometries
+  geometry::pybind::define_binary_predicate_with_strategy<
+      decltype(pairwise_distance_impl), Spheroid, StrategyMethod,
+      std::pair<geographic::MultiPoint, geographic::MultiPoint>,
+      std::pair<geographic::Ring, geographic::Ring>,
+      std::pair<geographic::LineString, geographic::LineString>>(
+      m, "for_each_point_pairwise_distance", kForEachPointPairwiseDistanceDoc,
+      std::move(pairwise_distance_impl));
 }
 
 }  // namespace pyinterp::geometry::geographic::pybind

docs/source/api/geometry.rst

@@ -122,6 +122,7 @@ Vectorized spatial predicates for geographic geometries.
    algorithms.for_each_point_within
    algorithms.for_each_point_covered_by
    algorithms.for_each_point_distance
+   algorithms.for_each_point_pairwise_distance
 
 Distance Strategies
 ^^^^^^^^^^^^^^^^^^^
@@ -220,6 +221,7 @@ Vectorized spatial predicates for Cartesian geometries.
    algorithms.for_each_point_within
    algorithms.for_each_point_covered_by
    algorithms.for_each_point_distance
+   algorithms.for_each_point_pairwise_distance
 
 Buffer Strategies
 ^^^^^^^^^^^^^^^^^

pyinterp/core/geometry/cartesian/algorithms.pyi

@@ -334,6 +334,18 @@ def for_each_point_distance(
     source: MultiPoint | LineString | Ring,
     container: Box | Ring | Polygon | MultiPolygon,
 ) -> NDArray1DFloat64: ...
+@overload
+def for_each_point_pairwise_distance(
+    geometry1: Ring, geometry2: Ring
+) -> NDArray1DFloat64: ...
+@overload
+def for_each_point_pairwise_distance(
+    geometry1: MultiPoint, geometry2: MultiPoint
+) -> NDArray1DFloat64: ...
+@overload
+def for_each_point_pairwise_distance(
+    geometry1: LineString, geometry2: LineString
+) -> NDArray1DFloat64: ...
 def for_each_point_within(
     source: MultiPoint | LineString | Ring,
     container: Box | Ring | Polygon | MultiPolygon,

pyinterp/core/geometry/geographic/algorithms.pyi

@@ -651,6 +651,27 @@ def for_each_point_distance(
     spheroid: Spheroid | None = None,
     strategy: Strategy = ...,
 ) -> NDArray1DFloat64: ...
+@overload
+def for_each_point_pairwise_distance(
+    geometry1: Ring,
+    geometry2: Ring,
+    spheroid: Spheroid | None = None,
+    strategy: Strategy = Strategy.ANDOYER,
+) -> NDArray1DFloat64: ...
+@overload
+def for_each_point_pairwise_distance(
+    geometry1: MultiPoint,
+    geometry2: MultiPoint,
+    spheroid: Spheroid | None = None,
+    strategy: Strategy = Strategy.ANDOYER,
+) -> NDArray1DFloat64: ...
+@overload
+def for_each_point_pairwise_distance(
+    geometry1: LineString,
+    geometry2: LineString,
+    spheroid: Spheroid | None = None,
+    strategy: Strategy = Strategy.ANDOYER,
+) -> NDArray1DFloat64: ...
 def for_each_point_within(
     source: MultiPoint | LineString | Ring,
     container: Box | Ring | Polygon | MultiPolygon,

pyinterp/tests/core/geometry/cartesian/algorithm/for_each_points/test_distance.py

@@ -16,7 +16,10 @@
     Polygon,
     Ring,
 )
-from .......core.geometry.cartesian.algorithms import for_each_point_distance
+from .......core.geometry.cartesian.algorithms import (
+    for_each_point_distance,
+    for_each_point_pairwise_distance,
+)
 
 
 class TestForEachPointDistance:
@@ -141,3 +144,65 @@ def test_for_each_point_distance_all_inside(self) -> None:
         assert len(result) == 4
         # All distances should be 0
         assert np.allclose(result, 0.0)
+
+
+class TestForEachPointPairwiseDistance:
+    """Tests for for_each_point_pairwise_distance algorithm."""
+
+    def test_for_each_point_pairwise_distance_multipoint(self) -> None:
+        """Test pairwise distances for multipoints."""
+        geometry1 = MultiPoint(
+            np.array([0.0, 1.0, 4.0]), np.array([0.0, 1.0, 4.0])
+        )
+        geometry2 = MultiPoint(
+            np.array([0.0, 2.0, 4.0]), np.array([0.0, 2.0, 7.0])
+        )
+
+        result = for_each_point_pairwise_distance(geometry1, geometry2)
+
+        assert isinstance(result, np.ndarray)
+        assert result.dtype == np.float64
+        assert len(result) == 3
+        assert np.allclose(result, np.array([0.0, np.sqrt(2.0), 3.0]))
+
+    def test_for_each_point_pairwise_distance_linestring(self) -> None:
+        """Test pairwise distances for linestrings."""
+        geometry1 = LineString(np.array([0.0, 2.0]), np.array([0.0, 2.0]))
+        geometry2 = LineString(np.array([0.0, 2.0]), np.array([0.0, 5.0]))
+
+        result = for_each_point_pairwise_distance(geometry1, geometry2)
+
+        assert isinstance(result, np.ndarray)
+        assert result.dtype == np.float64
+        assert len(result) == 2
+        assert np.allclose(result, np.array([0.0, 3.0]))
+
+    def test_for_each_point_pairwise_distance_ring(self) -> None:
+        """Test pairwise distances for rings."""
+        geometry1 = Ring(
+            np.array([0.0, 0.0, 1.0, 1.0, 0.0]),
+            np.array([0.0, 1.0, 1.0, 0.0, 0.0]),
+        )
+        geometry2 = Ring(
+            np.array([0.0, 0.0, 2.0, 2.0, 0.0]),
+            np.array([0.0, 1.0, 1.0, 0.0, 0.0]),
+        )
+
+        result = for_each_point_pairwise_distance(geometry1, geometry2)
+
+        assert isinstance(result, np.ndarray)
+        assert result.dtype == np.float64
+        assert len(result) == 5
+        assert np.allclose(result, np.array([0.0, 0.0, 1.0, 1.0, 0.0]))
+
+    def test_for_each_point_pairwise_distance_size_mismatch(self) -> None:
+        """Test pairwise distances with different number of points."""
+        geometry1 = MultiPoint(np.array([0.0, 1.0]), np.array([0.0, 1.0]))
+        geometry2 = MultiPoint(np.array([0.0]), np.array([0.0]))
+
+        with pytest.raises(
+            ValueError,
+            match="Source and target geometries must have the same "
+            "number of points",
+        ):
+            for_each_point_pairwise_distance(geometry1, geometry2)