diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 951dfd3cb..bbc25701b 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -3,7 +3,7 @@ ## Unreleased - BREAKING: `Densify` and `Length` are now defined on the metric space, rather than a generic method on the geometry. - ``` + ```rust // before line_string.length::() line_string.densify::() @@ -12,8 +12,16 @@ Euclidean.length(&line_string) Euclidean.densify(&line_string) ``` -- Add `CustomHaversine` for doing calculations on spheres with a custom radius. +- `Haversine` can be configured with a custom radius for doing calculations on custom sphere. Use `HAVERSINE` for the default earth radius. - + ```rust + // before + Haversine::distance(point1, point2) + + // after + Haversine::new(3_389_500.0).distance(point1, point2) + HAVERSINE.distance(point_1, point_2) + ``` - Docs: Fix page location of citation for mean earth radius used in Haversine calculations - - Docs: Add top-level doc link for `InteriorPoint` diff --git a/geo/src/algorithm/cross_track_distance.rs b/geo/src/algorithm/cross_track_distance.rs index 0aac1f899..37d9da75c 100644 --- a/geo/src/algorithm/cross_track_distance.rs +++ b/geo/src/algorithm/cross_track_distance.rs @@ -1,4 +1,4 @@ -use crate::{Bearing, Distance, Haversine, MEAN_EARTH_RADIUS}; +use crate::{Bearing, Distance, HAVERSINE, MEAN_EARTH_RADIUS}; use geo_types::{CoordFloat, Point}; use num_traits::FromPrimitive; @@ -43,9 +43,9 @@ where { fn cross_track_distance(&self, line_point_a: &Point, line_point_b: &Point) -> T { let mean_earth_radius = T::from(MEAN_EARTH_RADIUS).unwrap(); - let l_delta_13: T = Haversine.distance(*line_point_a, *self) / mean_earth_radius; - let theta_13: T = Haversine.bearing(*line_point_a, *self).to_radians(); - let theta_12: T = Haversine.bearing(*line_point_a, *line_point_b).to_radians(); + let l_delta_13: T = HAVERSINE.distance(*line_point_a, *self) / mean_earth_radius; + let theta_13: T = HAVERSINE.bearing(*line_point_a, *self).to_radians(); + let theta_12: T = HAVERSINE.bearing(*line_point_a, *line_point_b).to_radians(); let l_delta_xt: T = (l_delta_13.sin() * (theta_12 - theta_13).sin()).asin(); mean_earth_radius * l_delta_xt.abs() } @@ -55,7 +55,7 @@ where mod test { use crate::CrossTrackDistance; use crate::Point; - use crate::{Distance, Haversine}; + use crate::{Distance, HAVERSINE}; #[test] fn distance1_test() { @@ -90,13 +90,13 @@ mod test { assert_relative_eq!( p.cross_track_distance(&line_point_a, &line_point_b), - Haversine.distance(p, Point::new(1., 0.)), + HAVERSINE.distance(p, Point::new(1., 0.)), epsilon = 1.0e-6 ); assert_relative_eq!( p.cross_track_distance(&line_point_b, &line_point_a), - Haversine.distance(p, Point::new(1., 0.)), + HAVERSINE.distance(p, Point::new(1., 0.)), epsilon = 1.0e-6 ); } diff --git a/geo/src/algorithm/densify_haversine.rs b/geo/src/algorithm/densify_haversine.rs index 03ad9ff7e..f2689d59d 100644 --- a/geo/src/algorithm/densify_haversine.rs +++ b/geo/src/algorithm/densify_haversine.rs @@ -1,6 +1,6 @@ use num_traits::FromPrimitive; -use crate::line_measures::Haversine; +use crate::line_measures::HAVERSINE; // Densify will soon be deprecated too, so let's just allow deprecated for now #[allow(deprecated)] use crate::HaversineLength; @@ -10,7 +10,7 @@ use crate::{ #[deprecated( since = "0.29.0", - note = "Please use the `Haversine.densify(&line)` via the `Densify` trait instead." + note = "Please use the `HAVERSINE.densify(&line)` via the `Densify` trait instead." )] /// Returns a new spherical geometry containing both existing and new interpolated coordinates with /// a maximum distance of `max_distance` between them. @@ -50,7 +50,7 @@ where type Output = MultiPolygon; fn densify_haversine(&self, max_distance: T) -> Self::Output { - Haversine.densify(self, max_distance) + HAVERSINE.densify(self, max_distance) } } @@ -64,7 +64,7 @@ where type Output = Polygon; fn densify_haversine(&self, max_distance: T) -> Self::Output { - Haversine.densify(self, max_distance) + HAVERSINE.densify(self, max_distance) } } @@ -78,7 +78,7 @@ where type Output = MultiLineString; fn densify_haversine(&self, max_distance: T) -> Self::Output { - Haversine.densify(self, max_distance) + HAVERSINE.densify(self, max_distance) } } @@ -92,7 +92,7 @@ where type Output = LineString; fn densify_haversine(&self, max_distance: T) -> Self::Output { - Haversine.densify(self, max_distance) + HAVERSINE.densify(self, max_distance) } } @@ -106,7 +106,7 @@ where type Output = LineString; fn densify_haversine(&self, max_distance: T) -> Self::Output { - Haversine.densify(self, max_distance) + HAVERSINE.densify(self, max_distance) } } @@ -120,7 +120,7 @@ where type Output = Polygon; fn densify_haversine(&self, max_distance: T) -> Self::Output { - Haversine.densify(self, max_distance) + HAVERSINE.densify(self, max_distance) } } @@ -134,7 +134,7 @@ where type Output = Polygon; fn densify_haversine(&self, max_distance: T) -> Self::Output { - Haversine.densify(self, max_distance) + HAVERSINE.densify(self, max_distance) } } diff --git a/geo/src/algorithm/haversine_bearing.rs b/geo/src/algorithm/haversine_bearing.rs index 3549553a3..23db4056c 100644 --- a/geo/src/algorithm/haversine_bearing.rs +++ b/geo/src/algorithm/haversine_bearing.rs @@ -2,7 +2,7 @@ use crate::{CoordFloat, Point}; #[deprecated( since = "0.29.0", - note = "Please use the `Haversine.bearing` method from the `Bearing` trait instead" + note = "Please use the `HAVERSINE.bearing` method from the `Bearing` trait instead" )] /// Returns the bearing to another Point in degrees. /// diff --git a/geo/src/algorithm/haversine_closest_point.rs b/geo/src/algorithm/haversine_closest_point.rs index 60b7fa31d..e1afb145f 100644 --- a/geo/src/algorithm/haversine_closest_point.rs +++ b/geo/src/algorithm/haversine_closest_point.rs @@ -1,4 +1,4 @@ -use crate::line_measures::{Bearing, Destination, Distance, Haversine}; +use crate::line_measures::{Bearing, Destination, Distance, HAVERSINE}; use crate::{Closest, Contains}; use crate::{CoordsIter, GeoFloat, Point, MEAN_EARTH_RADIUS}; use geo_types::{ @@ -92,7 +92,7 @@ where } // This can probably be done cheaper - let d3 = Haversine.distance(p2, p1); + let d3 = HAVERSINE.distance(p2, p1); if d3 <= T::epsilon() { // I think here it should be return Closest::SinglePoint(p1) // If the line segment is degenerated to a point, that point is still the closest @@ -102,18 +102,18 @@ where } let pi = T::from(std::f64::consts::PI).unwrap(); - let crs_ad = Haversine.bearing(p1, *from).to_radians(); - let crs_ab = Haversine.bearing(p1, p2).to_radians(); + let crs_ad = HAVERSINE.bearing(p1, *from).to_radians(); + let crs_ab = HAVERSINE.bearing(p1, p2).to_radians(); let crs_ba = if crs_ab > T::zero() { crs_ab - pi } else { crs_ab + pi }; - let crs_bd = Haversine.bearing(p2, *from).to_radians(); + let crs_bd = HAVERSINE.bearing(p2, *from).to_radians(); let d_crs1 = crs_ad - crs_ab; let d_crs2 = crs_bd - crs_ba; - let d1 = Haversine.distance(p1, *from); + let d1 = HAVERSINE.distance(p1, *from); // d1, d2, d3 are in principle not needed, only the sign matters let projection1 = d_crs1.cos(); @@ -127,13 +127,13 @@ where if xtd < T::epsilon() { return Closest::Intersection(*from); } else { - return Closest::SinglePoint(Haversine.destination(p1, crs_ab.to_degrees(), atd)); + return Closest::SinglePoint(HAVERSINE.destination(p1, crs_ab.to_degrees(), atd)); } } // Projected falls outside the GC Arc // Return shortest distance pt, project either on point sp1 or sp2 - let d2 = Haversine.distance(p2, *from); + let d2 = HAVERSINE.distance(p2, *from); if d1 < d2 { return Closest::SinglePoint(p1); } @@ -166,7 +166,7 @@ where return intersect; } Closest::SinglePoint(pt) => { - let dist = Haversine.distance(pt, *from); + let dist = HAVERSINE.distance(pt, *from); if dist < min_distance { min_distance = dist; rv = Closest::SinglePoint(pt); @@ -198,7 +198,7 @@ where return (intersect, T::zero()); } Closest::SinglePoint(pt) => { - let dist = Haversine.distance(pt, *from); + let dist = HAVERSINE.distance(pt, *from); if dist < min_distance { min_distance = dist; rv = Closest::SinglePoint(pt); @@ -301,7 +301,7 @@ where // This mean on top of the line. Closest::Intersection(pt) => return Closest::Intersection(pt), Closest::SinglePoint(pt) => { - let dist = Haversine.distance(pt, *from); + let dist = HAVERSINE.distance(pt, *from); if dist < min_distance { min_distance = dist; rv = Closest::SinglePoint(pt); diff --git a/geo/src/algorithm/haversine_destination.rs b/geo/src/algorithm/haversine_destination.rs index f915569d0..fd48eeee2 100644 --- a/geo/src/algorithm/haversine_destination.rs +++ b/geo/src/algorithm/haversine_destination.rs @@ -1,9 +1,9 @@ -use crate::{CoordFloat, Destination, Haversine, Point}; +use crate::{CoordFloat, Destination, Point, HAVERSINE}; use num_traits::FromPrimitive; #[deprecated( since = "0.29.0", - note = "Please use the `Haversine.destination` method from the `Destination` trait instead" + note = "Please use the `HAVERSINE.destination` method from the `Destination` trait instead" )] /// Returns a new Point using the distance to the existing Point and a bearing for the direction /// @@ -39,7 +39,7 @@ where T: CoordFloat + FromPrimitive, { fn haversine_destination(&self, bearing: T, distance: T) -> Point { - Haversine.destination(*self, bearing, distance) + HAVERSINE.destination(*self, bearing, distance) } } diff --git a/geo/src/algorithm/haversine_distance.rs b/geo/src/algorithm/haversine_distance.rs index df1cf7a40..9ccc3f890 100644 --- a/geo/src/algorithm/haversine_distance.rs +++ b/geo/src/algorithm/haversine_distance.rs @@ -1,9 +1,9 @@ -use crate::{CoordFloat, Distance, Haversine, Point}; +use crate::{CoordFloat, Distance, Point, HAVERSINE}; use num_traits::FromPrimitive; #[deprecated( since = "0.29.0", - note = "Please use the `Haversine.distance` method from the `Distance` trait instead" + note = "Please use the `HAVERSINE.distance` method from the `Distance` trait instead" )] /// Determine the distance between two geometries using the [haversine formula]. /// @@ -55,7 +55,7 @@ where T: CoordFloat + FromPrimitive, { fn haversine_distance(&self, rhs: &Point) -> T { - Haversine.distance(*self, *rhs) + HAVERSINE.distance(*self, *rhs) } } diff --git a/geo/src/algorithm/haversine_intermediate.rs b/geo/src/algorithm/haversine_intermediate.rs index fce21e8c2..f2410752b 100644 --- a/geo/src/algorithm/haversine_intermediate.rs +++ b/geo/src/algorithm/haversine_intermediate.rs @@ -1,4 +1,4 @@ -use crate::{CoordFloat, Haversine, InterpolatePoint, Point}; +use crate::{CoordFloat, InterpolatePoint, Point, HAVERSINE}; use num_traits::FromPrimitive; #[deprecated( @@ -9,7 +9,7 @@ use num_traits::FromPrimitive; pub trait HaversineIntermediate { #[deprecated( since = "0.29.0", - note = "Please use `Haversine.point_at_ratio_between` from the `InterpolatePoint` trait instead" + note = "Please use `HAVERSINE.point_at_ratio_between` from the `InterpolatePoint` trait instead" )] /// Returns a new `Point` along a great circle route between `self` and `other`. /// @@ -40,7 +40,7 @@ pub trait HaversineIntermediate { #[deprecated( since = "0.29.0", - note = "Please use `Haversine.points_along_line` from the `InterpolatePoint` trait instead" + note = "Please use `HAVERSINE.points_along_line` from the `InterpolatePoint` trait instead" )] /// Interpolates `Point`s along a great circle route between self and `other`. /// @@ -62,7 +62,7 @@ where T: CoordFloat + FromPrimitive, { fn haversine_intermediate(&self, other: &Point, ratio: T) -> Point { - Haversine.point_at_ratio_between(*self, *other, ratio) + HAVERSINE.point_at_ratio_between(*self, *other, ratio) } fn haversine_intermediate_fill( @@ -71,7 +71,7 @@ where max_dist: T, include_ends: bool, ) -> Vec> { - Haversine + HAVERSINE .points_along_line(*self, *other, max_dist, include_ends) .collect() } diff --git a/geo/src/algorithm/haversine_length.rs b/geo/src/algorithm/haversine_length.rs index bac48e593..e999a25c5 100644 --- a/geo/src/algorithm/haversine_length.rs +++ b/geo/src/algorithm/haversine_length.rs @@ -1,11 +1,11 @@ use num_traits::FromPrimitive; use crate::{CoordFloat, Line, LineString, MultiLineString}; -use crate::{Haversine, Length}; +use crate::{Length, HAVERSINE}; #[deprecated( since = "0.29.0", - note = "Please use the `Haversine.length(&line)` via the `Length` trait instead." + note = "Please use the `HAVERSINE.length(&line)` via the `Length` trait instead." )] /// Determine the length of a geometry using the [haversine formula]. /// @@ -51,7 +51,7 @@ where T: CoordFloat + FromPrimitive, { fn haversine_length(&self) -> T { - Haversine.length(self) + HAVERSINE.length(self) } } @@ -61,7 +61,7 @@ where T: CoordFloat + FromPrimitive, { fn haversine_length(&self) -> T { - Haversine.length(self) + HAVERSINE.length(self) } } @@ -71,6 +71,6 @@ where T: CoordFloat + FromPrimitive, { fn haversine_length(&self) -> T { - Haversine.length(self) + HAVERSINE.length(self) } } diff --git a/geo/src/algorithm/line_measures/densify.rs b/geo/src/algorithm/line_measures/densify.rs index ca57a65cd..60b99b3d7 100644 --- a/geo/src/algorithm/line_measures/densify.rs +++ b/geo/src/algorithm/line_measures/densify.rs @@ -43,7 +43,7 @@ use num_traits::FromPrimitive; /// /// // For Haversine, the unit of distance is in meters /// let max_dist = 200_000.0; -/// let densified = Haversine.densify(&line_string, max_dist); +/// let densified = HAVERSINE.densify(&line_string, max_dist); /// // Haversine interprets coordinate points as lng/lat /// let expected_output = wkt!(LINESTRING( /// 0.0 0.0, @@ -319,7 +319,7 @@ impl Densifiable for Triangle { #[cfg(test)] mod tests { use super::*; - use crate::{coord, polygon, wkt, Euclidean, Geodesic, Haversine, Rhumb}; + use crate::{coord, polygon, wkt, Euclidean, Geodesic, Rhumb, HAVERSINE}; #[test] fn densify_line() { @@ -335,7 +335,7 @@ mod tests { let densified_rhumb = Rhumb.densify(&line, 100_000.0); assert!(densified_rhumb.coords_count() > 2); - let densified_haversine = Haversine.densify(&line, 100_000.0); + let densified_haversine = HAVERSINE.densify(&line, 100_000.0); assert!(densified_haversine.coords_count() > 2); } @@ -353,7 +353,7 @@ mod tests { let densified_rhumb_ls = Rhumb.densify(&line_string, 500_000.0); assert!(densified_rhumb_ls.coords_count() > line_string.coords_count()); - let densified_haversine_ls = Haversine.densify(&line_string, 500_000.0); + let densified_haversine_ls = HAVERSINE.densify(&line_string, 500_000.0); assert!(densified_haversine_ls.coords_count() > line_string.coords_count()); } @@ -459,7 +459,7 @@ mod tests { 4.925 45.804 ))); - let actual_haversine = Haversine.densify(&polygon, 50000.0); + let actual_haversine = HAVERSINE.densify(&polygon, 50000.0); assert_relative_eq!(actual_haversine, exepcted_haversine); let expected_geodesic = wkt!(POLYGON(( @@ -504,7 +504,7 @@ mod tests { -3.1944 55.949 )); - let dense = Haversine.densify(&linestring, 110.0); + let dense = HAVERSINE.densify(&linestring, 110.0); assert_relative_eq!(dense, expected); } @@ -512,7 +512,7 @@ mod tests { fn test_line_densify() { let output = wkt!(LINESTRING(0.0 0.0, 0.0 0.5, 0.0 1.0)); let line = Line::new(coord! {x: 0.0, y: 0.0}, coord! { x: 0.0, y: 1.0 }); - let dense = Haversine.densify(&line, 100000.0); + let dense = HAVERSINE.densify(&line, 100000.0); assert_relative_eq!(dense, output); } } diff --git a/geo/src/algorithm/line_measures/interpolate_point.rs b/geo/src/algorithm/line_measures/interpolate_point.rs index c60b8c6ba..6e0d1d5f4 100644 --- a/geo/src/algorithm/line_measures/interpolate_point.rs +++ b/geo/src/algorithm/line_measures/interpolate_point.rs @@ -42,7 +42,7 @@ pub trait InterpolatePoint { #[cfg(test)] mod tests { - use crate::{Euclidean, Geodesic, Haversine, InterpolatePoint, Point, Rhumb}; + use crate::{Euclidean, Geodesic, InterpolatePoint, Point, Rhumb, HAVERSINE}; #[test] fn point_at_ratio_between_line_ends() { @@ -50,13 +50,13 @@ mod tests { let end = Point::new(1.0, 1.0); let ratio = 0.0; - assert_eq!(Haversine.point_at_ratio_between(start, end, ratio), start); + assert_eq!(HAVERSINE.point_at_ratio_between(start, end, ratio), start); assert_eq!(Euclidean.point_at_ratio_between(start, end, ratio), start); assert_eq!(Geodesic.point_at_ratio_between(start, end, ratio), start); assert_eq!(Rhumb.point_at_ratio_between(start, end, ratio), start); let ratio = 1.0; - assert_eq!(Haversine.point_at_ratio_between(start, end, ratio), end); + assert_eq!(HAVERSINE.point_at_ratio_between(start, end, ratio), end); assert_eq!(Euclidean.point_at_ratio_between(start, end, ratio), end); assert_eq!(Geodesic.point_at_ratio_between(start, end, ratio), end); assert_eq!(Rhumb.point_at_ratio_between(start, end, ratio), end); @@ -70,19 +70,19 @@ mod tests { let start = Point::new(1.0, 1.0); let ratio = 0.0; - assert_eq!(Haversine.point_at_ratio_between(start, start, ratio), start); + assert_eq!(HAVERSINE.point_at_ratio_between(start, start, ratio), start); assert_eq!(Euclidean.point_at_ratio_between(start, start, ratio), start); assert_eq!(Geodesic.point_at_ratio_between(start, start, ratio), start); assert_eq!(Rhumb.point_at_ratio_between(start, start, ratio), start); let ratio = 0.5; - assert_eq!(Haversine.point_at_ratio_between(start, start, ratio), start); + assert_eq!(HAVERSINE.point_at_ratio_between(start, start, ratio), start); assert_eq!(Euclidean.point_at_ratio_between(start, start, ratio), start); assert_eq!(Geodesic.point_at_ratio_between(start, start, ratio), start); assert_eq!(Rhumb.point_at_ratio_between(start, start, ratio), start); let ratio = 1.0; - assert_eq!(Haversine.point_at_ratio_between(start, start, ratio), start); + assert_eq!(HAVERSINE.point_at_ratio_between(start, start, ratio), start); assert_eq!(Euclidean.point_at_ratio_between(start, start, ratio), start); assert_eq!(Geodesic.point_at_ratio_between(start, start, ratio), start); assert_eq!(Rhumb.point_at_ratio_between(start, start, ratio), start); @@ -96,7 +96,7 @@ mod tests { let distance = 0.0; assert_eq!( - Haversine.point_at_distance_between(start, start, distance), + HAVERSINE.point_at_distance_between(start, start, distance), start ); @@ -116,7 +116,7 @@ mod tests { let due_north = Point::new(1.0, 1.9); let due_south = Point::new(1.0, 0.1); assert_relative_eq!( - Haversine.point_at_distance_between(start, start, distance), + HAVERSINE.point_at_distance_between(start, start, distance), due_north, epsilon = 1.0e-1 ); @@ -142,7 +142,7 @@ mod tests { let max_distance = 1.0; let include_ends = true; - let points: Vec<_> = Haversine + let points: Vec<_> = HAVERSINE .points_along_line(start, start, max_distance, include_ends) .collect(); assert_eq!(points, vec![start, start]); @@ -163,7 +163,7 @@ mod tests { assert_eq!(points, vec![start, start]); let include_ends = false; - let points: Vec<_> = Haversine + let points: Vec<_> = HAVERSINE .points_along_line(start, start, max_distance, include_ends) .collect(); assert_eq!(points, vec![]); diff --git a/geo/src/algorithm/line_measures/length.rs b/geo/src/algorithm/line_measures/length.rs index 5657485de..df1730fb4 100644 --- a/geo/src/algorithm/line_measures/length.rs +++ b/geo/src/algorithm/line_measures/length.rs @@ -87,7 +87,7 @@ impl LengthMeasurable for MultiLineString { #[cfg(test)] mod tests { use super::*; - use crate::{coord, Euclidean, Geodesic, Haversine, Rhumb}; + use crate::{coord, Euclidean, Geodesic, Rhumb, HAVERSINE}; #[test] fn lines() { @@ -107,7 +107,7 @@ mod tests { ); assert_eq!( 343_557., // meters - Haversine.length(&line).round() + HAVERSINE.length(&line).round() ); // computing Euclidean length of an unprojected (lng/lat) line gives a nonsense answer @@ -141,7 +141,7 @@ mod tests { ); assert_eq!( 6_304_387., // meters - Haversine.length(&line_string).round() + HAVERSINE.length(&line_string).round() ); // computing Euclidean length of an unprojected (lng/lat) gives a nonsense answer diff --git a/geo/src/algorithm/line_measures/metric_spaces/haversine.rs b/geo/src/algorithm/line_measures/metric_spaces/haversine.rs index bbf8c3bc7..3dd26db29 100644 --- a/geo/src/algorithm/line_measures/metric_spaces/haversine.rs +++ b/geo/src/algorithm/line_measures/metric_spaces/haversine.rs @@ -8,73 +8,62 @@ use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS}; /// /// Distances are considered [great circle] lengths and are measured in meters. /// -/// # References -/// -/// *Note*: this implementation uses a mean earth radius of 6371.0088 km (6_371_008.7714 m), based on the recommendation of -/// the IUGG: -/// -/// Moritz, H. (2000). Geodetic Reference System 1980. Journal of Geodesy, 74(1), 128–133. doi:10.1007/s001900050278 -/// "Derived Geometric Constants: **R1: mean radius**" (p131) -/// - -/// - -/// -/// If you'd like to use a different radius, see [`CustomHaversine`]. -/// -/// [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula// -/// [great circle]: https://en.wikipedia.org/wiki/Great_circle -pub struct Haversine; - -/// A spherical model using the [haversine formula] and a custom radius. -/// -/// Distances are considered [great circle] lengths and are measured in meters. -/// -/// ⚠️ For normal spherical measurements of the Earth, you probably want to use [`Haversine`] instead. -/// Use [`CustomHaversine`] only if you need to measure some non-Earth sphere, or if you want to -/// use a different value for Earth's radius. +/// # Examples /// /// ``` /// # use approx::assert_relative_eq; -/// use geo::{wkt, CustomHaversine, Haversine, Distance}; +/// use geo::{wkt, Haversine, Distance}; /// /// let start = wkt!(POINT(23.319941 42.698334)); // Sofia: Longitude, Latitude /// let finish = wkt!(POINT(24.742168 42.136097)); // Plovdiv: Longitude, Latitude /// -/// // Typically, you can just use `Haversine` for measuring on the Earth's surface +/// // Typically, you can just use `HAVERSINE` for measuring on the Earth's surface /// assert_relative_eq!( /// 132433.09929460194, -/// Haversine.distance(start, finish) +/// HAVERSINE.distance(start, finish) /// ); /// -/// // Beneath the hood, `Haversine` uses the mean radius of the GRS80 ellipsoid. +/// // `HAVERSINE` has a radius equal to the mean radius of the GRS80 ellipsoid. /// assert_relative_eq!( -/// Haversine.distance(start, finish), -/// CustomHaversine::GRS80_MEAN_RADIUS.distance(start, finish) +/// HAVERSINE.distance(start, finish), +/// Haversine::GRS80_MEAN_RADIUS.distance(start, finish) /// ); /// /// // You may choose to use one of the other well known estimations of the Earth's radius, /// // which may result in *slightly* different results. /// assert_relative_eq!( /// 132433.06564071847, -/// CustomHaversine::GRS80_EQUAL_AREA.distance(start, finish) +/// Haversine::GRS80_EQUAL_AREA.distance(start, finish) /// ); /// /// // Or you can specify whatever radius you want to get some "out of this world" results. -/// let mars_sphere = CustomHaversine::new(3_389_500.0); // 👽 Mars radius in meters +/// let mars_sphere = Haversine::new(3_389_500.0); // 👽 Mars radius in meters /// assert_relative_eq!( /// 70456.97222377927, /// mars_sphere.distance(start, finish) /// ); /// ``` /// +/// You may specify a custom radius for the Earth (or other sphere), but for normal spherical +/// measurements of the Earth, you probably want to just use [`HAVERSINE`] which uses the +/// earth radius of 6371.0088 km (6_371_008.7714 m), based on the recommendation of the IUGG. +/// /// # References /// +/// Moritz, H. (2000). Geodetic Reference System 1980. Journal of Geodesy, 74(1), 128–133. doi:10.1007/s001900050278 +/// "Derived Geometric Constants: **R1: mean radius**" (p131) +/// - +/// - +/// /// [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula// /// [great circle]: https://en.wikipedia.org/wiki/Great_circle -pub struct CustomHaversine { +pub struct Haversine { radius: f64, } -impl CustomHaversine { +pub const HAVERSINE: Haversine = Haversine::GRS80_MEAN_RADIUS; + +impl Haversine { /// ## Parameters /// - radius: The radius of the sphere, typically in meters. pub const fn new(radius: f64) -> Self { @@ -125,7 +114,7 @@ impl CustomHaversine { }; } -impl Default for CustomHaversine { +impl Default for Haversine { fn default() -> Self { Self::GRS80_MEAN_RADIUS } @@ -148,7 +137,7 @@ impl Bearing for Haversine { /// /// let origin = Point::new(9.0, 10.0); /// let destination = Point::new(9.5, 10.1); - /// let bearing = Haversine.bearing(origin, destination); + /// let bearing = HAVERSINE.bearing(origin, destination); /// // A little north of east /// assert_relative_eq!(bearing, 78.47, epsilon = 1.0e-2); /// ``` @@ -159,13 +148,6 @@ impl Bearing for Haversine { /// () /// /// [great circle]: https://en.wikipedia.org/wiki/Great_circle - fn bearing(&self, origin: Point, destination: Point) -> F { - CustomHaversine::default().bearing(origin, destination) - } -} - -impl Bearing for CustomHaversine { - /// Same as [`Haversine.bearing`], but on a sphere with a custom radius. fn bearing(&self, origin: Point, destination: Point) -> F { let three_sixty = F::from(360.0).expect("Numeric type to be constructable from primitive 360"); @@ -199,7 +181,7 @@ impl Destination for Haversine { /// use geo::Point; /// /// let origin = Point::new(9.177789688110352, 48.776781529534965); - /// let destination = Haversine.destination(origin, 45., 10000.); + /// let destination = HAVERSINE.destination(origin, 45., 10000.); /// assert_relative_eq!(Point::new(9.274409949623532, 48.84033274015048), destination); /// ``` /// @@ -209,13 +191,6 @@ impl Destination for Haversine { /// the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf) /// /// [great circle]: https://en.wikipedia.org/wiki/Great_circle - fn destination(&self, origin: Point, bearing: F, meters: F) -> Point { - CustomHaversine::default().destination(origin, bearing, meters) - } -} - -impl Destination for CustomHaversine { - /// Same as [`Haversine.destination`], but on a sphere with a custom radius. fn destination(&self, origin: Point, bearing: F, meters: F) -> Point { let center_lng = origin.x().to_radians(); let center_lat = origin.y().to_radians(); @@ -252,7 +227,7 @@ impl Distance, Point> for Haversin /// 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); + /// let distance = HAVERSINE.distance(new_york_city, london); /// /// assert_relative_eq!( /// 5_570_230., // meters @@ -266,13 +241,6 @@ impl Distance, Point> for Haversin /// the IUGG](ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf) /// /// [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula - fn distance(&self, origin: Point, destination: Point) -> F { - CustomHaversine::default().distance(origin, destination) - } -} - -impl Distance, Point> for CustomHaversine { - /// Same as [`Haversine.distance`], but on a sphere with a custom radius. fn distance(&self, origin: Point, destination: Point) -> F { let two = F::one() + F::one(); let theta1 = origin.y().to_radians(); @@ -302,10 +270,10 @@ impl InterpolatePoint for Haversine { /// 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); + /// 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); + /// 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); /// ``` /// @@ -316,7 +284,8 @@ impl InterpolatePoint for Haversine { end: Point, meters_from_start: F, ) -> Point { - CustomHaversine::default().point_at_distance_between(start, end, meters_from_start) + 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. @@ -331,13 +300,13 @@ impl InterpolatePoint for Haversine { /// 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); + /// 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); + /// 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); + /// 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); /// ``` /// @@ -348,7 +317,14 @@ impl InterpolatePoint for Haversine { end: Point, ratio_from_start: F, ) -> Point { - CustomHaversine::default().point_at_ratio_between(start, end, ratio_from_start) + if start == end || ratio_from_start == F::zero() { + return start; + } + if ratio_from_start == F::one() { + return end; + } + let calculation = HaversineIntermediateFillCalculation::new(start, end); + calculation.point_at_ratio(ratio_from_start) } /// Interpolates `Point`s along a [great circle] between `start` and `end`. @@ -361,50 +337,6 @@ impl InterpolatePoint for Haversine { /// /// [great circle]: https://en.wikipedia.org/wiki/Great_circle /// [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula - fn points_along_line( - &self, - start: Point, - end: Point, - max_distance: F, - include_ends: bool, - ) -> impl Iterator> { - CustomHaversine::GRS80_MEAN_RADIUS.points_along_line(start, end, max_distance, include_ends) - } -} - -/// Interpolate Point(s) along a [great circle]. -/// -/// [great circle]: https://en.wikipedia.org/wiki/Great_circle -impl InterpolatePoint for CustomHaversine { - /// Same as [`Haversine.point_at_distance_between`], but on a sphere with a custom radius. - fn point_at_distance_between( - &self, - start: Point, - end: Point, - meters_from_start: F, - ) -> Point { - let bearing = self.bearing(start, end); - self.destination(start, bearing, meters_from_start) - } - - /// Same as [`Haversine.point_at_ratio_between`], but on a sphere with a custom radius. - fn point_at_ratio_between( - &self, - start: Point, - end: Point, - ratio_from_start: F, - ) -> Point { - if start == end || ratio_from_start == F::zero() { - return start; - } - if ratio_from_start == F::one() { - return end; - } - let calculation = HaversineIntermediateFillCalculation::new(start, end); - calculation.point_at_ratio(ratio_from_start) - } - - /// Same as [`Haversine.points_along_line`], but on a sphere with a custom radius. fn points_along_line( &self, start: Point, @@ -534,28 +466,28 @@ mod tests { fn north() { let origin = Point::new(0.0, 0.0); let destination = Point::new(0.0, 1.0); - assert_relative_eq!(0.0, Haversine.bearing(origin, destination)); + assert_relative_eq!(0.0, HAVERSINE.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, Haversine.bearing(origin, destination)); + assert_relative_eq!(90.0, HAVERSINE.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, Haversine.bearing(origin, destination)); + assert_relative_eq!(180.0, HAVERSINE.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, Haversine.bearing(origin, destination)); + assert_relative_eq!(270.0, HAVERSINE.bearing(origin, destination)); } } @@ -568,7 +500,7 @@ mod tests { let bearing = 0.0; assert_relative_eq!( Point::new(0.0, 0.899320363724538), - Haversine.destination(origin, bearing, 100_000.0) + HAVERSINE.destination(origin, bearing, 100_000.0) ); } @@ -578,7 +510,7 @@ mod tests { let bearing = 90.0; assert_relative_eq!( Point::new(0.8993203637245415, 5.506522912913066e-17), - Haversine.destination(origin, bearing, 100_000.0) + HAVERSINE.destination(origin, bearing, 100_000.0) ); } @@ -588,7 +520,7 @@ mod tests { let bearing = 180.0; assert_relative_eq!( Point::new(0.0, -0.899320363724538), - Haversine.destination(origin, bearing, 100_000.0) + HAVERSINE.destination(origin, bearing, 100_000.0) ); } @@ -598,7 +530,7 @@ mod tests { let bearing = 270.0; assert_relative_eq!( Point::new(-0.8993203637245415, -1.6519568738739197e-16), - Haversine.destination(origin, bearing, 100_000.0) + HAVERSINE.destination(origin, bearing, 100_000.0) ); } } @@ -611,7 +543,7 @@ mod tests { 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); + let distance = HAVERSINE.distance(new_york_city, london); assert_relative_eq!( 5_570_230., // meters @@ -626,7 +558,7 @@ mod tests { fn point_at_ratio_between_midpoint() { let start = Point::new(10.0, 20.0); let end = Point::new(125.0, 25.0); - let midpoint = Haversine.point_at_ratio_between(start, end, 0.5); + let midpoint = HAVERSINE.point_at_ratio_between(start, end, 0.5); assert_relative_eq!(midpoint, Point::new(65.87394172511485, 37.61809316888599)); } #[test] @@ -634,7 +566,7 @@ mod tests { let start = Point::new(10.0, 20.0); let end = Point::new(125.0, 25.0); let max_dist = 1000000.0; // meters - let route = Haversine + let route = HAVERSINE .points_along_line(start, end, max_dist, true) .collect::>(); assert_eq!(route.len(), 13); @@ -647,7 +579,7 @@ mod tests { let start = Point::new(10.0, 20.0); let end = Point::new(125.0, 25.0); let max_dist = 1000000.0; // meters - let route = Haversine + let route = HAVERSINE .points_along_line(start, end, max_dist, false) .collect::>(); assert_eq!(route.len(), 11); diff --git a/geo/src/algorithm/line_measures/metric_spaces/mod.rs b/geo/src/algorithm/line_measures/metric_spaces/mod.rs index d9c582a33..c90a39375 100644 --- a/geo/src/algorithm/line_measures/metric_spaces/mod.rs +++ b/geo/src/algorithm/line_measures/metric_spaces/mod.rs @@ -5,7 +5,7 @@ mod geodesic; pub use geodesic::Geodesic; mod haversine; -pub use haversine::{CustomHaversine, Haversine}; +pub use haversine::{Haversine, HAVERSINE}; mod rhumb; pub use rhumb::Rhumb; diff --git a/geo/src/algorithm/line_measures/mod.rs b/geo/src/algorithm/line_measures/mod.rs index 2cee63c47..e3bdb66b1 100644 --- a/geo/src/algorithm/line_measures/mod.rs +++ b/geo/src/algorithm/line_measures/mod.rs @@ -19,4 +19,4 @@ mod densify; pub use densify::{Densifiable, Densify}; pub mod metric_spaces; -pub use metric_spaces::{Euclidean, Geodesic, Haversine, Rhumb}; +pub use metric_spaces::{Euclidean, Geodesic, Haversine, Rhumb, HAVERSINE}; diff --git a/geo/src/algorithm/linestring_segment.rs b/geo/src/algorithm/linestring_segment.rs index 4ab85bb38..11d33bfb4 100644 --- a/geo/src/algorithm/linestring_segment.rs +++ b/geo/src/algorithm/linestring_segment.rs @@ -1,6 +1,7 @@ use crate::algorithm::{Densify, Length, LineInterpolatePoint, LinesIter}; use crate::geometry::{Coord, LineString, MultiLineString}; -use crate::line_measures::{Euclidean, Haversine}; +use crate::line_measures::Euclidean; +use crate::HAVERSINE; /// Segments a LineString into `segment_count` equal length LineStrings as a MultiLineString /// using Euclidean distance calculations. See `LineStringSegmentizeHaversine` @@ -114,7 +115,7 @@ implement_segmentize!(LineStringSegmentize, line_segmentize, Euclidean); implement_segmentize!( LineStringSegmentizeHaversine, line_segmentize_haversine, - Haversine + HAVERSINE ); #[cfg(test)] @@ -326,7 +327,7 @@ mod test { let lens = segments .0 .iter() - .map(|li| Haversine.length(li)) + .map(|li| HAVERSINE.length(li)) .collect::>(); let epsilon = 1e-6; // 6th decimal place which is micrometers @@ -342,7 +343,7 @@ mod test { ] .into(); - assert_relative_eq!(Haversine.length(&linestring), 83.3523000093029); + assert_relative_eq!(HAVERSINE.length(&linestring), 83.3523000093029); let n = 8; @@ -350,8 +351,8 @@ mod test { // different at 12th decimal which is a picometer assert_relative_eq!( - Haversine.length(&linestring), - Haversine.length(&segments), + HAVERSINE.length(&linestring), + HAVERSINE.length(&segments), epsilon = 1e-11 ); } diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index 05da29174..57f977420 100644 --- a/geo/src/algorithm/mod.rs +++ b/geo/src/algorithm/mod.rs @@ -187,7 +187,7 @@ pub mod lines_iter; pub use lines_iter::LinesIter; pub mod line_measures; -pub use line_measures::metric_spaces::{CustomHaversine, Euclidean, Geodesic, Haversine, Rhumb}; +pub use line_measures::metric_spaces::{Euclidean, Geodesic, Haversine, Rhumb, HAVERSINE}; pub use line_measures::{Bearing, Densify, Destination, Distance, InterpolatePoint, Length}; /// Split a LineString into n segments