From 5f6bc012c364ef35ad2c55a62ef6ca04cd3e8fa8 Mon Sep 17 00:00:00 2001 From: Michael Kirk Date: Mon, 7 Oct 2024 15:22:10 -0700 Subject: [PATCH] Geodesic, Rhumb, Haversine line measures: deprecate legacy traits after moving implementation to modern traits Moved some implementation details from HaversineIntermediateFill to Haversine. --- geo/CHANGES.md | 5 + geo/benches/geodesic_distance.rs | 6 +- geo/src/algorithm/cross_track_distance.rs | 14 +- geo/src/algorithm/densify_haversine.rs | 7 +- geo/src/algorithm/geodesic_bearing.rs | 11 ++ geo/src/algorithm/geodesic_destination.rs | 17 +- geo/src/algorithm/geodesic_distance.rs | 11 +- geo/src/algorithm/geodesic_intermediate.rs | 67 +++----- geo/src/algorithm/geodesic_length.rs | 5 +- geo/src/algorithm/haversine_bearing.rs | 16 +- geo/src/algorithm/haversine_closest_point.rs | 24 +-- geo/src/algorithm/haversine_destination.rs | 33 ++-- geo/src/algorithm/haversine_distance.rs | 47 +++-- geo/src/algorithm/haversine_intermediate.rs | 155 +++++------------ geo/src/algorithm/haversine_length.rs | 4 +- .../line_measures/metric_spaces/geodesic.rs | 67 ++++++-- .../line_measures/metric_spaces/haversine.rs | 160 ++++++++++++++++-- .../line_measures/metric_spaces/rhumb.rs | 33 ++-- geo/src/algorithm/mod.rs | 11 +- geo/src/algorithm/rhumb/bearing.rs | 25 ++- geo/src/algorithm/rhumb/destination.rs | 25 ++- geo/src/algorithm/rhumb/distance.rs | 42 ++--- geo/src/algorithm/rhumb/intermediate.rs | 45 +++-- geo/src/algorithm/rhumb/length.rs | 5 +- geo/src/algorithm/rhumb/mod.rs | 18 +- 25 files changed, 514 insertions(+), 339 deletions(-) diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 88ef46d917..33f725dbd1 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -38,6 +38,11 @@ Haversine::distance(p1, p2) ``` * +* Deprecated legacy line measure traits in favor of those added in the previous changelog entry: + * `GeodesicBearing`, `GeodesicDistance`, `GeodesicDestination`, `GeodesicIntermediate` + * `RhumbBearing`, `RhumbDistance`, `RhumbDestination`, `RhumbIntermediate` + * `HaversineBearing`, `HaversineDistance`, `HaversineDestination`, `HaversineIntermediate` + * * Change IntersectionMatrix::is_equal_topo to now consider empty geometries as equal. * * Fix `(LINESTRING EMPTY).contains(LINESTRING EMPTY)` and `(MULTIPOLYGON EMPTY).contains(MULTIPOINT EMPTY)` which previously diff --git a/geo/benches/geodesic_distance.rs b/geo/benches/geodesic_distance.rs index e4cfb5f227..cee6d4f193 100644 --- a/geo/benches/geodesic_distance.rs +++ b/geo/benches/geodesic_distance.rs @@ -1,5 +1,5 @@ use criterion::{criterion_group, criterion_main}; -use geo::algorithm::GeodesicDistance; +use geo::{Distance, Geodesic}; fn criterion_benchmark(c: &mut criterion::Criterion) { c.bench_function("geodesic distance f64", |bencher| { @@ -7,9 +7,7 @@ fn criterion_benchmark(c: &mut criterion::Criterion) { let b = geo::Point::new(16.372477, 48.208810); bencher.iter(|| { - criterion::black_box( - criterion::black_box(&a).geodesic_distance(criterion::black_box(&b)), - ); + criterion::black_box(criterion::black_box(Geodesic::distance(a, b))); }); }); } diff --git a/geo/src/algorithm/cross_track_distance.rs b/geo/src/algorithm/cross_track_distance.rs index 37cfb8bc1a..7088d4b6b8 100644 --- a/geo/src/algorithm/cross_track_distance.rs +++ b/geo/src/algorithm/cross_track_distance.rs @@ -1,4 +1,4 @@ -use crate::{HaversineBearing, HaversineDistance, 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 = line_point_a.haversine_distance(self) / mean_earth_radius; - let theta_13: T = line_point_a.haversine_bearing(*self).to_radians(); - let theta_12: T = line_point_a.haversine_bearing(*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() } @@ -54,8 +54,8 @@ where #[cfg(test)] mod test { use crate::CrossTrackDistance; - use crate::HaversineDistance; use crate::Point; + 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), - p.haversine_distance(&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), - p.haversine_distance(&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 bd7cb94360..4f2465443f 100644 --- a/geo/src/algorithm/densify_haversine.rs +++ b/geo/src/algorithm/densify_haversine.rs @@ -1,12 +1,12 @@ use num_traits::FromPrimitive; +use crate::line_measures::{Haversine, InterpolatePoint}; +use crate::HaversineLength; use crate::{ CoordFloat, CoordsIter, Line, LineString, MultiLineString, MultiPolygon, Point, Polygon, Rect, Triangle, }; -use crate::{HaversineIntermediate, HaversineLength}; - /// Returns a new spherical geometry containing both existing and new interpolated coordinates with /// a maximum distance of `max_distance` between them. /// @@ -52,8 +52,7 @@ fn densify_line( let ratio = frac * T::from(segment_idx).unwrap(); let start = line.start; let end = line.end; - let interpolated_point = - Point::from(start).haversine_intermediate(&Point::from(end), ratio); + let interpolated_point = Haversine::point_at_ratio_between(Point(start), Point(end), ratio); container.push(interpolated_point); } } diff --git a/geo/src/algorithm/geodesic_bearing.rs b/geo/src/algorithm/geodesic_bearing.rs index d3ad16ec46..a1666506b3 100644 --- a/geo/src/algorithm/geodesic_bearing.rs +++ b/geo/src/algorithm/geodesic_bearing.rs @@ -8,17 +8,23 @@ use geographiclib_rs::{Geodesic, InverseGeodesic}; /// /// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf pub trait GeodesicBearing { + #[deprecated( + since = "0.29.0", + note = "Please use the `Geodesic::bearing` method from the `Bearing` trait instead" + )] /// Returns the bearing to another Point in degrees, where North is 0° and East is 90°. /// /// # Examples /// /// ``` /// # use approx::assert_relative_eq; + /// # #[allow(deprecated)] /// use geo::GeodesicBearing; /// use geo::Point; /// /// let p_1 = Point::new(9.177789688110352, 48.776781529534965); /// let p_2 = Point::new(9.27411867078536, 48.8403266058781); + /// # #[allow(deprecated)] /// let bearing = p_1.geodesic_bearing(p_2); /// assert_relative_eq!(bearing, 45., epsilon = 1.0e-6); /// ``` @@ -69,6 +75,7 @@ mod test { fn north_bearing() { let p_1 = point!(x: 9., y: 47.); let p_2 = point!(x: 9., y: 48.); + #[allow(deprecated)] let bearing = p_1.geodesic_bearing(p_2); assert_relative_eq!(bearing, 0.); } @@ -77,6 +84,7 @@ mod test { fn east_bearing() { let p_1 = point!(x: 9., y: 10.); let p_2 = point!(x: 18.118501133357412, y: 9.875322179340463); + #[allow(deprecated)] let bearing = p_1.geodesic_bearing(p_2); assert_relative_eq!(bearing, 90.); } @@ -85,14 +93,17 @@ mod test { fn northeast_bearing() { let p_1 = point!(x: 9.177789688110352f64, y: 48.776781529534965); let p_2 = point!(x: 9.27411867078536, y: 48.8403266058781); + #[allow(deprecated)] let bearing = p_1.geodesic_bearing(p_2); assert_relative_eq!(bearing, 45., epsilon = 1.0e-11); } #[test] fn consistent_with_destination() { + #[allow(deprecated)] use crate::algorithm::GeodesicDestination; let p_1 = point!(x: 9.177789688110352, y: 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.geodesic_destination(45., 10000.); let (bearing, distance) = p_1.geodesic_bearing_distance(p_2); assert_relative_eq!(bearing, 45., epsilon = 1.0e-11); diff --git a/geo/src/algorithm/geodesic_destination.rs b/geo/src/algorithm/geodesic_destination.rs index eb0f614d6a..25261a9bcb 100644 --- a/geo/src/algorithm/geodesic_destination.rs +++ b/geo/src/algorithm/geodesic_destination.rs @@ -1,7 +1,11 @@ +use crate::line_measures::{Destination, Geodesic}; use crate::Point; use geo_types::CoordNum; -use geographiclib_rs::{DirectGeodesic, Geodesic}; +#[deprecated( + since = "0.29.0", + note = "Please use the `Geodesic::destination` method from the `Destination` trait instead" +)] /// Returns a new Point using the distance to the existing Point and a bearing for the direction on a geodesic. /// /// This uses the geodesic methods given by [Karney (2013)]. @@ -18,6 +22,7 @@ pub trait GeodesicDestination { /// # Examples /// /// ```rust + /// # #[allow(deprecated)] /// use geo::GeodesicDestination; /// use geo::Point; /// @@ -26,6 +31,7 @@ pub trait GeodesicDestination { /// let northeast_bearing = 45.0; /// let distance = 10e6; /// + /// # #[allow(deprecated)] /// let p_1 = jfk.geodesic_destination(northeast_bearing, distance); /// use approx::assert_relative_eq; /// assert_relative_eq!(p_1.x(), 49.052487092959836); @@ -34,21 +40,23 @@ pub trait GeodesicDestination { fn geodesic_destination(&self, bearing: T, distance: T) -> Point; } +#[allow(deprecated)] impl GeodesicDestination for Point { fn geodesic_destination(&self, bearing: f64, distance: f64) -> Point { - let (lat, lon) = Geodesic::wgs84().direct(self.y(), self.x(), bearing, distance); - Point::new(lon, lat) + Geodesic::destination(*self, bearing, distance) } } #[cfg(test)] mod test { use super::*; + #[allow(deprecated)] use crate::GeodesicDistance; #[test] fn returns_a_new_point() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.geodesic_destination(45., 10000.); assert_relative_eq!( @@ -57,6 +65,7 @@ mod test { epsilon = 1.0e-6 ); + #[allow(deprecated)] let distance = p_1.geodesic_distance(&p_2); assert_relative_eq!(distance, 10000., epsilon = 1.0e-6) } @@ -64,6 +73,7 @@ mod test { #[test] fn bearing_zero_is_north() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.geodesic_destination(0., 1000.); assert_relative_eq!(p_1.x(), p_2.x(), epsilon = 1.0e-6); assert!(p_2.y() > p_1.y()) @@ -72,6 +82,7 @@ mod test { #[test] fn bearing_90_is_east() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.geodesic_destination(90., 1000.); assert_relative_eq!(p_1.y(), p_2.y(), epsilon = 1.0e-6); assert!(p_2.x() > p_1.x()) diff --git a/geo/src/algorithm/geodesic_distance.rs b/geo/src/algorithm/geodesic_distance.rs index c1ea101e2f..87d1755cc5 100644 --- a/geo/src/algorithm/geodesic_distance.rs +++ b/geo/src/algorithm/geodesic_distance.rs @@ -1,6 +1,9 @@ -use crate::Point; -use geographiclib_rs::{Geodesic, InverseGeodesic}; +use crate::{Distance, Geodesic, Point}; +#[deprecated( + since = "0.29.0", + note = "Please use the `Geodesic::distance` method from the `Distance` trait instead" +)] /// Determine the distance between two geometries on an ellipsoidal model of the earth. /// /// This uses the geodesic measurement methods given by [Karney (2013)]. As opposed to older methods @@ -28,6 +31,7 @@ pub trait GeodesicDistance { /// // London /// let p2 = point!(x: -0.1278, y: 51.5074); /// + /// # #[allow(deprecated)] /// let distance = p1.geodesic_distance(&p2); /// /// assert_eq!( @@ -39,8 +43,9 @@ pub trait GeodesicDistance { fn geodesic_distance(&self, rhs: &Rhs) -> T; } +#[allow(deprecated)] impl GeodesicDistance for Point { fn geodesic_distance(&self, rhs: &Point) -> f64 { - Geodesic::wgs84().inverse(self.y(), self.x(), rhs.y(), rhs.x()) + Geodesic::distance(*self, *rhs) } } diff --git a/geo/src/algorithm/geodesic_intermediate.rs b/geo/src/algorithm/geodesic_intermediate.rs index 3670fbba5e..109cdb9e2a 100644 --- a/geo/src/algorithm/geodesic_intermediate.rs +++ b/geo/src/algorithm/geodesic_intermediate.rs @@ -1,21 +1,32 @@ -use crate::{CoordFloat, Point}; -use geographiclib_rs::{DirectGeodesic, Geodesic, InverseGeodesic}; +use crate::{CoordFloat, Geodesic, InterpolatePoint, Point}; +#[deprecated( + since = "0.29.0", + note = "Please use the `InterpolatePoint` trait instead" +)] /// Returns a new Point along a route between two existing points on an ellipsoidal model of the earth pub trait GeodesicIntermediate { + #[deprecated( + since = "0.29.0", + note = "Please use `Geodesic::point_at_ratio_between` from the `InterpolatePoint` trait instead" + )] /// Returns a new Point along a route between two existing points on an ellipsoidal model of the earth /// /// # Examples /// /// ``` /// # use approx::assert_relative_eq; + /// # #[allow(deprecated)] /// use geo::GeodesicIntermediate; /// use geo::Point; /// /// let p1 = Point::new(10.0, 20.0); /// let p2 = Point::new(125.0, 25.0); + /// # #[allow(deprecated)] /// let i20 = p1.geodesic_intermediate(&p2, 0.2); + /// # #[allow(deprecated)] /// let i50 = p1.geodesic_intermediate(&p2, 0.5); + /// # #[allow(deprecated)] /// let i80 = p1.geodesic_intermediate(&p2, 0.8); /// let i20_should = Point::new(29.842907, 29.951445); /// let i50_should = Point::new(65.879360, 37.722253); @@ -25,6 +36,11 @@ pub trait GeodesicIntermediate { /// assert_relative_eq!(i80, i80_should, epsilon = 1.0e-6); /// ``` fn geodesic_intermediate(&self, other: &Point, f: T) -> Point; + + #[deprecated( + since = "0.29.0", + note = "Please use `Geodesic::points_along_line` from the `InterpolatePoint` trait instead" + )] fn geodesic_intermediate_fill( &self, other: &Point, @@ -33,15 +49,10 @@ pub trait GeodesicIntermediate { ) -> Vec>; } +#[allow(deprecated)] impl GeodesicIntermediate for Point { fn geodesic_intermediate(&self, other: &Point, f: f64) -> Point { - let g = Geodesic::wgs84(); - let (total_distance, azi1, _azi2, _a12) = - g.inverse(self.y(), self.x(), other.y(), other.x()); - let distance = total_distance * f; - let (lat2, lon2) = g.direct(self.y(), self.x(), azi1, distance); - - Point::new(lon2, lat2) + Geodesic::point_at_ratio_between(*self, *other, f) } fn geodesic_intermediate_fill( @@ -50,36 +61,7 @@ impl GeodesicIntermediate for Point { max_dist: f64, include_ends: bool, ) -> Vec { - let g = Geodesic::wgs84(); - let (total_distance, azi1, _azi2, _a12) = - g.inverse(self.y(), self.x(), other.y(), other.x()); - - if total_distance <= max_dist { - return if include_ends { - vec![*self, *other] - } else { - vec![] - }; - } - - let number_of_points = (total_distance / max_dist).ceil(); - let interval = 1.0 / number_of_points; - - let mut current_step = interval; - let mut points = if include_ends { vec![*self] } else { vec![] }; - - while current_step < 1.0 { - let (lat2, lon2) = g.direct(self.y(), self.x(), azi1, total_distance * current_step); - let point = Point::new(lon2, lat2); - points.push(point); - current_step += interval; - } - - if include_ends { - points.push(*other); - } - - points + Geodesic::points_along_line(*self, *other, max_dist, include_ends).collect() } } @@ -92,7 +74,9 @@ mod tests { fn f_is_zero_or_one_test() { let p1 = Point::new(10.0, 20.0); let p2 = Point::new(15.0, 25.0); + #[allow(deprecated)] let i0 = p1.geodesic_intermediate(&p2, 0.0); + #[allow(deprecated)] let i100 = p1.geodesic_intermediate(&p2, 1.0); assert_relative_eq!(i0, p1, epsilon = 1.0e-6); assert_relative_eq!(i100, p2, epsilon = 1.0e-6); @@ -102,8 +86,11 @@ mod tests { fn various_f_values_test() { let p1 = Point::new(10.0, 20.0); let p2 = Point::new(125.0, 25.0); + #[allow(deprecated)] let i20 = p1.geodesic_intermediate(&p2, 0.2); + #[allow(deprecated)] let i50 = p1.geodesic_intermediate(&p2, 0.5); + #[allow(deprecated)] let i80 = p1.geodesic_intermediate(&p2, 0.8); let i20_should = Point::new(29.842907, 29.951445); let i50_should = Point::new(65.879360, 37.722253); @@ -119,7 +106,9 @@ mod tests { let p2 = Point::new(40.0, 50.0); let max_dist = 1000000.0; // meters let include_ends = true; + #[allow(deprecated)] let i50 = p1.geodesic_intermediate(&p2, 0.5); + #[allow(deprecated)] let route = p1.geodesic_intermediate_fill(&p2, max_dist, include_ends); assert_eq!(route, vec![p1, i50, p2]); } diff --git a/geo/src/algorithm/geodesic_length.rs b/geo/src/algorithm/geodesic_length.rs index b9a7e22aef..41f09a1b19 100644 --- a/geo/src/algorithm/geodesic_length.rs +++ b/geo/src/algorithm/geodesic_length.rs @@ -1,5 +1,4 @@ -use crate::GeodesicDistance; -use crate::{Line, LineString, MultiLineString}; +use crate::{Distance, Geodesic, Line, LineString, MultiLineString}; /// Determine the length of a geometry on an ellipsoidal model of the earth. /// @@ -49,7 +48,7 @@ impl GeodesicLength for Line { /// The units of the returned value is meters. fn geodesic_length(&self) -> f64 { let (start, end) = self.points(); - start.geodesic_distance(&end) + Geodesic::distance(start, end) } } diff --git a/geo/src/algorithm/haversine_bearing.rs b/geo/src/algorithm/haversine_bearing.rs index 37b7cd6113..031665270c 100644 --- a/geo/src/algorithm/haversine_bearing.rs +++ b/geo/src/algorithm/haversine_bearing.rs @@ -1,10 +1,13 @@ use crate::{CoordFloat, Point}; +#[deprecated( + since = "0.29.0", + note = "Please use the `Haversine::bearing` method from the `Bearing` trait instead" +)] /// Returns the bearing to another Point in degrees. /// /// Bullock, R.: Great Circle Distances and Bearings Between Two Locations, 2007. /// () - pub trait HaversineBearing { /// Returns the bearing to another Point in degrees, where North is 0° and East is 90°. /// @@ -12,17 +15,20 @@ pub trait HaversineBearing { /// /// ``` /// # use approx::assert_relative_eq; + /// # #[allow(deprecated)] /// use geo::HaversineBearing; /// use geo::Point; /// /// let p_1 = Point::new(9.177789688110352, 48.776781529534965); /// let p_2 = Point::new(9.274410083250379, 48.84033282787534); + /// # #[allow(deprecated)] /// let bearing = p_1.haversine_bearing(p_2); /// assert_relative_eq!(bearing, 45., epsilon = 1.0e-6); /// ``` fn haversine_bearing(&self, point: Point) -> T; } +#[allow(deprecated)] impl HaversineBearing for Point where T: CoordFloat, @@ -41,13 +47,16 @@ where #[cfg(test)] mod test { use crate::point; + #[allow(deprecated)] use crate::HaversineBearing; + #[allow(deprecated)] use crate::HaversineDestination; #[test] fn north_bearing() { let p_1 = point!(x: 9., y: 47.); let p_2 = point!(x: 9., y: 48.); + #[allow(deprecated)] let bearing = p_1.haversine_bearing(p_2); assert_relative_eq!(bearing, 0.); } @@ -56,6 +65,7 @@ mod test { fn equatorial_east_bearing() { let p_1 = point!(x: 9., y: 0.); let p_2 = point!(x: 10., y: 0.); + #[allow(deprecated)] let bearing = p_1.haversine_bearing(p_2); assert_relative_eq!(bearing, 90.); } @@ -65,6 +75,7 @@ mod test { let p_1 = point!(x: 9., y: 10.); let p_2 = point!(x: 18.12961917258341, y: 9.875828894123304); + #[allow(deprecated)] let bearing = p_1.haversine_bearing(p_2); assert_relative_eq!(bearing, 90.); } @@ -73,6 +84,7 @@ mod test { fn northeast_bearing() { let p_1 = point!(x: 9.177789688110352f64, y: 48.776781529534965); let p_2 = point!(x: 9.274409949623548, y: 48.84033274015048); + #[allow(deprecated)] let bearing = p_1.haversine_bearing(p_2); assert_relative_eq!(bearing, 45., epsilon = 1.0e-6); } @@ -80,8 +92,10 @@ mod test { #[test] fn consistent_with_destination() { let p_1 = point!(x: 9.177789688110352f64, y: 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.haversine_destination(45., 10000.); + #[allow(deprecated)] let b_1 = p_1.haversine_bearing(p_2); assert_relative_eq!(b_1, 45., epsilon = 1.0e-6); } diff --git a/geo/src/algorithm/haversine_closest_point.rs b/geo/src/algorithm/haversine_closest_point.rs index ac896d4b1e..b722e2c12a 100644 --- a/geo/src/algorithm/haversine_closest_point.rs +++ b/geo/src/algorithm/haversine_closest_point.rs @@ -1,6 +1,6 @@ -use crate::{haversine_distance::HaversineDistance, HaversineBearing}; +use crate::line_measures::{Bearing, Destination, Distance, Haversine}; use crate::{Closest, Contains}; -use crate::{CoordsIter, GeoFloat, HaversineDestination, Point, MEAN_EARTH_RADIUS}; +use crate::{CoordsIter, GeoFloat, Point, MEAN_EARTH_RADIUS}; use geo_types::{ Coord, Geometry, GeometryCollection, Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Polygon, Rect, Triangle, @@ -92,7 +92,7 @@ where } // This can probably be done cheaper - let d3 = p2.haversine_distance(&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 = p1.haversine_bearing(*from).to_radians(); - let crs_ab = p1.haversine_bearing(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 = p2.haversine_bearing(*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 = p1.haversine_distance(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(p1.haversine_destination(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 = p2.haversine_distance(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 = pt.haversine_distance(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 = pt.haversine_distance(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 = pt.haversine_distance(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 0f000f3dc3..b3f9a088ca 100644 --- a/geo/src/algorithm/haversine_destination.rs +++ b/geo/src/algorithm/haversine_destination.rs @@ -1,6 +1,10 @@ -use crate::{utils::normalize_longitude, CoordFloat, Point, MEAN_EARTH_RADIUS}; +use crate::{CoordFloat, Destination, Haversine, Point}; use num_traits::FromPrimitive; +#[deprecated( + since = "0.29.0", + 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 /// /// *Note*: this implementation uses a mean earth radius of 6371.088 km, based on the [recommendation of @@ -16,54 +20,47 @@ pub trait HaversineDestination { /// # Examples /// /// ```rust + /// # #[allow(deprecated)] /// use geo::HaversineDestination; /// use geo::Point; /// use approx::assert_relative_eq; /// /// let p_1 = Point::new(9.177789688110352, 48.776781529534965); + /// # #[allow(deprecated)] /// let p_2 = p_1.haversine_destination(45., 10000.); /// assert_relative_eq!(p_2, Point::new(9.274409949623548, 48.84033274015048), epsilon = 1e-6) /// ``` fn haversine_destination(&self, bearing: T, distance: T) -> Point; } +#[allow(deprecated)] impl HaversineDestination for Point where T: CoordFloat + FromPrimitive, { fn haversine_destination(&self, bearing: T, distance: T) -> Point { - let center_lng = self.x().to_radians(); - let center_lat = self.y().to_radians(); - let bearing_rad = bearing.to_radians(); - - let rad = distance / T::from(MEAN_EARTH_RADIUS).unwrap(); - - let lat = - { center_lat.sin() * rad.cos() + center_lat.cos() * rad.sin() * bearing_rad.cos() } - .asin(); - let lng = { bearing_rad.sin() * rad.sin() * center_lat.cos() } - .atan2(rad.cos() - center_lat.sin() * lat.sin()) - + center_lng; - - Point::new(normalize_longitude(lng.to_degrees()), lat.to_degrees()) + Haversine::destination(*self, bearing, distance) } } #[cfg(test)] mod test { use super::*; + #[allow(deprecated)] use crate::{HaversineBearing, HaversineDistance}; use num_traits::pow; #[test] fn returns_a_new_point() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.haversine_destination(45., 10000.); assert_relative_eq!( p_2, Point::new(9.274409949623548, 48.84033274015048), epsilon = 1.0e-6 ); + #[allow(deprecated)] let distance = p_1.haversine_distance(&p_2); assert_relative_eq!(distance, 10000., epsilon = 1.0e-6) } @@ -71,9 +68,12 @@ mod test { #[test] fn direct_and_indirect_destinations_are_close() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.haversine_destination(45., 10000.); let square_edge = { pow(10000., 2) / 2f64 }.sqrt(); + #[allow(deprecated)] let p_3 = p_1.haversine_destination(0., square_edge); + #[allow(deprecated)] let p_4 = p_3.haversine_destination(90., square_edge); assert_relative_eq!(p_4, p_2, epsilon = 1.0e-6); } @@ -81,6 +81,7 @@ mod test { #[test] fn bearing_zero_is_north() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.haversine_destination(0., 1000.); assert_relative_eq!(p_1.x(), p_2.x(), epsilon = 1.0e-6); assert!(p_2.y() > p_1.y()) @@ -92,7 +93,9 @@ mod test { let pt2 = Point::new(-170.0, -30.0); for (start, end) in [(pt1, pt2), (pt2, pt1)] { + #[allow(deprecated)] let bearing = start.haversine_bearing(end); + #[allow(deprecated)] let results: Vec<_> = (0..8) .map(|n| start.haversine_destination(bearing, n as f64 * 250_000.)) .collect(); diff --git a/geo/src/algorithm/haversine_distance.rs b/geo/src/algorithm/haversine_distance.rs index f912a56198..b92031de2b 100644 --- a/geo/src/algorithm/haversine_distance.rs +++ b/geo/src/algorithm/haversine_distance.rs @@ -1,6 +1,10 @@ -use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS}; +use crate::{CoordFloat, Distance, Haversine, Point}; use num_traits::FromPrimitive; +#[deprecated( + since = "0.29.0", + note = "Please use the `Haversine::distance` method from the `Distance` trait instead" +)] /// Determine the distance between two geometries using the [haversine formula]. /// /// [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula @@ -27,6 +31,7 @@ pub trait HaversineDistance { /// // London /// let p2 = point!(x: -0.1278f64, y: 51.5074f64); /// + /// # #[allow(deprecated)] /// let distance = p1.haversine_distance(&p2); /// /// assert_eq!( @@ -39,25 +44,19 @@ pub trait HaversineDistance { fn haversine_distance(&self, rhs: &Rhs) -> T; } +#[allow(deprecated)] impl HaversineDistance> for Point where T: CoordFloat + FromPrimitive, { fn haversine_distance(&self, rhs: &Point) -> T { - let two = T::one() + T::one(); - let theta1 = self.y().to_radians(); - let theta2 = rhs.y().to_radians(); - let delta_theta = (rhs.y() - self.y()).to_radians(); - let delta_lambda = (rhs.x() - self.x()).to_radians(); - let a = (delta_theta / two).sin().powi(2) - + theta1.cos() * theta2.cos() * (delta_lambda / two).sin().powi(2); - let c = two * a.sqrt().asin(); - T::from(MEAN_EARTH_RADIUS).unwrap() * c + Haversine::distance(*self, *rhs) } } #[cfg(test)] mod test { + #[allow(deprecated)] use crate::HaversineDistance; use crate::Point; @@ -65,22 +64,18 @@ mod test { fn distance1_test() { let a = Point::new(0., 0.); let b = Point::new(1., 0.); - assert_relative_eq!( - a.haversine_distance(&b), - 111195.0802335329_f64, - epsilon = 1.0e-6 - ); + #[allow(deprecated)] + let distance = a.haversine_distance(&b); + assert_relative_eq!(distance, 111195.0802335329_f64, epsilon = 1.0e-6); } #[test] fn distance2_test() { let a = Point::new(-72.1235, 42.3521); let b = Point::new(72.1260, 70.612); - assert_relative_eq!( - a.haversine_distance(&b), - 7130580.307935911_f64, - epsilon = 1.0e-6 - ); + #[allow(deprecated)] + let distance = a.haversine_distance(&b); + assert_relative_eq!(distance, 7130580.307935911_f64, epsilon = 1.0e-6); } #[test] @@ -88,11 +83,9 @@ mod test { // this input comes from issue #100 let a = Point::new(-77.036585, 38.897448); let b = Point::new(-77.009080, 38.889825); - assert_relative_eq!( - a.haversine_distance(&b), - 2526.823504306046_f64, - epsilon = 1.0e-6 - ); + #[allow(deprecated)] + let distance = a.haversine_distance(&b); + assert_relative_eq!(distance, 2526.823504306046_f64, epsilon = 1.0e-6); } #[test] @@ -100,6 +93,8 @@ mod test { // this input comes from issue #100 let a = Point::::new(-77.03658, 38.89745); let b = Point::::new(-77.00908, 38.889825); - assert_relative_eq!(a.haversine_distance(&b), 2526.8354_f32, epsilon = 1.0e-6); + #[allow(deprecated)] + let distance = a.haversine_distance(&b); + assert_relative_eq!(distance, 2526.8354_f32, epsilon = 1.0e-6); } } diff --git a/geo/src/algorithm/haversine_intermediate.rs b/geo/src/algorithm/haversine_intermediate.rs index a735495b4a..55cf861a71 100644 --- a/geo/src/algorithm/haversine_intermediate.rs +++ b/geo/src/algorithm/haversine_intermediate.rs @@ -1,9 +1,16 @@ -use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS}; +use crate::{CoordFloat, Haversine, InterpolatePoint, Point}; use num_traits::FromPrimitive; +#[deprecated( + since = "0.29.0", + note = "Please use the `InterpolatePoint` trait instead" +)] /// Returns a new Point along a great circle route between two existing points - pub trait HaversineIntermediate { + #[deprecated( + since = "0.29.0", + 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`. /// /// * `other` - The other point to interpolate towards. @@ -14,20 +21,27 @@ pub trait HaversineIntermediate { /// /// ``` /// # use approx::assert_relative_eq; + /// # #[allow(deprecated)] /// use geo::HaversineIntermediate; /// use geo::Point; /// /// let p1 = Point::new(10.0, 20.0); /// let p2 = Point::new(125.0, 25.0); /// + /// # #[allow(deprecated)] /// let i20 = p1.haversine_intermediate(&p2, 0.2); /// assert_relative_eq!(i20, Point::new(29.8, 29.9), epsilon = 0.2); /// + /// # #[allow(deprecated)] /// let i80 = p1.haversine_intermediate(&p2, 0.8); /// assert_relative_eq!(i80, Point::new(103.5, 33.5), epsilon = 0.2); /// ``` fn haversine_intermediate(&self, other: &Point, ratio: T) -> Point; + #[deprecated( + since = "0.29.0", + note = "Please use `Haversine::points_along_line` from the `InterpolatePoint` trait instead" + )] /// Interpolates `Point`s along a great circle route between self and `other`. /// /// As many points as necessary will be added such that the distance between points @@ -42,13 +56,13 @@ pub trait HaversineIntermediate { ) -> Vec>; } +#[allow(deprecated)] impl HaversineIntermediate for Point where T: CoordFloat + FromPrimitive, { fn haversine_intermediate(&self, other: &Point, ratio: T) -> Point { - let params = get_params(self, other); - get_point(¶ms, ratio) + Haversine::point_at_ratio_between(*self, *other, ratio) } fn haversine_intermediate_fill( @@ -57,124 +71,23 @@ where max_dist: T, include_ends: bool, ) -> Vec> { - let params = get_params(self, other); - let HaversineParams { d, .. } = params; - - let total_distance = d * T::from(MEAN_EARTH_RADIUS).unwrap(); - - if total_distance <= max_dist { - return if include_ends { - vec![*self, *other] - } else { - vec![] - }; - } - - let number_of_points = (total_distance / max_dist).ceil(); - let interval = T::one() / number_of_points; - - let mut current_step = interval; - let mut points = if include_ends { vec![*self] } else { vec![] }; - - while current_step < T::one() { - let point = get_point(¶ms, current_step); - points.push(point); - current_step = current_step + interval; - } - - if include_ends { - points.push(*other); - } - - points - } -} - -#[allow(clippy::many_single_char_names)] -struct HaversineParams { - d: T, - n: T, - o: T, - p: T, - q: T, - r: T, - s: T, -} - -#[allow(clippy::many_single_char_names)] -fn get_point(params: &HaversineParams, ratio: T) -> Point { - let one = T::one(); - - let HaversineParams { - d, - n, - o, - p, - q, - r, - s, - } = *params; - - let a = ((one - ratio) * d).sin() / d.sin(); - let b = (ratio * d).sin() / d.sin(); - - let x = a * n + b * o; - let y = a * p + b * q; - let z = a * r + b * s; - - let lat = z.atan2(x.hypot(y)); - let lon = y.atan2(x); - - Point::new(lon.to_degrees(), lat.to_degrees()) -} - -#[allow(clippy::many_single_char_names)] -fn get_params(p1: &Point, p2: &Point) -> HaversineParams { - let one = T::one(); - let two = one + one; - - let lat1 = p1.y().to_radians(); - let lon1 = p1.x().to_radians(); - let lat2 = p2.y().to_radians(); - let lon2 = p2.x().to_radians(); - - let (lat1_sin, lat1_cos) = lat1.sin_cos(); - let (lat2_sin, lat2_cos) = lat2.sin_cos(); - let (lon1_sin, lon1_cos) = lon1.sin_cos(); - let (lon2_sin, lon2_cos) = lon2.sin_cos(); - - let m = lat1_cos * lat2_cos; - - let n = lat1_cos * lon1_cos; - let o = lat2_cos * lon2_cos; - let p = lat1_cos * lon1_sin; - let q = lat2_cos * lon2_sin; - - let k = (((lat1 - lat2) / two).sin().powi(2) + m * ((lon1 - lon2) / two).sin().powi(2)).sqrt(); - - let d = two * k.asin(); - - HaversineParams { - d, - n, - o, - p, - q, - r: lat1_sin, - s: lat2_sin, + Haversine::points_along_line(*self, *other, max_dist, include_ends).collect() } } #[cfg(test)] mod test { use super::*; + #[allow(deprecated)] use crate::HaversineIntermediate; #[test] fn f_is_zero_or_one_test() { let p1 = Point::new(10.0, 20.0); let p2 = Point::new(15.0, 25.0); + #[allow(deprecated)] let i0 = p1.haversine_intermediate(&p2, 0.0); + #[allow(deprecated)] let i100 = p1.haversine_intermediate(&p2, 1.0); assert_relative_eq!(i0.x(), p1.x(), epsilon = 1.0e-6); assert_relative_eq!(i0.y(), p1.y(), epsilon = 1.0e-6); @@ -186,8 +99,11 @@ mod test { fn various_f_values_test() { let p1 = Point::new(10.0, 20.0); let p2 = Point::new(125.0, 25.0); + #[allow(deprecated)] let i20 = p1.haversine_intermediate(&p2, 0.2); + #[allow(deprecated)] let i50 = p1.haversine_intermediate(&p2, 0.5); + #[allow(deprecated)] let i80 = p1.haversine_intermediate(&p2, 0.8); let i20_should = Point::new(29.83519, 29.94841); let i50_should = Point::new(65.87471, 37.72201); @@ -204,6 +120,7 @@ mod test { fn should_be_north_pole_test() { let p1 = Point::new(0.0, 10.0); let p2 = Point::new(180.0, 10.0); + #[allow(deprecated)] let i50 = p1.haversine_intermediate(&p2, 0.5); let i50_should = Point::new(90.0, 90.0); assert_relative_eq!(i50.x(), i50_should.x(), epsilon = 1.0e-6); @@ -215,6 +132,7 @@ mod test { let p1 = Point::new(30.0, 40.0); let p2 = Point::new(40.0, 50.0); let max_dist = 1500000.0; // meters + #[allow(deprecated)] let route = p1.haversine_intermediate_fill(&p2, max_dist, true); assert_eq!(route, vec![p1, p2]); } @@ -224,15 +142,14 @@ mod test { let p1 = Point::new(30.0, 40.0); let p2 = Point::new(40.0, 50.0); let max_dist = 1000000.0; // meters + #[allow(deprecated)] let i50 = p1.clone().haversine_intermediate(&p2, 0.5); - assert_eq!( - p1.haversine_intermediate_fill(&p2, max_dist, true), - vec![p1, i50, p2] - ); - assert_eq!( - p1.haversine_intermediate_fill(&p2, max_dist, false), - vec![i50] - ); + #[allow(deprecated)] + let fill = p1.haversine_intermediate_fill(&p2, max_dist, true); + assert_eq!(fill, vec![p1, i50, p2]); + #[allow(deprecated)] + let fill = p1.haversine_intermediate_fill(&p2, max_dist, false); + assert_eq!(fill, vec![i50]); } #[test] @@ -240,9 +157,13 @@ mod test { let p1 = Point::new(30.0, 40.0); let p2 = Point::new(40.0, 50.0); let max_dist = 400000.0; // meters + #[allow(deprecated)] let i25 = p1.clone().haversine_intermediate(&p2, 0.25); + #[allow(deprecated)] let i50 = p1.clone().haversine_intermediate(&p2, 0.5); + #[allow(deprecated)] let i75 = p1.clone().haversine_intermediate(&p2, 0.75); + #[allow(deprecated)] let route = p1.haversine_intermediate_fill(&p2, max_dist, true); assert_eq!(route, vec![p1, i25, i50, i75, p2]); } diff --git a/geo/src/algorithm/haversine_length.rs b/geo/src/algorithm/haversine_length.rs index b92f6c508c..573e21547e 100644 --- a/geo/src/algorithm/haversine_length.rs +++ b/geo/src/algorithm/haversine_length.rs @@ -1,7 +1,7 @@ use num_traits::FromPrimitive; -use crate::HaversineDistance; use crate::{CoordFloat, Line, LineString, MultiLineString}; +use crate::{Distance, Haversine}; /// Determine the length of a geometry using the [haversine formula]. /// @@ -47,7 +47,7 @@ where { fn haversine_length(&self) -> T { let (start, end) = self.points(); - start.haversine_distance(&end) + Haversine::distance(start, end) } } diff --git a/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs b/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs index a9b924a973..fbcfc7804b 100644 --- a/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs +++ b/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs @@ -1,5 +1,6 @@ use super::super::{Bearing, Destination, Distance, InterpolatePoint}; use crate::Point; +use geographiclib_rs::{DirectGeodesic, InverseGeodesic}; /// An ellipsoidal model of the earth, using methods given by [Karney (2013)]. /// @@ -36,7 +37,13 @@ impl Bearing for Geodesic { /// [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 + let (azi1, _, _) = geographiclib_rs::Geodesic::wgs84().inverse( + origin.y(), + origin.x(), + destination.y(), + destination.x(), + ); + (azi1 + 360.0) % 360.0 } } @@ -75,7 +82,9 @@ impl Destination for Geodesic { /// [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) + let (lat, lon) = + geographiclib_rs::Geodesic::wgs84().direct(origin.y(), origin.x(), bearing, distance); + Point::new(lon, lat) } } @@ -112,7 +121,12 @@ impl Distance, Point> for Geodesic { /// [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) + geographiclib_rs::Geodesic::wgs84().inverse( + origin.y(), + origin.x(), + destination.y(), + destination.x(), + ) } } @@ -153,11 +167,12 @@ impl InterpolatePoint for Geodesic { end: Point, ratio_from_start: f64, ) -> Point { - crate::algorithm::GeodesicIntermediate::geodesic_intermediate( - &start, - &end, - ratio_from_start, - ) + 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) } /// Interpolates `Point`s along a [geodesic line] between `start` and `end`. @@ -180,13 +195,35 @@ impl InterpolatePoint for Geodesic { max_distance: f64, include_ends: bool, ) -> impl Iterator> { - crate::algorithm::GeodesicIntermediate::geodesic_intermediate_fill( - &start, - &end, - max_distance, - include_ends, - ) - .into_iter() + let g = geographiclib_rs::Geodesic::wgs84(); + let (total_distance, azi1, _azi2, _a12) = g.inverse(start.y(), start.x(), end.y(), end.x()); + + if total_distance <= max_distance { + return if include_ends { + vec![start, end].into_iter() + } else { + vec![].into_iter() + }; + } + + let number_of_points = (total_distance / max_distance).ceil(); + let interval = 1.0 / number_of_points; + + let mut current_step = interval; + let mut points = if include_ends { vec![start] } else { vec![] }; + + while current_step < 1.0 { + let (lat2, lon2) = g.direct(start.y(), start.x(), azi1, total_distance * current_step); + let point = Point::new(lon2, lat2); + points.push(point); + current_step += interval; + } + + if include_ends { + points.push(end); + } + + points.into_iter() } } diff --git a/geo/src/algorithm/line_measures/metric_spaces/haversine.rs b/geo/src/algorithm/line_measures/metric_spaces/haversine.rs index 50186a2a5c..869d6ad25b 100644 --- a/geo/src/algorithm/line_measures/metric_spaces/haversine.rs +++ b/geo/src/algorithm/line_measures/metric_spaces/haversine.rs @@ -1,7 +1,8 @@ use num_traits::FromPrimitive; use super::super::{Bearing, Destination, Distance, InterpolatePoint}; -use crate::{CoordFloat, Point}; +use crate::utils::normalize_longitude; +use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS}; /// A spherical model of the earth using the [haversine formula]. /// @@ -47,8 +48,14 @@ impl Bearing for Haversine { 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 + let (lng_a, lat_a) = (origin.x().to_radians(), origin.y().to_radians()); + let (lng_b, lat_b) = (destination.x().to_radians(), destination.y().to_radians()); + let delta_lng = lng_b - lng_a; + let s = lat_b.cos() * delta_lng.sin(); + let c = lat_a.cos() * lat_b.sin() - lat_a.sin() * lat_b.cos() * delta_lng.cos(); + + let degrees = F::atan2(s, c).to_degrees(); + (degrees + three_sixty) % three_sixty } } @@ -82,7 +89,20 @@ impl Destination for Haversine { /// /// [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) + 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 lat = + { center_lat.sin() * rad.cos() + center_lat.cos() * rad.sin() * bearing_rad.cos() } + .asin(); + let lng = { bearing_rad.sin() * rad.sin() * center_lat.cos() } + .atan2(rad.cos() - center_lat.sin() * lat.sin()) + + center_lng; + + Point::new(normalize_longitude(lng.to_degrees()), lat.to_degrees()) } } @@ -119,7 +139,15 @@ impl Distance, Point> for Haversin /// /// [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula fn distance(origin: Point, destination: Point) -> F { - crate::algorithm::HaversineDistance::haversine_distance(&origin, &destination) + let two = F::one() + F::one(); + let theta1 = origin.y().to_radians(); + let theta2 = destination.y().to_radians(); + let delta_theta = (destination.y() - origin.y()).to_radians(); + let delta_lambda = (destination.x() - origin.x()).to_radians(); + let a = (delta_theta / two).sin().powi(2) + + theta1.cos() * theta2.cos() * (delta_lambda / two).sin().powi(2); + let c = two * a.sqrt().asin(); + F::from(MEAN_EARTH_RADIUS).unwrap() * c } } @@ -151,11 +179,8 @@ impl InterpolatePoint for Haversine { /// /// [great circle]: https://en.wikipedia.org/wiki/Great_circle fn point_at_ratio_between(start: Point, end: Point, ratio_from_start: F) -> Point { - crate::algorithm::HaversineIntermediate::haversine_intermediate( - &start, - &end, - ratio_from_start, - ) + let calculation = HaversineIntermediateFillCalculation::new(start, end); + calculation.point_at_ratio(ratio_from_start) } /// Interpolates `Point`s along a [great circle] between `start` and `end`. @@ -174,13 +199,114 @@ impl InterpolatePoint for Haversine { max_distance: F, include_ends: bool, ) -> impl Iterator> { - crate::algorithm::HaversineIntermediate::haversine_intermediate_fill( - &start, - &end, - max_distance, - include_ends, - ) - .into_iter() + let calculation = HaversineIntermediateFillCalculation::new(start, end); + let HaversineIntermediateFillCalculation { d, .. } = calculation; + + let total_distance = d * F::from(MEAN_EARTH_RADIUS).unwrap(); + + if total_distance <= max_distance { + return if include_ends { + vec![start, end].into_iter() + } else { + vec![].into_iter() + }; + } + + let number_of_points = (total_distance / max_distance).ceil(); + let interval = F::one() / number_of_points; + + let mut current_step = interval; + let mut points = if include_ends { vec![start] } else { vec![] }; + + while current_step < F::one() { + let point = calculation.point_at_ratio(current_step); + points.push(point); + current_step = current_step + interval; + } + + if include_ends { + points.push(end); + } + + points.into_iter() + } +} + +#[allow(clippy::many_single_char_names)] +struct HaversineIntermediateFillCalculation { + d: T, + n: T, + o: T, + p: T, + q: T, + r: T, + s: T, +} + +impl HaversineIntermediateFillCalculation { + #[allow(clippy::many_single_char_names)] + fn new(p1: Point, p2: Point) -> Self { + let one = T::one(); + let two = one + one; + + let lat1 = p1.y().to_radians(); + let lon1 = p1.x().to_radians(); + let lat2 = p2.y().to_radians(); + let lon2 = p2.x().to_radians(); + + let (lat1_sin, lat1_cos) = lat1.sin_cos(); + let (lat2_sin, lat2_cos) = lat2.sin_cos(); + let (lon1_sin, lon1_cos) = lon1.sin_cos(); + let (lon2_sin, lon2_cos) = lon2.sin_cos(); + + let m = lat1_cos * lat2_cos; + + let n = lat1_cos * lon1_cos; + let o = lat2_cos * lon2_cos; + let p = lat1_cos * lon1_sin; + let q = lat2_cos * lon2_sin; + + let k = + (((lat1 - lat2) / two).sin().powi(2) + m * ((lon1 - lon2) / two).sin().powi(2)).sqrt(); + + let d = two * k.asin(); + + Self { + d, + n, + o, + p, + q, + r: lat1_sin, + s: lat2_sin, + } + } + + #[allow(clippy::many_single_char_names)] + fn point_at_ratio(&self, ratio_from_start: T) -> Point { + let one = T::one(); + + let HaversineIntermediateFillCalculation { + d, + n, + o, + p, + q, + r, + s, + } = *self; + + let a = ((one - ratio_from_start) * d).sin() / d.sin(); + let b = (ratio_from_start * d).sin() / d.sin(); + + let x = a * n + b * o; + let y = a * p + b * q; + let z = a * r + b * s; + + let lat = z.atan2(x.hypot(y)); + let lon = y.atan2(x); + + Point::new(lon.to_degrees(), lat.to_degrees()) } } diff --git a/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs b/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs index a674a376bd..6c550311df 100644 --- a/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs +++ b/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs @@ -1,7 +1,8 @@ use num_traits::FromPrimitive; use super::super::{Bearing, Destination, Distance, InterpolatePoint}; -use crate::{CoordFloat, Point}; +use crate::rhumb::RhumbCalculations; +use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS}; /// Provides [rhumb line] (a.k.a. loxodrome) geometry operations. A rhumb line appears as a straight /// line on a Mercator projection map. @@ -47,7 +48,10 @@ impl Bearing for Rhumb { /// 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) + let three_sixty = F::from(360.0f64).unwrap(); + + let calculations = RhumbCalculations::new(&origin, &destination); + (calculations.theta().to_degrees() + three_sixty) % three_sixty } } @@ -76,7 +80,12 @@ impl Destination for Rhumb { /// /// [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) + let delta = distance / F::from(MEAN_EARTH_RADIUS).unwrap(); // angular distance in radians + let lambda1 = origin.x().to_radians(); + let phi1 = origin.y().to_radians(); + let theta = bearing.to_radians(); + + crate::algorithm::rhumb::calculate_destination(delta, lambda1, phi1, theta) } } @@ -110,7 +119,8 @@ impl Distance, Point> for Rhumb { /// /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line fn distance(origin: Point, destination: Point) -> F { - crate::algorithm::RhumbDistance::rhumb_distance(&origin, &destination) + let calculations = RhumbCalculations::new(&origin, &destination); + calculations.delta() * F::from(MEAN_EARTH_RADIUS).unwrap() } } @@ -142,7 +152,8 @@ impl InterpolatePoint for Rhumb { /// /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line fn point_at_ratio_between(start: Point, end: Point, ratio_from_start: F) -> Point { - crate::algorithm::RhumbIntermediate::rhumb_intermediate(&start, &end, ratio_from_start) + let calculations = RhumbCalculations::new(&start, &end); + calculations.intermediate(ratio_from_start) } /// Interpolates `Point`s along a [rhumb line] between `start` and `end`. @@ -160,13 +171,11 @@ impl InterpolatePoint for Rhumb { max_distance: F, include_ends: bool, ) -> impl Iterator> { - crate::algorithm::RhumbIntermediate::rhumb_intermediate_fill( - &start, - &end, - max_distance, - include_ends, - ) - .into_iter() + let max_delta = max_distance / F::from(MEAN_EARTH_RADIUS).unwrap(); + let calculations = RhumbCalculations::new(&start, &end); + calculations + .intermediate_fill(max_delta, include_ends) + .into_iter() } } diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index 8edaa9e308..24c146b6f7 100644 --- a/geo/src/algorithm/mod.rs +++ b/geo/src/algorithm/mod.rs @@ -100,10 +100,12 @@ pub use geodesic_bearing::GeodesicBearing; /// Returns a new Point using a distance and bearing on a geodesic. pub mod geodesic_destination; +#[allow(deprecated)] pub use geodesic_destination::GeodesicDestination; /// Calculate the Geodesic distance between two `Point`s. pub mod geodesic_distance; +#[allow(deprecated)] pub use geodesic_distance::GeodesicDistance; /// Calculate the Geodesic area and perimeter of polygons. @@ -112,6 +114,7 @@ pub use geodesic_area::GeodesicArea; /// Calculate a new `Point` lying on a Geodesic arc between two `Point`s. pub mod geodesic_intermediate; +#[allow(deprecated)] pub use geodesic_intermediate::GeodesicIntermediate; /// Calculate the Geodesic length of a line. @@ -124,18 +127,22 @@ pub use hausdorff_distance::HausdorffDistance; /// Calculate the bearing to another `Point`, in degrees. pub mod haversine_bearing; +#[allow(deprecated)] pub use haversine_bearing::HaversineBearing; /// Calculate a destination `Point`, given a distance and a bearing. pub mod haversine_destination; +#[allow(deprecated)] pub use haversine_destination::HaversineDestination; /// Calculate the Haversine distance between two `Geometries`. pub mod haversine_distance; +#[allow(deprecated)] pub use haversine_distance::HaversineDistance; /// Calculate a new `Point` lying on a Great Circle arc between two `Point`s. pub mod haversine_intermediate; +#[allow(deprecated)] pub use haversine_intermediate::HaversineIntermediate; /// Calculate the Haversine length of a Line. @@ -291,4 +298,6 @@ pub use monotone::{monotone_subdivision, MonoPoly, MonotonicPolygons}; /// Rhumb-line-related algorithms and utils pub mod rhumb; -pub use rhumb::{RhumbBearing, RhumbDestination, RhumbDistance, RhumbIntermediate, RhumbLength}; +pub use rhumb::RhumbLength; +#[allow(deprecated)] +pub use rhumb::{RhumbBearing, RhumbDestination, RhumbDistance, RhumbIntermediate}; diff --git a/geo/src/algorithm/rhumb/bearing.rs b/geo/src/algorithm/rhumb/bearing.rs index 4c3be743f8..c67c5421da 100644 --- a/geo/src/algorithm/rhumb/bearing.rs +++ b/geo/src/algorithm/rhumb/bearing.rs @@ -1,14 +1,15 @@ use num_traits::FromPrimitive; -use crate::{CoordFloat, Point}; - -use super::RhumbCalculations; +use crate::{Bearing, CoordFloat, Point, Rhumb}; +#[deprecated( + since = "0.29.0", + note = "Please use the `Rhumb::bearing` method from the `Bearing` trait instead" +)] /// Returns the bearing to another Point in degrees. /// /// Bullock, R.: Great Circle Distances and Bearings Between Two Locations, 2007. /// () - pub trait RhumbBearing { /// Returns the bearing to another Point in degrees along a [rhumb line], where North is 0° and East is 90°. /// @@ -16,11 +17,13 @@ pub trait RhumbBearing { /// /// ``` /// # use approx::assert_relative_eq; + /// # #[allow(deprecated)] /// use geo::RhumbBearing; /// use geo::Point; /// /// let p_1 = Point::new(9.177789688110352, 48.776781529534965); /// let p_2 = Point::new(9.274348757829898, 48.84037308229984); + /// # #[allow(deprecated)] /// let bearing = p_1.rhumb_bearing(p_2); /// assert_relative_eq!(bearing, 45., epsilon = 1.0e-6); /// ``` @@ -28,28 +31,29 @@ pub trait RhumbBearing { fn rhumb_bearing(&self, point: Point) -> T; } +#[allow(deprecated)] impl RhumbBearing for Point where T: CoordFloat + FromPrimitive, { fn rhumb_bearing(&self, point: Point) -> T { - let three_sixty = T::from(360.0f64).unwrap(); - - let calculations = RhumbCalculations::new(self, &point); - (calculations.theta().to_degrees() + three_sixty) % three_sixty + Rhumb::bearing(*self, point) } } #[cfg(test)] mod test { use crate::point; + #[allow(deprecated)] use crate::RhumbBearing; + #[allow(deprecated)] use crate::RhumbDestination; #[test] fn north_bearing() { let p_1 = point!(x: 9., y: 47.); let p_2 = point!(x: 9., y: 48.); + #[allow(deprecated)] let bearing = p_1.rhumb_bearing(p_2); assert_relative_eq!(bearing, 0.); } @@ -58,6 +62,7 @@ mod test { fn equatorial_east_bearing() { let p_1 = point!(x: 9., y: 0.); let p_2 = point!(x: 10., y: 0.); + #[allow(deprecated)] let bearing = p_1.rhumb_bearing(p_2); assert_relative_eq!(bearing, 90.); } @@ -67,6 +72,7 @@ mod test { let p_1 = point!(x: 9., y: 10.); let p_2 = point!(x: 18.131938299366652, y: 10.); + #[allow(deprecated)] let bearing = p_1.rhumb_bearing(p_2); assert_relative_eq!(bearing, 90.); } @@ -75,6 +81,7 @@ mod test { fn northeast_bearing() { let p_1 = point!(x: 9.177789688110352f64, y: 48.776781529534965); let p_2 = point!(x: 9.274348757829898, y: 48.84037308229984); + #[allow(deprecated)] let bearing = p_1.rhumb_bearing(p_2); assert_relative_eq!(bearing, 45., epsilon = 1.0e-6); } @@ -82,8 +89,10 @@ mod test { #[test] fn consistent_with_destination() { let p_1 = point!(x: 9.177789688110352f64, y: 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.rhumb_destination(45., 10000.); + #[allow(deprecated)] let b_1 = p_1.rhumb_bearing(p_2); assert_relative_eq!(b_1, 45., epsilon = 1.0e-6); } diff --git a/geo/src/algorithm/rhumb/destination.rs b/geo/src/algorithm/rhumb/destination.rs index a9b9be85e0..a4e1bdbc1b 100644 --- a/geo/src/algorithm/rhumb/destination.rs +++ b/geo/src/algorithm/rhumb/destination.rs @@ -1,8 +1,10 @@ -use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS}; +use crate::{CoordFloat, Destination, Point, Rhumb}; use num_traits::FromPrimitive; -use super::calculate_destination; - +#[deprecated( + since = "0.29.0", + note = "Please use the `Rhumb::destination` method from the `Destination` trait instead" +)] /// Returns the destination Point having travelled the given distance along a [rhumb line] /// from the origin geometry with the given bearing /// @@ -20,10 +22,12 @@ pub trait RhumbDestination { /// # Examples /// /// ``` + /// # #[allow(deprecated)] /// use geo::RhumbDestination; /// use geo::Point; /// /// let p_1 = Point::new(9.177789688110352, 48.776781529534965); + /// # #[allow(deprecated)] /// let p_2 = p_1.rhumb_destination(45., 10000.); /// assert_eq!(p_2, Point::new(9.274348757829898, 48.84037308229984)) /// ``` @@ -31,31 +35,30 @@ pub trait RhumbDestination { fn rhumb_destination(&self, bearing: T, distance: T) -> Point; } +#[allow(deprecated)] impl RhumbDestination for Point where T: CoordFloat + FromPrimitive, { fn rhumb_destination(&self, bearing: T, distance: T) -> Point { - let delta = distance / T::from(MEAN_EARTH_RADIUS).unwrap(); // angular distance in radians - let lambda1 = self.x().to_radians(); - let phi1 = self.y().to_radians(); - let theta = bearing.to_radians(); - - calculate_destination(delta, lambda1, phi1, theta) + Rhumb::destination(*self, bearing, distance) } } #[cfg(test)] mod test { use super::*; + #[allow(deprecated)] use crate::RhumbDistance; use num_traits::pow; #[test] fn returns_a_new_point() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.rhumb_destination(45., 10000.); assert_eq!(p_2, Point::new(9.274348757829898, 48.84037308229984)); + #[allow(deprecated)] let distance = p_1.rhumb_distance(&p_2); assert_relative_eq!(distance, 10000., epsilon = 1.0e-6) } @@ -63,9 +66,12 @@ mod test { #[test] fn direct_and_indirect_destinations_are_close() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.rhumb_destination(45., 10000.); let square_edge = { pow(10000., 2) / 2f64 }.sqrt(); + #[allow(deprecated)] let p_3 = p_1.rhumb_destination(0., square_edge); + #[allow(deprecated)] let p_4 = p_3.rhumb_destination(90., square_edge); assert_relative_eq!(p_4, p_2, epsilon = 1.0e-3); } @@ -73,6 +79,7 @@ mod test { #[test] fn bearing_zero_is_north() { let p_1 = Point::new(9.177789688110352, 48.776781529534965); + #[allow(deprecated)] let p_2 = p_1.rhumb_destination(0., 1000.); assert_relative_eq!(p_1.x(), p_2.x(), epsilon = 1.0e-6); assert!(p_2.y() > p_1.y()) diff --git a/geo/src/algorithm/rhumb/distance.rs b/geo/src/algorithm/rhumb/distance.rs index a0fc8fd6b3..0032b288ac 100644 --- a/geo/src/algorithm/rhumb/distance.rs +++ b/geo/src/algorithm/rhumb/distance.rs @@ -1,8 +1,10 @@ -use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS}; +use crate::{CoordFloat, Distance, Point, Rhumb}; use num_traits::FromPrimitive; -use super::RhumbCalculations; - +#[deprecated( + since = "0.29.0", + note = "Please use the `Rhumb::distance` method from the `Distance` trait instead" +)] /// Determine the distance between two geometries along a [rhumb line]. /// /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line @@ -28,6 +30,7 @@ pub trait RhumbDistance { /// // London /// let p2 = point!(x: -0.1278f64, y: 51.5074f64); /// + /// # #[allow(deprecated)] /// let distance = p1.rhumb_distance(&p2); /// /// assert_eq!( @@ -40,41 +43,38 @@ pub trait RhumbDistance { fn rhumb_distance(&self, rhs: &Rhs) -> T; } +#[allow(deprecated)] impl RhumbDistance> for Point where T: CoordFloat + FromPrimitive, { fn rhumb_distance(&self, rhs: &Point) -> T { - let calculations = RhumbCalculations::new(self, rhs); - calculations.delta() * T::from(MEAN_EARTH_RADIUS).unwrap() + Rhumb::distance(*self, *rhs) } } #[cfg(test)] mod test { use crate::Point; + #[allow(deprecated)] use crate::RhumbDistance; #[test] fn distance1_test() { let a = Point::new(0., 0.); let b = Point::new(1., 0.); - assert_relative_eq!( - a.rhumb_distance(&b), - 111195.0802335329_f64, - epsilon = 1.0e-6 - ); + #[allow(deprecated)] + let distance = a.rhumb_distance(&b); + assert_relative_eq!(distance, 111195.0802335329_f64, epsilon = 1.0e-6); } #[test] fn distance2_test() { let a = Point::new(-72.1235, 42.3521); let b = Point::new(72.1260, 70.612); - assert_relative_eq!( - a.rhumb_distance(&b), - 8903668.508603323_f64, - epsilon = 1.0e-6 - ); + #[allow(deprecated)] + let distance = a.rhumb_distance(&b); + assert_relative_eq!(distance, 8903668.508603323_f64, epsilon = 1.0e-6); } #[test] @@ -82,11 +82,9 @@ mod test { // this input comes from issue #100 let a = Point::new(-77.036585, 38.897448); let b = Point::new(-77.009080, 38.889825); - assert_relative_eq!( - a.rhumb_distance(&b), - 2526.7031699343006_f64, - epsilon = 1.0e-6 - ); + #[allow(deprecated)] + let distance = a.rhumb_distance(&b); + assert_relative_eq!(distance, 2526.7031699343006_f64, epsilon = 1.0e-6); } #[test] @@ -94,6 +92,8 @@ mod test { // this input comes from issue #100 let a = Point::::new(-77.03658, 38.89745); let b = Point::::new(-77.00908, 38.889825); - assert_relative_eq!(a.rhumb_distance(&b), 2526.7273_f32, epsilon = 1.0e-6); + #[allow(deprecated)] + let distance = a.rhumb_distance(&b); + assert_relative_eq!(distance, 2526.7273_f32, epsilon = 1.0e-6); } } diff --git a/geo/src/algorithm/rhumb/intermediate.rs b/geo/src/algorithm/rhumb/intermediate.rs index 54651c6ec1..198be20abe 100644 --- a/geo/src/algorithm/rhumb/intermediate.rs +++ b/geo/src/algorithm/rhumb/intermediate.rs @@ -1,24 +1,33 @@ -use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS}; +use crate::{CoordFloat, InterpolatePoint, Point, Rhumb}; use num_traits::FromPrimitive; -use super::RhumbCalculations; - +#[deprecated( + since = "0.29.0", + note = "Please use the `InterpolatePoint` trait instead" +)] /// Returns a new Point along a rhumb line between two existing points - pub trait RhumbIntermediate { + #[deprecated( + since = "0.29.0", + note = "Please use `Rhumb::point_at_ratio_between` from the `InterpolatePoint` trait instead" + )] /// Returns a new Point along a [rhumb line] between two existing points. /// /// # Examples /// /// ```rust /// # use approx::assert_relative_eq; + /// # #[allow(deprecated)] /// use geo::RhumbIntermediate; /// use geo::Point; /// /// let p1 = Point::new(10.0, 20.0); /// let p2 = Point::new(125.0, 25.0); + /// # #[allow(deprecated)] /// let i20 = p1.rhumb_intermediate(&p2, 0.2); + /// # #[allow(deprecated)] /// let i50 = p1.rhumb_intermediate(&p2, 0.5); + /// # #[allow(deprecated)] /// let i80 = p1.rhumb_intermediate(&p2, 0.8); /// let i20_should = Point::new(32.7, 21.0); /// let i50_should = Point::new(67.0, 22.5); @@ -31,8 +40,12 @@ pub trait RhumbIntermediate { /// assert_relative_eq!(i80.y(), i80_should.y(), epsilon = 0.2); /// ``` /// [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line - fn rhumb_intermediate(&self, other: &Point, f: T) -> Point; + + #[deprecated( + since = "0.29.0", + note = "Please use `Rhumb::points_along_line` from the `InterpolatePoint` trait instead" + )] fn rhumb_intermediate_fill( &self, other: &Point, @@ -41,13 +54,13 @@ pub trait RhumbIntermediate { ) -> Vec>; } +#[allow(deprecated)] impl RhumbIntermediate for Point where T: CoordFloat + FromPrimitive, { fn rhumb_intermediate(&self, other: &Point, f: T) -> Point { - let calculations = RhumbCalculations::new(self, other); - calculations.intermediate(f) + Rhumb::point_at_ratio_between(*self, *other, f) } fn rhumb_intermediate_fill( @@ -56,22 +69,23 @@ where max_dist: T, include_ends: bool, ) -> Vec> { - let max_delta = max_dist / T::from(MEAN_EARTH_RADIUS).unwrap(); - let calculations = RhumbCalculations::new(self, other); - calculations.intermediate_fill(max_delta, include_ends) + Rhumb::points_along_line(*self, *other, max_dist, include_ends).collect() } } #[cfg(test)] mod test { use super::*; + #[allow(deprecated)] use crate::RhumbIntermediate; #[test] fn f_is_zero_or_one_test() { let p1 = Point::new(10.0, 20.0); let p2 = Point::new(15.0, 25.0); + #[allow(deprecated)] let i0 = p1.rhumb_intermediate(&p2, 0.0); + #[allow(deprecated)] let i100 = p1.rhumb_intermediate(&p2, 1.0); assert_relative_eq!(i0.x(), p1.x(), epsilon = 1.0e-6); assert_relative_eq!(i0.y(), p1.y(), epsilon = 1.0e-6); @@ -83,8 +97,11 @@ mod test { fn various_f_values_test() { let p1 = Point::new(10.0, 20.0); let p2 = Point::new(125.0, 25.0); + #[allow(deprecated)] let i20 = p1.rhumb_intermediate(&p2, 0.2); + #[allow(deprecated)] let i50 = p1.rhumb_intermediate(&p2, 0.5); + #[allow(deprecated)] let i80 = p1.rhumb_intermediate(&p2, 0.8); let i20_should = Point::new(32.6766, 21.0); let i50_should = Point::new(66.9801, 22.5); @@ -101,6 +118,7 @@ mod test { fn should_be_straight_across_test() { let p1 = Point::new(0.0, 10.0); let p2 = Point::new(180.0, 10.0); + #[allow(deprecated)] let i50 = p1.rhumb_intermediate(&p2, 0.5); let i50_should = Point::new(90.0, 10.0); assert_relative_eq!(i50.x(), i50_should.x(), epsilon = 1.0e-6); @@ -113,6 +131,7 @@ mod test { let p2 = Point::new(40.0, 50.0); let max_dist = 1500000.0; // meters let include_ends = true; + #[allow(deprecated)] let route = p1.rhumb_intermediate_fill(&p2, max_dist, include_ends); assert_eq!(route, vec![p1, p2]); } @@ -123,7 +142,9 @@ mod test { let p2 = Point::new(40.0, 50.0); let max_dist = 1000000.0; // meters let include_ends = true; + #[allow(deprecated)] let i50 = p1.clone().rhumb_intermediate(&p2, 0.5); + #[allow(deprecated)] let route = p1.rhumb_intermediate_fill(&p2, max_dist, include_ends); assert_eq!(route, vec![p1, i50, p2]); } @@ -134,9 +155,13 @@ mod test { let p2 = Point::new(40.0, 50.0); let max_dist = 400000.0; // meters let include_ends = true; + #[allow(deprecated)] let i25 = p1.clone().rhumb_intermediate(&p2, 0.25); + #[allow(deprecated)] let i50 = p1.clone().rhumb_intermediate(&p2, 0.5); + #[allow(deprecated)] let i75 = p1.clone().rhumb_intermediate(&p2, 0.75); + #[allow(deprecated)] let route = p1.rhumb_intermediate_fill(&p2, max_dist, include_ends); assert_eq!(route, vec![p1, i25, i50, i75, p2]); } diff --git a/geo/src/algorithm/rhumb/length.rs b/geo/src/algorithm/rhumb/length.rs index ce544ece42..8f761fe93b 100644 --- a/geo/src/algorithm/rhumb/length.rs +++ b/geo/src/algorithm/rhumb/length.rs @@ -1,7 +1,6 @@ use num_traits::FromPrimitive; -use crate::RhumbDistance; -use crate::{CoordFloat, Line, LineString, MultiLineString}; +use crate::{CoordFloat, Distance, Line, LineString, MultiLineString, Rhumb}; /// Determine the length of a geometry assuming each segment is a [rhumb line]. /// @@ -47,7 +46,7 @@ where { fn rhumb_length(&self) -> T { let (start, end) = self.points(); - start.rhumb_distance(&end) + Rhumb::distance(start, end) } } diff --git a/geo/src/algorithm/rhumb/mod.rs b/geo/src/algorithm/rhumb/mod.rs index 7ffe9f8b7d..8c581b9400 100644 --- a/geo/src/algorithm/rhumb/mod.rs +++ b/geo/src/algorithm/rhumb/mod.rs @@ -10,21 +10,25 @@ use crate::{point, utils::normalize_longitude, CoordFloat, Point}; use num_traits::FromPrimitive; mod distance; +#[allow(deprecated)] pub use distance::RhumbDistance; mod bearing; +#[allow(deprecated)] pub use bearing::RhumbBearing; mod destination; +#[allow(deprecated)] pub use destination::RhumbDestination; mod intermediate; +#[allow(deprecated)] pub use intermediate::RhumbIntermediate; mod length; pub use length::RhumbLength; -struct RhumbCalculations { +pub(crate) struct RhumbCalculations { from: Point, to: Point, phi1: T, @@ -34,7 +38,7 @@ struct RhumbCalculations { } impl RhumbCalculations { - fn new(from: &Point, to: &Point) -> Self { + pub(crate) fn new(from: &Point, to: &Point) -> Self { let pi = T::from(std::f64::consts::PI).unwrap(); let two = T::one() + T::one(); let four = two + two; @@ -63,11 +67,11 @@ impl RhumbCalculations { } } - fn theta(&self) -> T { + pub(crate) fn theta(&self) -> T { self.delta_lambda.atan2(self.delta_psi) } - fn delta(&self) -> T { + pub(crate) fn delta(&self) -> T { let threshold = T::from(10.0e-12).unwrap(); let q = if self.delta_psi > threshold { self.delta_phi / self.delta_psi @@ -78,14 +82,14 @@ impl RhumbCalculations { (self.delta_phi * self.delta_phi + q * q * self.delta_lambda * self.delta_lambda).sqrt() } - fn intermediate(&self, fraction: T) -> Point { + pub(crate) fn intermediate(&self, fraction: T) -> Point { let delta = fraction * self.delta(); let theta = self.theta(); let lambda1 = self.from.x().to_radians(); calculate_destination(delta, lambda1, self.phi1, theta) } - fn intermediate_fill(&self, max_delta: T, include_ends: bool) -> Vec> { + pub(crate) fn intermediate_fill(&self, max_delta: T, include_ends: bool) -> Vec> { let theta = self.theta(); let lambda1 = self.from.x().to_radians(); @@ -124,7 +128,7 @@ impl RhumbCalculations { } } -fn calculate_destination( +pub(crate) fn calculate_destination( delta: T, lambda1: T, phi1: T,