Add InterpolatePoint::point_at_distance_between for line_measures.

georust · Oct 28, 2024 · a524161 · a524161
1 parent aa6964d
commit a524161
Show file tree

Hide file tree

Showing 6 changed files with 177 additions and 9 deletions.
diff --git a/geo/CHANGES.md b/geo/CHANGES.md
@@ -71,6 +71,8 @@
   line_string.densify::<Geodesic>();
   line_string.densify::<Rhumb>();
    ```
+* Added `InterpolatePoint::point_at_distance_between` for line_measures.
+  * <https://github.com/georust/geo/pull/1235>
 * Change IntersectionMatrix::is_equal_topo to now consider empty geometries as equal.
   * <https://github.com/georust/geo/pull/1223>
 * Fix `(LINESTRING EMPTY).contains(LINESTRING EMPTY)` and `(MULTIPOLYGON EMPTY).contains(MULTIPOINT EMPTY)` which previously

diff --git a/geo/src/algorithm/line_measures/interpolate_point.rs b/geo/src/algorithm/line_measures/interpolate_point.rs
@@ -2,13 +2,19 @@ use crate::{CoordFloat, Point};
 
 /// Interpolate a `Point` along a line between two existing points
 pub trait InterpolatePoint<F: CoordFloat> {
-    /// Returns a new Point along a line between two existing points
+    /// Returns a new Point along a line between two existing points.
     ///
     /// See [specific implementations](#implementors) for details.
-    fn point_at_ratio_between(start: Point<F>, end: Point<F>, ratio_from_start: F) -> Point<F>;
+    fn point_at_distance_between(
+        start: Point<F>,
+        end: Point<F>,
+        distance_from_start: F,
+    ) -> Point<F>;
 
-    // TODO:
-    // fn point_at_distance_between(start: Point<F>, end: Point<F>, distance_from_start: F) -> Point<F>;
+    /// Returns a new Point along a line between two existing points.
+    ///
+    /// See [specific implementations](#implementors) for details.
+    fn point_at_ratio_between(start: Point<F>, end: Point<F>, ratio_from_start: F) -> Point<F>;
 
     /// Interpolates `Point`s along a line between `start` and `end`.
     ///

diff --git a/geo/src/algorithm/line_measures/metric_spaces/euclidean/mod.rs b/geo/src/algorithm/line_measures/metric_spaces/euclidean/mod.rs
@@ -22,12 +22,68 @@ use num_traits::FromPrimitive;
 /// [metric spaces]: super
 pub struct Euclidean;
 
+/// Interpolate Point(s) along a line on the [Euclidean plane].
+///
+/// [Euclidean plane]: https://en.wikipedia.org/wiki/Euclidean_plane
 impl<F: CoordFloat + FromPrimitive> InterpolatePoint<F> for Euclidean {
+    /// Returns the point at the given distance along the line between `start` and `end`.
+    ///
+    /// # Units
+    /// - `distance`: Measured in whatever units your `start` and `end` points use.
+    ///
+    ///   `distance` and your `start` and `end` points should have non-agular
+    ///   units, like meters or miles, **not** lon/lat.
+    ///   For lon/lat points, use the [`Haversine`] or [`Geodesic`] [metric spaces].
+    ///
+    /// [`Haversine`]: crate::line_measures::Haversine
+    /// [`Geodesic`]: crate::line_measures::Geodesic
+    /// [metric spaces]: crate::line_measures::metric_spaces
+    fn point_at_distance_between(
+        start: Point<F>,
+        end: Point<F>,
+        distance_from_start: F,
+    ) -> Point<F> {
+        let diff = end - start;
+        let total_distance = diff.x().hypot(diff.y());
+        let offset = diff * distance_from_start / total_distance;
+        start + offset
+    }
+
+    /// Returns the point at the given ratio along the line between `start` and `end`.
+    ///
+    /// # Units
+    /// - `distance`: Measured in whatever units your `start` and `end` points use.
+    ///
+    ///   `distance` and your `start` and `end` points should have non-agular
+    ///   units, like meters or miles, **not** lon/lat.
+    ///   For lon/lat points, use the [`Haversine`] or [`Geodesic`] [metric spaces].
+    ///
+    /// [`Haversine`]: crate::line_measures::Haversine
+    /// [`Geodesic`]: crate::line_measures::Geodesic
+    /// [metric spaces]: crate::line_measures::metric_spaces
     fn point_at_ratio_between(start: Point<F>, end: Point<F>, ratio_from_start: F) -> Point<F> {
         let diff = end - start;
         start + diff * ratio_from_start
     }
 
+    /// Interpolates `Point`s along a line between `start` and `end`.
+    ///
+    /// As many points as necessary will be added such that the distance between points
+    /// never exceeds `max_distance`. If the distance between start and end is less than
+    /// `max_distance`, no additional points will be included in the output.
+    ///
+    /// `include_ends`: Should the start and end points be included in the output?
+    ///
+    /// # Units
+    /// - `max_distance`: Measured in whatever units your `start` and `end` points use.
+    ///
+    ///   `max_distance` and your `start` and `end` points should have non-agular
+    ///   units, like meters or miles, **not** lon/lat.
+    ///   For lon/lat points, use the [`Haversine`] or [`Geodesic`] [metric spaces].
+    ///
+    /// [`Haversine`]: crate::line_measures::Haversine
+    /// [`Geodesic`]: crate::line_measures::Geodesic
+    /// [metric spaces]: crate::line_measures::metric_spaces
     fn points_along_line(
         start: Point<F>,
         end: Point<F>,
@@ -68,5 +124,21 @@ mod tests {
                 distance.round()
             );
         }
+
+        #[test]
+        fn test_point_at_distance_between() {
+            let new_york_city = Point::new(-8_238_310.24, 4_942_194.78);
+            // web mercator
+            let london = Point::new(-14_226.63, 6_678_077.70);
+            let start = MetricSpace::point_at_distance_between(new_york_city, london, 0.0);
+            assert_relative_eq!(new_york_city, start);
+
+            let midway =
+                MetricSpace::point_at_distance_between(new_york_city, london, 8_405_286.0 / 2.0);
+            assert_relative_eq!(Point::new(-4_126_268., 5_810_136.), midway, epsilon = 1.0);
+
+            let end = MetricSpace::point_at_distance_between(new_york_city, london, 8_405_286.0);
+            assert_relative_eq!(london, end, epsilon = 1.0);
+        }
     }
 }
diff --git a/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs b/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs
@@ -134,6 +134,44 @@ impl Distance<f64, Point<f64>, Point<f64>> for Geodesic {
 ///
 /// [geodesic line]: https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid
 impl InterpolatePoint<f64> for Geodesic {
+    /// Returns a new Point along a [geodesic line] between two existing points on an ellipsoidal model of the earth.
+    ///
+    /// # Units
+    /// `meters_from_start`: meters
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// # use approx::assert_relative_eq;
+    /// use geo::{Geodesic, InterpolatePoint};
+    /// use geo::Point;
+    ///
+    ///
+    /// let p1 = Point::new(10.0, 20.0);
+    /// let p2 = Point::new(125.0, 25.0);
+    ///
+    /// let closer_to_p1 = Geodesic::point_at_distance_between(p1, p2, 100_000.0);
+    /// assert_relative_eq!(closer_to_p1, Point::new(10.81, 20.49), epsilon = 1.0e-2);
+    ///
+    /// let closer_to_p2 = Geodesic::point_at_distance_between(p1, p2, 10_000_000.0);
+    /// assert_relative_eq!(closer_to_p2, Point::new(112.20, 30.67), epsilon = 1.0e-2);
+    /// ```
+    ///
+    /// # References
+    ///
+    /// This uses the geodesic methods given by [Karney (2013)].
+    ///
+    /// [geodesic line]: https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid
+    /// [Karney (2013)]:  https://arxiv.org/pdf/1109.4448.pdf
+    fn point_at_distance_between(
+        start: Point<f64>,
+        end: Point<f64>,
+        meters_from_start: f64,
+    ) -> Point<f64> {
+        let bearing = Self::bearing(start, end);
+        Self::destination(start, bearing, meters_from_start)
+    }
+
     /// Returns a new Point along a [geodesic line] between two existing points on an ellipsoidal model of the earth.
     ///
     /// # Examples
@@ -170,9 +208,7 @@ impl InterpolatePoint<f64> for Geodesic {
         let g = geographiclib_rs::Geodesic::wgs84();
         let (total_distance, azi1, _azi2, _a12) = g.inverse(start.y(), start.x(), end.y(), end.x());
         let distance = total_distance * ratio_from_start;
-        let (lat2, lon2) = g.direct(start.y(), start.x(), azi1, distance);
-
-        Point::new(lon2, lat2)
+        Self::destination(start, azi1, distance)
     }
 
     /// Interpolates `Point`s along a [geodesic line] between `start` and `end`.
@@ -335,6 +371,7 @@ mod tests {
                 let midpoint = MetricSpace::point_at_ratio_between(start, end, 0.5);
                 assert_relative_eq!(midpoint, Point::new(65.87936072133309, 37.72225378005785));
             }
+
             #[test]
             fn points_along_line_with_endpoints() {
                 let start = Point::new(10.0, 20.0);
@@ -347,6 +384,7 @@ mod tests {
                 assert_eq!(route.last().unwrap(), &end);
                 assert_relative_eq!(route[1], Point::new(17.878754355562464, 24.466667836189565));
             }
+
             #[test]
             fn points_along_line_without_endpoints() {
                 let start = Point::new(10.0, 20.0);

diff --git a/geo/src/algorithm/line_measures/metric_spaces/haversine.rs b/geo/src/algorithm/line_measures/metric_spaces/haversine.rs
@@ -88,12 +88,12 @@ impl<F: CoordFloat + FromPrimitive> Destination<F> for Haversine {
     /// the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf)
     ///
     /// [great circle]: https://en.wikipedia.org/wiki/Great_circle
-    fn destination(origin: Point<F>, bearing: F, distance: F) -> Point<F> {
+    fn destination(origin: Point<F>, bearing: F, meters: F) -> Point<F> {
         let center_lng = origin.x().to_radians();
         let center_lat = origin.y().to_radians();
         let bearing_rad = bearing.to_radians();
 
-        let rad = distance / F::from(MEAN_EARTH_RADIUS).unwrap();
+        let rad = meters / F::from(MEAN_EARTH_RADIUS).unwrap();
 
         let lat =
             { center_lat.sin() * rad.cos() + center_lat.cos() * rad.sin() * bearing_rad.cos() }
@@ -155,6 +155,31 @@ impl<F: CoordFloat + FromPrimitive> Distance<F, Point<F>, Point<F>> for Haversin
 ///
 /// [great circle]: https://en.wikipedia.org/wiki/Great_circle
 impl<F: CoordFloat + FromPrimitive> InterpolatePoint<F> for Haversine {
+    /// Returns a new Point along a [great circle] between two existing points.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// # use approx::assert_relative_eq;
+    /// use geo::{Haversine, InterpolatePoint};
+    /// use geo::Point;
+    ///
+    /// let p1 = Point::new(10.0, 20.0);
+    /// let p2 = Point::new(125.0, 25.0);
+    ///
+    /// let closer_to_p1 = Haversine::point_at_distance_between(p1, p2, 100_000.0);
+    /// assert_relative_eq!(closer_to_p1, Point::new(10.81, 20.49), epsilon = 1.0e-2);
+    ///
+    /// let closer_to_p2 = Haversine::point_at_distance_between(p1, p2, 10_000_000.0);
+    /// assert_relative_eq!(closer_to_p2, Point::new(112.33, 30.57), epsilon = 1.0e-2);
+    /// ```
+    ///
+    /// [great circle]: https://en.wikipedia.org/wiki/Great_circle
+    fn point_at_distance_between(start: Point<F>, end: Point<F>, meters_from_start: F) -> Point<F> {
+        let bearing = Self::bearing(start, end);
+        Self::destination(start, bearing, meters_from_start)
+    }
+
     /// Returns a new Point along a [great circle] between two existing points.
     ///
     /// # Examples

diff --git a/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs b/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs
@@ -130,6 +130,31 @@ impl<F: CoordFloat + FromPrimitive> Distance<F, Point<F>, Point<F>> for Rhumb {
 ///
 /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line
 impl<F: CoordFloat + FromPrimitive> InterpolatePoint<F> for Rhumb {
+    /// Returns a new Point along a [rhumb line] between two existing points.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// # use approx::assert_relative_eq;
+    /// use geo::{Rhumb, InterpolatePoint};
+    /// use geo::Point;
+    ///
+    /// let p1 = Point::new(10.0, 20.0);
+    /// let p2 = Point::new(125.0, 25.0);
+    ///
+    /// let closer_to_p1 = Rhumb::point_at_distance_between(p1, p2, 100_000.0);
+    /// assert_relative_eq!(closer_to_p1, Point::new(10.96, 20.04), epsilon = 1.0e-2);
+    ///
+    /// let closer_to_p2 = Rhumb::point_at_distance_between(p1, p2, 10_000_000.0);
+    /// assert_relative_eq!(closer_to_p2, Point::new(107.00, 24.23), epsilon = 1.0e-2);
+    /// ```
+    ///
+    /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line
+    fn point_at_distance_between(start: Point<F>, end: Point<F>, meters_from_start: F) -> Point<F> {
+        let bearing = Self::bearing(start, end);
+        Self::destination(start, bearing, meters_from_start)
+    }
+
     /// Returns a new Point along a [rhumb line] between two existing points.
     ///
     /// # Examples