Skip to content

Commit

Permalink
Geodesic, Rhumb, Haversine line measures: deprecate legacy traits aft…
Browse files Browse the repository at this point in the history
…er moving implementation to modern traits

Moved some implementation details from HaversineIntermediateFill to
Haversine.
  • Loading branch information
michaelkirk committed Oct 8, 2024
1 parent ccf2ffd commit 414de5a
Show file tree
Hide file tree
Showing 24 changed files with 509 additions and 339 deletions.
6 changes: 2 additions & 4 deletions geo/benches/geodesic_distance.rs
Original file line number Diff line number Diff line change
@@ -1,15 +1,13 @@
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| {
let a = geo::Point::new(17.107558, 48.148636);
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)));
});
});
}
Expand Down
14 changes: 7 additions & 7 deletions geo/src/algorithm/cross_track_distance.rs
Original file line number Diff line number Diff line change
@@ -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;

Expand Down Expand Up @@ -43,9 +43,9 @@ where
{
fn cross_track_distance(&self, line_point_a: &Point<T>, line_point_b: &Point<T>) -> 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()
}
Expand All @@ -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() {
Expand Down Expand Up @@ -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
);
}
Expand Down
7 changes: 3 additions & 4 deletions geo/src/algorithm/densify_haversine.rs
Original file line number Diff line number Diff line change
@@ -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.
///
Expand Down Expand Up @@ -52,8 +52,7 @@ fn densify_line<T: CoordFloat + FromPrimitive>(
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);
}
}
Expand Down
11 changes: 11 additions & 0 deletions geo/src/algorithm/geodesic_bearing.rs
Original file line number Diff line number Diff line change
Expand Up @@ -8,17 +8,23 @@ use geographiclib_rs::{Geodesic, InverseGeodesic};
///
/// [Karney (2013)]: https://arxiv.org/pdf/1109.4448.pdf
pub trait GeodesicBearing<T: CoordNum> {
#[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);
/// ```
Expand Down Expand Up @@ -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.);
}
Expand All @@ -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.);
}
Expand All @@ -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);
Expand Down
17 changes: 14 additions & 3 deletions geo/src/algorithm/geodesic_destination.rs
Original file line number Diff line number Diff line change
@@ -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)].
Expand All @@ -18,6 +22,7 @@ pub trait GeodesicDestination<T: CoordNum> {
/// # Examples
///
/// ```rust
/// # #[allow(deprecated)]
/// use geo::GeodesicDestination;
/// use geo::Point;
///
Expand All @@ -26,6 +31,7 @@ pub trait GeodesicDestination<T: CoordNum> {
/// 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);
Expand All @@ -34,21 +40,23 @@ pub trait GeodesicDestination<T: CoordNum> {
fn geodesic_destination(&self, bearing: T, distance: T) -> Point<T>;
}

#[allow(deprecated)]
impl GeodesicDestination<f64> for Point<f64> {
fn geodesic_destination(&self, bearing: f64, distance: f64) -> Point<f64> {
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!(
Expand All @@ -57,13 +65,15 @@ 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)
}

#[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())
Expand All @@ -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())
Expand Down
11 changes: 8 additions & 3 deletions geo/src/algorithm/geodesic_distance.rs
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -28,6 +31,7 @@ pub trait GeodesicDistance<T, Rhs = Self> {
/// // London
/// let p2 = point!(x: -0.1278, y: 51.5074);
///
/// # #[allow(deprecated)]
/// let distance = p1.geodesic_distance(&p2);
///
/// assert_eq!(
Expand All @@ -39,8 +43,9 @@ pub trait GeodesicDistance<T, Rhs = Self> {
fn geodesic_distance(&self, rhs: &Rhs) -> T;
}

#[allow(deprecated)]
impl GeodesicDistance<f64> for Point {
fn geodesic_distance(&self, rhs: &Point) -> f64 {
Geodesic::wgs84().inverse(self.y(), self.x(), rhs.y(), rhs.x())
Geodesic::distance(*self, *rhs)
}
}
67 changes: 28 additions & 39 deletions geo/src/algorithm/geodesic_intermediate.rs
Original file line number Diff line number Diff line change
@@ -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<T: CoordFloat> {
#[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);
Expand All @@ -25,6 +36,11 @@ pub trait GeodesicIntermediate<T: CoordFloat> {
/// assert_relative_eq!(i80, i80_should, epsilon = 1.0e-6);
/// ```
fn geodesic_intermediate(&self, other: &Point<T>, f: T) -> Point<T>;

#[deprecated(
since = "0.29.0",
note = "Please use `Geodesic::points_along_line` from the `InterpolatePoint` trait instead"
)]
fn geodesic_intermediate_fill(
&self,
other: &Point<T>,
Expand All @@ -33,15 +49,10 @@ pub trait GeodesicIntermediate<T: CoordFloat> {
) -> Vec<Point<T>>;
}

#[allow(deprecated)]
impl GeodesicIntermediate<f64> 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(
Expand All @@ -50,36 +61,7 @@ impl GeodesicIntermediate<f64> for Point {
max_dist: f64,
include_ends: bool,
) -> Vec<Point> {
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()
}
}

Expand All @@ -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);
Expand All @@ -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);
Expand All @@ -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]);
}
Expand Down
Loading

0 comments on commit 414de5a

Please sign in to comment.