geometry_strategy/geographic/
length.rs

1//! Geographic geodesic length / perimeter.
2//!
3//! Mirrors the per-CS dispatch arm of
4//! `boost/geometry/algorithms/length.hpp` reached when the geometry's
5//! CS family is `Geographic`. The Boost geographic length strategy is
6//! `strategies::length::geographic` from
7//! `boost/geometry/strategies/length/geographic.hpp`, which hands a
8//! geodesic inverse formula to the per-segment walk; we mirror that —
9//! [`GeographicLength`] / [`GeographicPerimeter`] are thin shells
10//! around the existing [`Andoyer`] distance
11//! strategy.
12
13use alloc::vec::Vec;
14
15use geometry_cs::{CoordinateSystem, GeographicFamily};
16use geometry_tag::SameAs;
17use geometry_trait::{Closure, Linestring, Point, Ring};
18
19use crate::distance::DistanceStrategy;
20use crate::geographic::Andoyer;
21use crate::length::LengthStrategy;
22
23/// Geographic length: sum of Andoyer geodesic distances between
24/// consecutive points.
25///
26/// Carries the same `spheroid` parameter as
27/// [`Andoyer`] so callers can pick WGS84 /
28/// GDA / custom. `Default::default()` produces
29/// `GeographicLength { andoyer: Andoyer::WGS84 }`.
30///
31/// Mirrors `boost::geometry::strategies::length::geographic<>` from
32/// `strategies/length/geographic.hpp`.
33#[derive(Debug, Default, Clone, Copy)]
34pub struct GeographicLength {
35    /// The inner Andoyer geodesic kernel summed over each segment.
36    pub andoyer: Andoyer,
37}
38
39/// Geographic perimeter — same as [`GeographicLength`] but adds the
40/// closing edge for an open ring.
41///
42/// Separate from [`GeographicLength`] because Rust's coherence rules
43/// cannot prove that a single type does not implement both
44/// [`Linestring`] and [`Ring`]; splitting the strategy keeps the
45/// per-tag dispatch disjoint at the impl level.
46#[derive(Debug, Default, Clone, Copy)]
47pub struct GeographicPerimeter {
48    /// The inner Andoyer geodesic kernel summed over each segment.
49    pub andoyer: Andoyer,
50}
51
52// ---- Linestring ------------------------------------------------------
53
54impl<L> LengthStrategy<L> for GeographicLength
55where
56    L: Linestring,
57    <L::Point as Point>::Cs: CoordinateSystem,
58    <<L::Point as Point>::Cs as CoordinateSystem>::Family: SameAs<GeographicFamily>,
59    Andoyer: DistanceStrategy<L::Point, L::Point, Out = <L::Point as Point>::Scalar>,
60{
61    type Out = <L::Point as Point>::Scalar;
62
63    #[inline]
64    fn length(&self, g: &L) -> Self::Out {
65        let pts: Vec<&L::Point> = g.points().collect();
66        let mut acc = <Self::Out as geometry_coords::CoordinateScalar>::ZERO;
67        for w in pts.windows(2) {
68            acc = acc + self.andoyer.distance(w[0], w[1]);
69        }
70        acc
71    }
72}
73
74// ---- Ring ------------------------------------------------------------
75
76impl<R> LengthStrategy<R> for GeographicPerimeter
77where
78    R: Ring,
79    <R::Point as Point>::Cs: CoordinateSystem,
80    <<R::Point as Point>::Cs as CoordinateSystem>::Family: SameAs<GeographicFamily>,
81    Andoyer: DistanceStrategy<R::Point, R::Point, Out = <R::Point as Point>::Scalar>,
82{
83    type Out = <R::Point as Point>::Scalar;
84
85    #[inline]
86    fn length(&self, g: &R) -> Self::Out {
87        let pts: Vec<&R::Point> = g.points().collect();
88        let mut acc = <Self::Out as geometry_coords::CoordinateScalar>::ZERO;
89        for w in pts.windows(2) {
90            acc = acc + self.andoyer.distance(w[0], w[1]);
91        }
92        if matches!(g.closure(), Closure::Open) {
93            if let (Some(first), Some(last)) = (pts.first(), pts.last()) {
94                acc = acc + self.andoyer.distance(last, first);
95            }
96        }
97        acc
98    }
99}
100
101#[cfg(all(test, feature = "std"))]
102mod tests {
103    //! Reference from `boost/geometry/test/algorithms/length/length_geo.cpp`.
104    //! Values are cross-checked against a direct [`Andoyer`] call on the
105    //! same segments (the length strategy reuses the distance strategy
106    //! verbatim, so the two must agree exactly).
107    #![allow(
108        clippy::float_cmp,
109        reason = "lengths are compared with an explicit absolute tolerance, not `==`"
110    )]
111
112    use super::{GeographicLength, GeographicPerimeter};
113    use crate::distance::DistanceStrategy;
114    use crate::geographic::Andoyer;
115    use crate::length::LengthStrategy;
116    use geometry_adapt::{Adapt, WithCs};
117    use geometry_cs::{Degree, Geographic};
118    use geometry_model::{Linestring, Ring};
119
120    type Gg = WithCs<Adapt<[f64; 2]>, Geographic<Degree>>;
121
122    #[inline]
123    fn deg(lon: f64, lat: f64) -> Gg {
124        WithCs::new(Adapt([lon, lat]))
125    }
126
127    /// A single-segment linestring's length equals a direct Andoyer
128    /// call — `(4, 52) → (3, 40)` on WGS84 ≈ 1336.040 km
129    /// (`andoyer.cpp:230-231`).
130    #[test]
131    fn ams_geodesic_length() {
132        let ls: Linestring<Gg> = Linestring(vec![deg(4.0, 52.0), deg(3.0, 40.0)]);
133        let got = GeographicLength::default().length(&ls);
134        let direct = Andoyer::WGS84.distance(&deg(4.0, 52.0), &deg(3.0, 40.0));
135        assert!((got - direct).abs() < 1e-6);
136        assert!((got / 1000.0 - 1_336.039_890).abs() < 0.01);
137    }
138
139    /// A three-point linestring sums two geodesic hops.
140    #[test]
141    fn three_point_geodesic() {
142        let ls: Linestring<Gg> = Linestring(vec![deg(0., 0.), deg(0., 10.), deg(0., 20.)]);
143        let got = GeographicLength::default().length(&ls);
144        let a = Andoyer::WGS84;
145        let expected =
146            a.distance(&deg(0., 0.), &deg(0., 10.)) + a.distance(&deg(0., 10.), &deg(0., 20.));
147        assert!((got - expected).abs() < 1e-6);
148    }
149
150    /// An open ring's perimeter adds the closing edge back to the first
151    /// vertex.
152    #[test]
153    fn open_ring_closes_itself() {
154        let mut r = Ring::<Gg, true, false>::new();
155        r.push(deg(0., 0.));
156        r.push(deg(10., 0.));
157        r.push(deg(10., 10.));
158        r.push(deg(0., 10.));
159        let got = GeographicPerimeter::default().length(&r);
160        let a = Andoyer::WGS84;
161        let expected = a.distance(&deg(0., 0.), &deg(10., 0.))
162            + a.distance(&deg(10., 0.), &deg(10., 10.))
163            + a.distance(&deg(10., 10.), &deg(0., 10.))
164            + a.distance(&deg(0., 10.), &deg(0., 0.));
165        assert!((got - expected).abs() < 1e-6);
166    }
167}
geometry_strategy/geographic/length.rs

geometry_strategy/geographic/
length.rs
