From a34590a716af441051ec54a7efb5ba1707abfae0 Mon Sep 17 00:00:00 2001 From: Michael Kirk Date: Tue, 24 Sep 2024 11:17:22 -0700 Subject: [PATCH] Unify line measures All of these algorithm implementations already existed, but I'm proposing a new way to organize them which I hope will be more consistent and discoverable. Note: The new Haversine::bearing and Geodesic::bearing return 0..360, the legacy traits HaversineBearing and GeodesicBearing returned -180..180 Additional changes: Deleted the deprecated `Bearing` trait which was previously superceeded by the unambiguous `HaversineBearing` trait, but now is re-defined as `Haversine::bearing` = Future Work = In an effort to minimize this PR, while keeping the change reasonably coherent, I've left some things out == Methods on Euclidean == -[ ] bearing (doesn't currently exist) -[ ] destination (doesn't currently exist) -[ ] intermediate (exists, but InterpolatePoint trait also needs intermediate_fill) -[ ] intermediate_fill (doesn't currently exist) == Deprecate Legacy Traits == -[ ] Deprecate Legacy impls -[ ] Switcheroo the actual implementation: move the actual implementation to the new traits, and have the legacy traits delegate to the new traits. -[ ] Move over any tests from the legacy implementation to the new home == Methods on Geoms (Future PR) == -[ ] Length -[ ] Haversine -[ ] Rhumb -[ ] Geodesic -[ ] Euclidean -[ ] Densify -[ ] Haversine -[ ] Rhumb -[ ] Geodesic -[ ] Euclidean FIXES https://github.com/georust/geo/issues/1210 FIXES https://github.com/georust/geo/issues/1181 --- geo/src/algorithm/bearing.rs | 33 -- geo/src/algorithm/geodesic_intermediate.rs | 2 - geo/src/algorithm/line_measures/bearing.rs | 13 + .../algorithm/line_measures/destination.rs | 19 + geo/src/algorithm/line_measures/distance.rs | 12 + .../line_measures/interpolate_point.rs | 34 ++ .../line_measures/metric_spaces/euclidean.rs | 71 ++++ .../line_measures/metric_spaces/geodesic.rs | 325 ++++++++++++++++++ .../line_measures/metric_spaces/haversine.rs | 322 +++++++++++++++++ .../line_measures/metric_spaces/mod.rs | 11 + .../line_measures/metric_spaces/rhumb.rs | 309 +++++++++++++++++ geo/src/algorithm/line_measures/mod.rs | 16 + geo/src/algorithm/mod.rs | 14 +- geo/src/algorithm/rhumb/bearing.rs | 1 - geo/src/lib.rs | 3 +- 15 files changed, 1138 insertions(+), 47 deletions(-) delete mode 100644 geo/src/algorithm/bearing.rs create mode 100644 geo/src/algorithm/line_measures/bearing.rs create mode 100644 geo/src/algorithm/line_measures/destination.rs create mode 100644 geo/src/algorithm/line_measures/distance.rs create mode 100644 geo/src/algorithm/line_measures/interpolate_point.rs create mode 100644 geo/src/algorithm/line_measures/metric_spaces/euclidean.rs create mode 100644 geo/src/algorithm/line_measures/metric_spaces/geodesic.rs create mode 100644 geo/src/algorithm/line_measures/metric_spaces/haversine.rs create mode 100644 geo/src/algorithm/line_measures/metric_spaces/mod.rs create mode 100644 geo/src/algorithm/line_measures/metric_spaces/rhumb.rs create mode 100644 geo/src/algorithm/line_measures/mod.rs diff --git a/geo/src/algorithm/bearing.rs b/geo/src/algorithm/bearing.rs deleted file mode 100644 index 40c53591a8..0000000000 --- a/geo/src/algorithm/bearing.rs +++ /dev/null @@ -1,33 +0,0 @@ -use crate::{haversine_bearing, CoordFloat, Point}; - -pub use haversine_bearing::HaversineBearing; - -#[deprecated(since = "0.24.1", note = "renamed to `HaversineBearing`")] -pub trait Bearing { - /// Returns the bearing to another Point in degrees, where North is 0° and East is 90°. - /// - /// # Examples - /// - /// ``` - /// # use approx::assert_relative_eq; - /// use geo::Bearing; - /// use geo::Point; - /// - /// let p_1 = Point::new(9.177789688110352, 48.776781529534965); - /// let p_2 = Point::new(9.274410083250379, 48.84033282787534); - /// let bearing = p_1.bearing(p_2); - /// assert_relative_eq!(bearing, 45., epsilon = 1.0e-6); - /// ``` - #[deprecated( - since = "0.24.1", - note = "renamed to `HaversineBearing::haversine_bearing`" - )] - fn bearing(&self, point: Point) -> T; -} - -#[allow(deprecated)] -impl> Bearing for B { - fn bearing(&self, point: Point) -> T { - self.haversine_bearing(point) - } -} diff --git a/geo/src/algorithm/geodesic_intermediate.rs b/geo/src/algorithm/geodesic_intermediate.rs index a4d43836cb..3670fbba5e 100644 --- a/geo/src/algorithm/geodesic_intermediate.rs +++ b/geo/src/algorithm/geodesic_intermediate.rs @@ -2,7 +2,6 @@ use crate::{CoordFloat, Point}; use geographiclib_rs::{DirectGeodesic, Geodesic, InverseGeodesic}; /// Returns a new Point along a route between two existing points on an ellipsoidal model of the earth - pub trait GeodesicIntermediate { /// Returns a new Point along a route between two existing points on an ellipsoidal model of the earth /// @@ -25,7 +24,6 @@ pub trait GeodesicIntermediate { /// assert_relative_eq!(i50, i50_should, epsilon = 1.0e-6); /// assert_relative_eq!(i80, i80_should, epsilon = 1.0e-6); /// ``` - fn geodesic_intermediate(&self, other: &Point, f: T) -> Point; fn geodesic_intermediate_fill( &self, diff --git a/geo/src/algorithm/line_measures/bearing.rs b/geo/src/algorithm/line_measures/bearing.rs new file mode 100644 index 0000000000..d092ec9e88 --- /dev/null +++ b/geo/src/algorithm/line_measures/bearing.rs @@ -0,0 +1,13 @@ +use geo_types::{CoordFloat, Point}; + +/// Calculate the bearing between two points +pub trait Bearing { + /// Calculate the bearing from `origin` to `destination` in degrees. + /// + /// See [specific implementations](#implementors) for details. + /// + /// # Units + /// - `origin`, `destination`: Point where the units of x/y depend on the [trait implementation](#implementors). + /// - returns: degrees, where: North: 0°, East: 90°, South: 180°, West: 270° + fn bearing(origin: Point, destination: Point) -> F; +} diff --git a/geo/src/algorithm/line_measures/destination.rs b/geo/src/algorithm/line_measures/destination.rs new file mode 100644 index 0000000000..91ffb73ef4 --- /dev/null +++ b/geo/src/algorithm/line_measures/destination.rs @@ -0,0 +1,19 @@ +use geo_types::{CoordFloat, Point}; + +/// Calculate the destination point from an origin point, a bearing and a distance. +pub trait Destination { + /// Returns a new point having travelled the `distance` along a line + /// from the `origin` point with the given `bearing`. + /// + /// See [specific implementations](#implementors) for details. + /// + /// # Units + /// + /// - `origin`: Point where the units of x/y depend on the [trait implementation](#implementors). + /// - `bearing`: degrees, where: North: 0°, East: 90°, South: 180°, West: 270° + /// - `distance`: depends on the [trait implementation](#implementors). + /// - returns: Point where the units of x/y depend on the [trait implementation](#implementors). + /// + /// [`metric_spaces`]: super::metric_spaces + fn destination(origin: Point, bearing: F, distance: F) -> Point; +} diff --git a/geo/src/algorithm/line_measures/distance.rs b/geo/src/algorithm/line_measures/distance.rs new file mode 100644 index 0000000000..e89749fb28 --- /dev/null +++ b/geo/src/algorithm/line_measures/distance.rs @@ -0,0 +1,12 @@ +/// Calculate the distance between the `Origin` and `Destination` geometry. +pub trait Distance { + /// Note that not all implementations support all geometry combinations, but at least `Point` to `Point` + /// is supported. + /// See [specific implementations](#implementors) for details. + /// + /// # Units + /// + /// - `origin`, `destination`: geometry where the units of x/y depend on the trait implementation. + /// - returns: depends on the trait implementation. + fn distance(origin: Origin, destination: Destination) -> F; +} diff --git a/geo/src/algorithm/line_measures/interpolate_point.rs b/geo/src/algorithm/line_measures/interpolate_point.rs new file mode 100644 index 0000000000..018bb0003e --- /dev/null +++ b/geo/src/algorithm/line_measures/interpolate_point.rs @@ -0,0 +1,34 @@ +use crate::{CoordFloat, Point}; + +// REVIEW: Naming alternatives: +// - LinearReferencing +// - PointAlongLine +// - LineInterpolatePoint (postgis) +// - Interpolate (shapely) +// - Position (geographiclib) +// - Intermediate (georust::geo) +pub trait InterpolatePoint { + /// Returns a new Point along a line between two existing points + /// + /// See [specific implementations](#implementors) for details. + fn point_at_ratio_between(start: Point, end: Point, ratio_from_start: F) -> Point; + + // TODO: + // fn point_at_distance_between(start: Point, end: Point, distance_from_start: F) -> Point; + + /// Interpolates `Point`s along a line between `start` and `end`. + /// + /// See [specific implementations](#implementors) for details. + /// + /// 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? + fn points_along_line( + start: Point, + end: Point, + max_distance: F, + include_ends: bool, + ) -> impl Iterator>; +} diff --git a/geo/src/algorithm/line_measures/metric_spaces/euclidean.rs b/geo/src/algorithm/line_measures/metric_spaces/euclidean.rs new file mode 100644 index 0000000000..a3e292fd33 --- /dev/null +++ b/geo/src/algorithm/line_measures/metric_spaces/euclidean.rs @@ -0,0 +1,71 @@ +use super::super::Distance; +use crate::{GeoFloat, Point}; + +/// Euclidean space measures distance with the pythagorean formula - what you'd measure with a ruler. +/// +/// You must use projected coordinates with Euclidean space — +/// for lon/lat points, use the [`Haversine`], [`Geodesic`], or other [metric spaces]. +/// +/// [`Haversine`]: super::Haversine +/// [`Geodesic`]: super::Geodesic +/// [metric spaces]: super +pub struct Euclidean; + +/// Calculate the Euclidean distance (a.k.a. pythagorean distance) between two Points +impl Distance, Point> for Euclidean { + /// Calculate the Euclidean distance (a.k.a. pythagorean distance) between two Points + /// + /// # Units + /// - `origin`, `destination`: Point where the units of x/y represent non-angular units + /// — e.g. meters or miles, not lon/lat. For lon/lat points, use the + /// [`Haversine`] or [`Geodesic`] [metric spaces]. + /// - returns: distance in the same units as the `origin` and `destination` points + /// + /// # Example + /// ``` + /// use geo::{Euclidean, Distance}; + /// use geo::Point; + /// // web mercator + /// let new_york_city = Point::new(-8238310.24, 4942194.78); + /// // web mercator + /// let london = Point::new(-14226.63, 6678077.70); + /// let distance: f64 = Euclidean::distance(new_york_city, london); + /// + /// assert_eq!( + /// 8_405_286., // meters in web mercator + /// distance.round() + /// ); + /// ``` + /// + /// [`Haversine`]: super::Haversine + /// [`Geodesic`]: super::Geodesic + /// [metric spaces]: super + fn distance(origin: Point, destination: Point) -> F { + crate::EuclideanDistance::euclidean_distance(&origin, &destination) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + type MetricSpace = Euclidean; + + mod distance { + use super::*; + + #[test] + fn new_york_to_london() { + // web mercator + let new_york_city = Point::new(-8238310.24, 4942194.78); + // web mercator + let london = Point::new(-14226.63, 6678077.70); + let distance: f64 = MetricSpace::distance(new_york_city, london); + + assert_relative_eq!( + 8_405_286., // meters in web mercator + distance.round() + ); + } + } +} diff --git a/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs b/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs new file mode 100644 index 0000000000..63af33c564 --- /dev/null +++ b/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs @@ -0,0 +1,325 @@ +use super::super::{Bearing, Destination, Distance, InterpolatePoint}; +use crate::Point; + +/// An ellipsoidal model of the earth, using methods given by [Karney (2013)]. +/// +/// Distances are computed using [geodesic lines]. +/// +/// [geodesic lines]: https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid +/// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf +pub struct Geodesic; + +impl Bearing for Geodesic { + /// Returns the bearing from `origin` to `destination` in degrees along a [geodesic line]. + /// + /// # Units + /// + /// - `origin`, `destination`: Point where x/y are lon/lat degree coordinates + /// - returns: degrees, where: North: 0°, East: 90°, South: 180°, West: 270° + /// + /// ``` + /// # use approx::assert_relative_eq; + /// use geo::{Geodesic, Bearing}; + /// use geo::Point; + /// + /// let origin = Point::new(9.0, 10.0); + /// let destination = Point::new(9.5, 10.1); + /// let bearing = Geodesic::bearing(origin, destination); + /// // A little north of east + /// assert_relative_eq!(bearing, 78.54, 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 bearing(origin: Point, destination: Point) -> f64 { + (crate::algorithm::GeodesicBearing::geodesic_bearing(&origin, destination) + 360.0) % 360.0 + } +} + +impl Destination for Geodesic { + /// Returns a new point having travelled the `distance` along a [geodesic line] + /// from the `origin` point with the given `bearing`. + /// + /// This uses the geodesic methods given by [Karney (2013)]. + /// + /// # Units + /// + /// - `bearing`: degrees, where: North: 0°, East: 90°, South: 180°, West: 270° + /// - `distance`: meters + /// - returns: Point where x/y are lon/lat degree coordinates + /// + /// # Examples + /// + /// ``` + /// # use approx::assert_relative_eq; + /// use geo::{Geodesic, Destination}; + /// use geo::Point; + /// + /// // Determine the point 100 km NE of JFK airport. + /// let jfk = Point::new(-73.78, 40.64); + /// let northeast_bearing = 45.0; + /// let distance = 100_000.0; + /// + /// let northeast_of_jfk = Geodesic::destination(jfk, northeast_bearing, distance); + /// assert_relative_eq!(Point::new(-72.94, 41.27), northeast_of_jfk, 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 destination(origin: Point, bearing: f64, distance: f64) -> Point { + crate::algorithm::GeodesicDestination::geodesic_destination(&origin, bearing, distance) + } +} + +impl Distance, Point> for Geodesic { + /// Determine the length of the [geodesic line] between two geometries on an ellipsoidal model of the earth. + /// + /// # Units + /// - `origin`, `destination`: Point where x/y are lon/lat degree coordinates/ + /// - returns: meters + /// + /// # Examples + /// ```rust + /// use geo::{Geodesic, Distance}; + /// use geo::Point; + /// + /// // New York City + /// let new_york_city = Point::new(-74.006, 40.7128); + /// + /// // London + /// let london = Point::new(-0.1278, 51.5074); + /// + /// let distance = Geodesic::distance(new_york_city, london); + /// + /// assert_eq!( + /// 5_585_234., // meters + /// distance.round() + /// ); + /// ``` + /// + /// # 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 distance(origin: Point, destination: Point) -> f64 { + crate::algorithm::GeodesicDistance::geodesic_distance(&origin, &destination) + } +} + +/// Interpolate Point(s) along a [geodesic line]. +/// +/// [geodesic line]: https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid +impl InterpolatePoint for Geodesic { + /// Returns a new Point along a [geodesic line] between two existing points on an ellipsoidal model of the earth. + /// + /// # 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_ratio_between(p1, p2, 0.1); + /// assert_relative_eq!(closer_to_p1, Point::new(19.52, 25.31), epsilon = 1.0e-2); + /// + /// let closer_to_p2 = Geodesic::point_at_ratio_between(p1, p2, 0.9); + /// assert_relative_eq!(closer_to_p2, Point::new(114.73, 29.69), epsilon = 1.0e-2); + /// + /// let midpoint = Geodesic::point_at_ratio_between(p1, p2, 0.5); + /// assert_relative_eq!(midpoint, Point::new(65.88, 37.72), 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_ratio_between( + start: Point, + end: Point, + ratio_from_start: f64, + ) -> Point { + crate::algorithm::GeodesicIntermediate::geodesic_intermediate( + &start, + &end, + ratio_from_start, + ) + } + + /// Interpolates `Point`s along a [geodesic line] between `start` and `end`. + /// + /// As many points as necessary will be added such that the geodesic 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? + /// + /// # 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 points_along_line( + start: Point, + end: Point, + max_distance: f64, + include_ends: bool, + ) -> impl Iterator> { + crate::algorithm::GeodesicIntermediate::geodesic_intermediate_fill( + &start, + &end, + max_distance, + include_ends, + ) + .into_iter() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + type MetricSpace = Geodesic; + + mod bearing { + use super::*; + + #[test] + fn north() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(0.0, 1.0); + assert_relative_eq!(0.0, MetricSpace::bearing(origin, destination)); + } + + #[test] + fn east() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(1.0, 0.0); + assert_relative_eq!(90.0, MetricSpace::bearing(origin, destination)); + } + + #[test] + fn south() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(0.0, -1.0); + assert_relative_eq!(180.0, MetricSpace::bearing(origin, destination)); + } + + #[test] + fn west() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(-1.0, 0.0); + assert_relative_eq!(270.0, MetricSpace::bearing(origin, destination)); + } + } + + mod destination { + use super::*; + + #[test] + fn north() { + let origin = Point::new(0.0, 0.0); + let bearing = 0.0; + assert_relative_eq!( + Point::new(0.0, 0.9043687229127633), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + #[test] + fn east() { + let origin = Point::new(0.0, 0.0); + let bearing = 90.0; + assert_relative_eq!( + Point::new(0.8983152841195217, 0.0), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + #[test] + fn south() { + let origin = Point::new(0.0, 0.0); + let bearing = 180.0; + assert_relative_eq!( + Point::new(0.0, -0.9043687229127633), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + #[test] + fn west() { + let origin = Point::new(0.0, 0.0); + let bearing = 270.0; + assert_relative_eq!( + Point::new(-0.8983152841195217, 0.0), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + mod distance { + use super::*; + + #[test] + fn new_york_to_london() { + let new_york_city = Point::new(-74.006f64, 40.7128f64); + let london = Point::new(-0.1278f64, 51.5074f64); + + let distance = MetricSpace::distance(new_york_city, london); + + assert_relative_eq!( + 5_585_234.0, // meters + distance.round() + ); + } + } + + mod interpolate_point { + use super::*; + + #[test] + fn point_at_ratio_between_midpoint() { + let start = Point::new(10.0, 20.0); + let end = Point::new(125.0, 25.0); + 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); + let end = Point::new(125.0, 25.0); + let max_dist = 1000000.0; // meters + let route = + MetricSpace::points_along_line(start, end, max_dist, true).collect::>(); + assert_eq!(route.len(), 13); + assert_eq!(route[0], start); + 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); + let end = Point::new(125.0, 25.0); + let max_dist = 1000000.0; // meters + let route = + MetricSpace::points_along_line(start, end, max_dist, false).collect::>(); + assert_eq!(route.len(), 11); + assert_eq!(route[0], Point::new(17.878754355562464, 24.466667836189565)); + } + } + } +} diff --git a/geo/src/algorithm/line_measures/metric_spaces/haversine.rs b/geo/src/algorithm/line_measures/metric_spaces/haversine.rs new file mode 100644 index 0000000000..e394e89422 --- /dev/null +++ b/geo/src/algorithm/line_measures/metric_spaces/haversine.rs @@ -0,0 +1,322 @@ +use num_traits::FromPrimitive; + +use super::super::{Bearing, Destination, Distance, InterpolatePoint}; +use crate::{CoordFloat, Point}; + +/// A spherical model of the earth using the [haversine formula]. +/// +/// Distances are considered [great circle] lengths. +/// +/// # References +/// +/// *Note*: this implementation uses a mean earth radius of 6371.088 km, based on the [recommendation of +/// the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf) +/// +/// [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula// +/// [great circle]: https://en.wikipedia.org/wiki/Great_circle +pub struct Haversine; + +impl Bearing for Haversine { + /// Returns the bearing from `origin` to `destination` in degrees along a [great circle]. + /// + /// # Units + /// + /// - `origin`, `destination`: Points where x/y are lon/lat degree coordinates + /// returns: degrees, where: North: 0°, East: 90°, South: 180°, West: 270° + /// + /// # Examples + /// + /// ``` + /// # use approx::assert_relative_eq; + /// use geo::{Haversine, Bearing}; + /// use geo::Point; + /// + /// let origin = Point::new(9.0, 10.0); + /// let destination = Point::new(9.5, 10.1); + /// let bearing = Haversine::bearing(origin, destination); + /// // A little north of east + /// assert_relative_eq!(bearing, 78.47, epsilon = 1.0e-2); + /// ``` + /// + /// # References + /// + /// Bullock, R.: Great Circle Distances and Bearings Between Two Locations, 2007. + /// () + /// + /// [great circle]: https://en.wikipedia.org/wiki/Great_circle + fn bearing(origin: Point, destination: Point) -> F { + let three_sixty = + F::from(360.0).expect("Numeric type to be constructable from primitive 360"); + (crate::algorithm::HaversineBearing::haversine_bearing(&origin, destination) + three_sixty) + % three_sixty + } +} + +impl Destination for Haversine { + /// Returns a new point having travelled the `distance` along a [great circle] + /// from the `origin` point with the given `bearing`. + /// + /// # Units + /// + /// - `origin`: Point where x/y are lon/lat degree coordinates + /// - `bearing`: degrees, where: North: 0°, East: 90°, South: 180°, West: 270° + /// - `distance`: meters + /// - returns: Point where x/y are lon/lat degree coordinates + /// + /// # Examples + /// + /// ``` + /// # use approx::assert_relative_eq; + /// use geo::{Haversine, Destination}; + /// use geo::Point; + /// + /// let origin = Point::new(9.177789688110352, 48.776781529534965); + /// let destination = Haversine::destination(origin, 45., 10000.); + /// assert_relative_eq!(Point::new(9.274409949623532, 48.84033274015048), destination); + /// ``` + /// + /// # References + /// + /// *Note*: this implementation uses a mean earth radius of 6371.088 km, based on the [recommendation of + /// the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf) + /// + /// [great circle]: https://en.wikipedia.org/wiki/Great_circle + fn destination(origin: Point, bearing: F, distance: F) -> Point { + crate::algorithm::HaversineDestination::haversine_destination(&origin, bearing, distance) + } +} + +impl Distance, Point> for Haversine { + /// Determine the distance between two points using the [haversine formula]. + /// + /// # Units + /// + /// - `origin`, `destination`: Points where x/y are lon/lat degree coordinates + /// - returns: meters + /// + /// # Examples + /// + /// ``` + /// # use approx::assert_relative_eq; + /// use geo::{Haversine, Distance}; + /// use geo::Point; + /// + /// let new_york_city = Point::new(-74.006f64, 40.7128f64); + /// let london = Point::new(-0.1278f64, 51.5074f64); + /// + /// let distance = Haversine::distance(new_york_city, london); + /// + /// assert_relative_eq!( + /// 5_570_230., // meters + /// distance.round() + /// ); + /// ``` + /// + /// # References + /// + /// *Note*: this implementation uses a mean earth radius of 6371.088 km, based on the [recommendation of + /// the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf) + /// + /// [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula + fn distance(origin: Point, destination: Point) -> F { + crate::algorithm::HaversineDistance::haversine_distance(&origin, &destination) + } +} + +/// Interpolate Point(s) along a [great circle]. +/// +/// [great circle]: https://en.wikipedia.org/wiki/Great_circle +impl InterpolatePoint 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_ratio_between(p1, p2, 0.1); + /// assert_relative_eq!(closer_to_p1, Point::new(19.52, 25.27), epsilon = 1.0e-2); + /// + /// let closer_to_p2 = Haversine::point_at_ratio_between(p1, p2, 0.9); + /// assert_relative_eq!(closer_to_p2, Point::new(114.72, 29.65), epsilon = 1.0e-2); + /// + /// let midpoint = Haversine::point_at_ratio_between(p1, p2, 0.5); + /// assert_relative_eq!(midpoint, Point::new(65.87, 37.62), epsilon = 1.0e-2); + /// ``` + /// + /// [great circle]: https://en.wikipedia.org/wiki/Great_circle + fn point_at_ratio_between( + start: Point, + end: Point, + ratio_from_start: f64, + ) -> Point { + crate::algorithm::HaversineIntermediate::haversine_intermediate( + &start, + &end, + ratio_from_start, + ) + } + + /// Interpolates `Point`s along a [great circle] between `start` and `end`. + /// + /// As many points as necessary will be added such that the [haversine 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? + /// + /// [great circle]: https://en.wikipedia.org/wiki/Great_circle + /// [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula + fn points_along_line( + start: Point, + end: Point, + max_distance: f64, + include_ends: bool, + ) -> impl Iterator> { + crate::algorithm::HaversineIntermediate::haversine_intermediate_fill( + &start, + &end, + max_distance, + include_ends, + ) + .into_iter() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + type MetricSpace = Haversine; + + mod bearing { + use super::*; + + #[test] + fn north() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(0.0, 1.0); + assert_relative_eq!(0.0, MetricSpace::bearing(origin, destination)); + } + + #[test] + fn east() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(1.0, 0.0); + assert_relative_eq!(90.0, MetricSpace::bearing(origin, destination)); + } + + #[test] + fn south() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(0.0, -1.0); + assert_relative_eq!(180.0, MetricSpace::bearing(origin, destination)); + } + + #[test] + fn west() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(-1.0, 0.0); + assert_relative_eq!(270.0, MetricSpace::bearing(origin, destination)); + } + } + + mod destination { + use super::*; + + #[test] + fn north() { + let origin = Point::new(0.0, 0.0); + let bearing = 0.0; + assert_relative_eq!( + Point::new(0.0, 0.899320363724538), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + #[test] + fn east() { + let origin = Point::new(0.0, 0.0); + let bearing = 90.0; + assert_relative_eq!( + Point::new(0.8993203637245415, 5.506522912913066e-17), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + #[test] + fn south() { + let origin = Point::new(0.0, 0.0); + let bearing = 180.0; + assert_relative_eq!( + Point::new(0.0, -0.899320363724538), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + #[test] + fn west() { + let origin = Point::new(0.0, 0.0); + let bearing = 270.0; + assert_relative_eq!( + Point::new(-0.8993203637245415, -1.6519568738739197e-16), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + } + + mod distance { + use super::*; + + #[test] + fn new_york_to_london() { + let new_york_city = Point::new(-74.006f64, 40.7128f64); + let london = Point::new(-0.1278f64, 51.5074f64); + + let distance = MetricSpace::distance(new_york_city, london); + + assert_relative_eq!( + 5_570_230., // meters + distance.round() + ); + } + } + mod interpolate_point { + use super::*; + + #[test] + fn point_at_ratio_between_midpoint() { + let start = Point::new(10.0, 20.0); + let end = Point::new(125.0, 25.0); + let midpoint = MetricSpace::point_at_ratio_between(start, end, 0.5); + assert_relative_eq!(midpoint, Point::new(65.87394172511485, 37.61809316888599)); + } + #[test] + fn points_along_line_with_endpoints() { + let start = Point::new(10.0, 20.0); + let end = Point::new(125.0, 25.0); + let max_dist = 1000000.0; // meters + let route = + MetricSpace::points_along_line(start, end, max_dist, true).collect::>(); + assert_eq!(route.len(), 13); + assert_eq!(route[0], start); + assert_eq!(route.last().unwrap(), &end); + assert_relative_eq!(route[1], Point::new(17.882467331860965, 24.435542998803793)); + } + #[test] + fn points_along_line_without_endpoints() { + let start = Point::new(10.0, 20.0); + let end = Point::new(125.0, 25.0); + let max_dist = 1000000.0; // meters + let route = + MetricSpace::points_along_line(start, end, max_dist, false).collect::>(); + assert_eq!(route.len(), 11); + assert_relative_eq!(route[0], Point::new(17.882467331860965, 24.435542998803793)); + } + } +} diff --git a/geo/src/algorithm/line_measures/metric_spaces/mod.rs b/geo/src/algorithm/line_measures/metric_spaces/mod.rs new file mode 100644 index 0000000000..7a1a153fe4 --- /dev/null +++ b/geo/src/algorithm/line_measures/metric_spaces/mod.rs @@ -0,0 +1,11 @@ +mod euclidean; +pub use euclidean::Euclidean; + +mod geodesic; +pub use geodesic::Geodesic; + +mod haversine; +pub use haversine::Haversine; + +mod rhumb; +pub use rhumb::Rhumb; diff --git a/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs b/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs new file mode 100644 index 0000000000..de12c93ca4 --- /dev/null +++ b/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs @@ -0,0 +1,309 @@ +use num_traits::FromPrimitive; + +use super::super::{Bearing, Destination, Distance, InterpolatePoint}; +use crate::{CoordFloat, Point}; + +/// Provides [rhumb line] (a.k.a. loxodrome) geometry operations. A rhumb line appears as a straight +/// line on a Mercator projection map. +/// +/// # References +/// +/// The distance, destination, and bearing implementations are adapted in part +/// from their equivalents in [Turf.js](https://turfjs.org/), which in turn are +/// adapted from the Movable Type +/// [spherical geodesy tools](https://www.movable-type.co.uk/scripts/latlong.html). +/// +/// Turf.js is copyright its authors and the geodesy tools are copyright Chris +/// Veness; both are available under an MIT license. +/// +/// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line +pub struct Rhumb; + +impl Bearing for Rhumb { + /// Returns the bearing from `origin` to `destination` in degrees along a [rhumb line]. + /// + /// # Units + /// + /// - `origin`, `destination`: Points where x/y are lon/lat degree coordinates + /// - returns: degrees, where: North: 0°, East: 90°, South: 180°, West: 270°/ + /// + /// # Examples + /// + /// ``` + /// # use approx::assert_relative_eq; + /// use geo::{Rhumb, Bearing}; + /// use geo::Point; + /// + /// let origin = Point::new(9.177789688110352, 48.776781529534965); + /// let destination = Point::new(9.274348757829898, 48.84037308229984); + /// let bearing = Rhumb::bearing(origin, destination); + /// assert_relative_eq!(bearing, 45., epsilon = 1.0e-6); + /// ``` + /// + /// # References + /// + /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line + /// + /// Bullock, R.: Great Circle Distances and Bearings Between Two Locations, 2007. + /// () + fn bearing(origin: Point, destination: Point) -> F { + crate::algorithm::RhumbBearing::rhumb_bearing(&origin, destination) + } +} + +impl Destination for Rhumb { + /// Returns a new point having travelled the `distance` along a [rhumb line] + /// from the `origin` point with the given `bearing`. + /// + /// # Units + /// + /// - `origin`: Point where x/y are lon/lat degree coordinates + /// - `bearing`: degrees, where: North: 0°, East: 90°, South: 180°, West: 270° + /// - `distance`: meters + /// - returns: Point where x/y are lon/lat degree coordinates + /// + /// # Examples + /// + /// ``` + /// # use approx::assert_relative_eq; + /// use geo::{Rhumb, Destination}; + /// use geo::Point; + /// + /// let p_1 = Point::new(9.177789688110352, 48.776781529534965); + /// let p_2 = Rhumb::destination(p_1, 45., 10000.); + /// assert_relative_eq!(p_2, Point::new(9.274348757829898, 48.84037308229984)) + /// ``` + /// + /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line + fn destination(origin: Point, bearing: F, distance: F) -> Point { + crate::algorithm::RhumbDestination::rhumb_destination(&origin, bearing, distance) + } +} + +impl Distance, Point> for Rhumb { + /// Determine the distance along the [rhumb line] between two points. + /// + /// # Units + /// + /// - `origin`, `destination`: Points where x/y are lon/lat degree coordinates + /// - returns: meters + /// + /// # Examples + /// + /// ``` + /// use geo::{Rhumb, Distance}; + /// use geo::point; + /// + /// // New York City + /// let p1 = point!(x: -74.006f64, y: 40.7128); + /// + /// // London + /// let p2 = point!(x: -0.1278, y: 51.5074); + /// + /// let distance = Rhumb::distance(p1, p2); + /// + /// assert_eq!( + /// 5_794_129., // meters + /// distance.round() + /// ); + /// ``` + /// + /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line + fn distance(origin: Point, destination: Point) -> F { + crate::algorithm::RhumbDistance::rhumb_distance(&origin, &destination) + } +} + +/// Interpolate Point(s) along a [rhumb line]. +/// +/// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line +impl InterpolatePoint 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_ratio_between(p1, p2, 0.1); + /// assert_relative_eq!(closer_to_p1, Point::new(21.32, 20.50), epsilon = 1.0e-2); + /// + /// let closer_to_p2 = Rhumb::point_at_ratio_between(p1, p2, 0.9); + /// assert_relative_eq!(closer_to_p2, Point::new(113.31, 24.50), epsilon = 1.0e-2); + /// + /// let midpoint = Rhumb::point_at_ratio_between(p1, p2, 0.5); + /// assert_relative_eq!(midpoint, Point::new(66.98, 22.50), epsilon = 1.0e-2); + /// ``` + /// + /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line + fn point_at_ratio_between( + start: Point, + end: Point, + ratio_from_start: f64, + ) -> Point { + crate::algorithm::RhumbIntermediate::rhumb_intermediate(&start, &end, ratio_from_start) + } + + /// Interpolates `Point`s along a [rhumb 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? + /// + /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line + fn points_along_line( + start: Point, + end: Point, + max_distance: f64, + include_ends: bool, + ) -> impl Iterator> { + crate::algorithm::RhumbIntermediate::rhumb_intermediate_fill( + &start, + &end, + max_distance, + include_ends, + ) + .into_iter() + } +} + +#[cfg(test)] +mod tests { + use super::*; + + type MetricSpace = Rhumb; + + mod bearing { + use super::*; + + #[test] + fn north() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(0.0, 1.0); + assert_relative_eq!(0.0, MetricSpace::bearing(origin, destination)); + } + + #[test] + fn east() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(1.0, 0.0); + assert_relative_eq!(90.0, MetricSpace::bearing(origin, destination)); + } + + #[test] + fn south() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(0.0, -1.0); + assert_relative_eq!(180.0, MetricSpace::bearing(origin, destination)); + } + + #[test] + fn west() { + let origin = Point::new(0.0, 0.0); + let destination = Point::new(-1.0, 0.0); + assert_relative_eq!(270.0, MetricSpace::bearing(origin, destination)); + } + } + + mod destination { + use super::*; + + #[test] + fn north() { + let origin = Point::new(0.0, 0.0); + let bearing = 0.0; + assert_relative_eq!( + Point::new(0.0, 0.899320363724538), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + #[test] + fn east() { + let origin = Point::new(0.0, 0.0); + let bearing = 90.0; + assert_relative_eq!( + Point::new(0.8993203637245415, 5.506522912913066e-17), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + #[test] + fn south() { + let origin = Point::new(0.0, 0.0); + let bearing = 180.0; + assert_relative_eq!( + Point::new(0.0, -0.899320363724538), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + + #[test] + fn west() { + let origin = Point::new(0.0, 0.0); + let bearing = 270.0; + assert_relative_eq!( + Point::new(-0.8993203637245415, -1.6520247072649334e-16), + MetricSpace::destination(origin, bearing, 100_000.0) + ); + } + } + + mod distance { + use super::*; + + #[test] + fn new_york_to_london() { + let new_york_city = Point::new(-74.006, 40.7128); + let london = Point::new(-0.1278, 51.5074); + + let distance: f64 = MetricSpace::distance(new_york_city, london); + + assert_relative_eq!( + 5_794_129., // meters + distance.round() + ); + } + } + + mod interpolate_point { + use super::*; + + #[test] + fn point_at_ratio_between_midpoint() { + let start = Point::new(10.0, 20.0); + let end = Point::new(125.0, 25.0); + let midpoint = MetricSpace::point_at_ratio_between(start, end, 0.5); + assert_relative_eq!(midpoint, Point::new(66.98011173721943, 22.500000000000007)); + } + #[test] + fn points_along_line_with_endpoints() { + let start = Point::new(10.0, 20.0); + let end = Point::new(125.0, 25.0); + let max_dist = 1000000.0; // meters + let route = + MetricSpace::points_along_line(start, end, max_dist, true).collect::>(); + assert_eq!(route.len(), 13); + assert_eq!(route[0], start); + assert_eq!(route.last().unwrap(), &end); + assert_relative_eq!(route[1], Point::new(19.43061818495096, 20.416666666666668)); + } + #[test] + fn points_along_line_without_endpoints() { + let start = Point::new(10.0, 20.0); + let end = Point::new(125.0, 25.0); + let max_dist = 1000000.0; // meters + let route = + MetricSpace::points_along_line(start, end, max_dist, false).collect::>(); + assert_eq!(route.len(), 11); + assert_relative_eq!(route[0], Point::new(19.43061818495096, 20.416666666666668)); + } + } +} diff --git a/geo/src/algorithm/line_measures/mod.rs b/geo/src/algorithm/line_measures/mod.rs new file mode 100644 index 0000000000..a38aa5e40d --- /dev/null +++ b/geo/src/algorithm/line_measures/mod.rs @@ -0,0 +1,16 @@ +//! Line measurements like [`Bearing`] and [`Distance`] for various metric spaces like [`Haversine`], [`Geodesic`], and [`Rhumb`]. + +mod bearing; +pub use bearing::Bearing; + +mod destination; +pub use destination::Destination; + +mod distance; +pub use distance::Distance; + +mod interpolate_point; +pub use interpolate_point::InterpolatePoint; + +pub mod metric_spaces; +pub use metric_spaces::{Geodesic, Haversine, Rhumb}; diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index 8eadc51e46..b0a6e4dd1f 100644 --- a/geo/src/algorithm/mod.rs +++ b/geo/src/algorithm/mod.rs @@ -6,16 +6,6 @@ pub use kernels::{Kernel, Orientation}; pub mod area; pub use area::Area; -/// Calculate the bearing to another `Point`, in degrees. -#[deprecated( - since = "0.24.1", - note = "renamed to `haversine_bearing::HaversineBearing`" -)] -pub mod bearing; -#[allow(deprecated)] -#[deprecated(since = "0.24.1", note = "renamed to `HaversineBearing`")] -pub use bearing::Bearing; - /// Boolean Ops such as union, xor, difference; pub mod bool_ops; pub use bool_ops::{BooleanOps, OpType}; @@ -295,6 +285,10 @@ pub use outlier_detection::OutlierDetection; pub mod monotone; pub use monotone::{monotone_subdivision, MonoPoly, MonotonicPolygons}; +pub mod line_measures; +pub use line_measures::metric_spaces::{Euclidean, Geodesic, Haversine, Rhumb}; +pub use line_measures::{Bearing, Destination, Distance, InterpolatePoint}; + /// Rhumb-line-related algorithms and utils pub mod rhumb; pub use rhumb::{RhumbBearing, RhumbDestination, RhumbDistance, RhumbIntermediate, RhumbLength}; diff --git a/geo/src/algorithm/rhumb/bearing.rs b/geo/src/algorithm/rhumb/bearing.rs index c1c6e3c853..4c3be743f8 100644 --- a/geo/src/algorithm/rhumb/bearing.rs +++ b/geo/src/algorithm/rhumb/bearing.rs @@ -25,7 +25,6 @@ pub trait RhumbBearing { /// assert_relative_eq!(bearing, 45., epsilon = 1.0e-6); /// ``` /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line - fn rhumb_bearing(&self, point: Point) -> T; } diff --git a/geo/src/lib.rs b/geo/src/lib.rs index f6ecfae665..f15d5bfd5f 100644 --- a/geo/src/lib.rs +++ b/geo/src/lib.rs @@ -220,6 +220,7 @@ extern crate serde; pub use crate::algorithm::*; pub use crate::types::Closest; +use num_traits::FromPrimitive; use std::cmp::Ordering; pub use crate::relate::PreparedGeometry; @@ -312,7 +313,7 @@ impl GeoFloat for T where } /// A trait for methods which work for both integers **and** floating point -pub trait GeoNum: CoordNum { +pub trait GeoNum: CoordNum + FromPrimitive { type Ker: Kernel; /// Return the ordering between self and other.