gistools/geometry/tools/lines/
length.rs

1use libm::{asin, cos, fmin, pow, sin, sqrt};
2use s2json::{
3    Feature, Geometry, GetXY, GetZ, MultiLineString, MultiLineString3D, MultiLineString3DGeometry,
4    MultiLineStringGeometry, MultiPoint, MultiPoint3D, MultiPoint3DGeometry, MultiPointGeometry,
5    MultiPolygon, MultiPolygon3D, MultiPolygon3DGeometry, MultiPolygonGeometry, Point, Point3D,
6    Point3DGeometry, PointGeometry, VectorFeature, VectorGeometry, VectorMultiLineString,
7    VectorMultiLineStringGeometry, VectorMultiPoint, VectorMultiPointGeometry, VectorMultiPolygon,
8    VectorMultiPolygonGeometry, VectorPoint, VectorPointGeometry,
9};
10
11/// Get the total euclidean distance of a line or lines
12///
13/// This trait is implemented for:
14/// - [`Feature`]
15/// - [`Geometry`]
16/// - [`PointGeometry`]
17/// - [`MultiPointGeometry`]
18/// - [`LineStringGeometry`]
19/// - [`MultiLineStringGeometry`]
20/// - [`MultiPolygonGeometry`]
21/// - [`Point3DGeometry`]
22/// - [`MultiPoint3DGeometry`]
23/// - [`LineString3DGeometry`]
24/// - [`MultiLineString3DGeometry`]
25/// - [`MultiPolygon3DGeometry`]
26/// - [`VectorFeature`]
27/// - [`VectorGeometry`]
28/// - [`VectorPointGeometry`]
29/// - [`VectorMultiPointGeometry`]
30/// - [`VectorLineStringGeometry`]
31/// - [`VectorMultiLineStringGeometry`]
32/// - [`VectorMultiPolygonGeometry`]
33/// - [`VectorMultiPoint`]
34/// - [`VectorLineString`]
35/// - [`VectorMultiLineString`]
36/// - [`VectorMultiPolygon`]
37/// - `&[P]` where P implements [`GetXY`] and [`GetZ`]
38///
39/// And all specific geometries of the above enums
40pub trait LengthOfLines {
41    /// Get the total euclidean distance of a line or lines
42    fn line_length(&self) -> f64;
43}
44
45// Feature and below
46
47impl<P: GetXY + GetZ> LengthOfLines for &[P] {
48    fn line_length(&self) -> f64 {
49        let mut res = 0.;
50        let mut prev: Option<&P> = None;
51        for p in self.iter() {
52            if let Some(prev) = prev {
53                res += euclidean_distance(p, prev);
54            }
55            prev = Some(p);
56        }
57        res
58    }
59}
60
61impl<M, P: Clone + Default, D: Clone + Default> LengthOfLines for Feature<M, P, D> {
62    fn line_length(&self) -> f64 {
63        self.geometry.line_length()
64    }
65}
66impl<M: Clone + Default> LengthOfLines for Geometry<M> {
67    fn line_length(&self) -> f64 {
68        match self {
69            Geometry::Point(g) => g.line_length(),
70            Geometry::MultiPoint(g) => g.line_length(),
71            Geometry::LineString(g) => g.line_length(),
72            Geometry::MultiLineString(g) => g.line_length(),
73            Geometry::Polygon(g) => g.line_length(),
74            Geometry::MultiPolygon(g) => g.line_length(),
75            Geometry::Point3D(g) => g.line_length(),
76            Geometry::MultiPoint3D(g) => g.line_length(),
77            Geometry::LineString3D(g) => g.line_length(),
78            Geometry::MultiLineString3D(g) => g.line_length(),
79            Geometry::Polygon3D(g) => g.line_length(),
80            Geometry::MultiPolygon3D(g) => g.line_length(),
81        }
82    }
83}
84impl<M: Clone + Default> LengthOfLines for PointGeometry<M> {
85    fn line_length(&self) -> f64 {
86        self.coordinates.line_length()
87    }
88}
89impl<M: Clone + Default> LengthOfLines for MultiPointGeometry<M> {
90    fn line_length(&self) -> f64 {
91        self.coordinates.line_length()
92    }
93}
94impl<M: Clone + Default> LengthOfLines for MultiLineStringGeometry<M> {
95    fn line_length(&self) -> f64 {
96        self.coordinates.line_length()
97    }
98}
99impl<M: Clone + Default> LengthOfLines for MultiPolygonGeometry<M> {
100    fn line_length(&self) -> f64 {
101        self.coordinates.line_length()
102    }
103}
104impl<M: Clone + Default> LengthOfLines for Point3DGeometry<M> {
105    fn line_length(&self) -> f64 {
106        self.coordinates.line_length()
107    }
108}
109impl<M: Clone + Default> LengthOfLines for MultiPoint3DGeometry<M> {
110    fn line_length(&self) -> f64 {
111        self.coordinates.line_length()
112    }
113}
114impl<M: Clone + Default> LengthOfLines for MultiLineString3DGeometry<M> {
115    fn line_length(&self) -> f64 {
116        self.coordinates.line_length()
117    }
118}
119impl<M: Clone + Default> LengthOfLines for MultiPolygon3DGeometry<M> {
120    fn line_length(&self) -> f64 {
121        self.coordinates.line_length()
122    }
123}
124
125// Feature Point types
126
127impl LengthOfLines for Point {
128    fn line_length(&self) -> f64 {
129        0.
130    }
131}
132impl LengthOfLines for MultiPoint {
133    fn line_length(&self) -> f64 {
134        let mut res = 0.;
135        let mut prev: Option<&Point> = None;
136        for p in self {
137            if let Some(prev) = prev {
138                res += euclidean_distance(p, prev);
139            }
140            prev = Some(p);
141        }
142        res
143    }
144}
145impl LengthOfLines for MultiLineString {
146    fn line_length(&self) -> f64 {
147        let mut res = 0.;
148        for p in self {
149            res += p.line_length();
150        }
151        res
152    }
153}
154impl LengthOfLines for MultiPolygon {
155    fn line_length(&self) -> f64 {
156        let mut res = 0.;
157        for p in self {
158            res += p.line_length();
159        }
160        res
161    }
162}
163impl LengthOfLines for Point3D {
164    fn line_length(&self) -> f64 {
165        0.
166    }
167}
168impl LengthOfLines for MultiPoint3D {
169    fn line_length(&self) -> f64 {
170        let mut res = 0.;
171        let mut prev: Option<&Point3D> = None;
172        for p in self {
173            if let Some(prev) = prev {
174                res += euclidean_distance(p, prev);
175            }
176            prev = Some(p);
177        }
178        res
179    }
180}
181impl LengthOfLines for MultiLineString3D {
182    fn line_length(&self) -> f64 {
183        let mut res = 0.;
184        for p in self {
185            res += p.line_length();
186        }
187        res
188    }
189}
190impl LengthOfLines for MultiPolygon3D {
191    fn line_length(&self) -> f64 {
192        let mut res = 0.;
193        for p in self {
194            res += p.line_length();
195        }
196        res
197    }
198}
199
200// Vector Feature and below
201
202impl<M, P: Clone + Default, D: Clone + Default> LengthOfLines for VectorFeature<M, P, D> {
203    fn line_length(&self) -> f64 {
204        self.geometry.line_length()
205    }
206}
207impl<M: Clone + Default> LengthOfLines for VectorGeometry<M> {
208    fn line_length(&self) -> f64 {
209        match self {
210            VectorGeometry::Point(g) => g.line_length(),
211            VectorGeometry::MultiPoint(g) => g.line_length(),
212            VectorGeometry::LineString(g) => g.line_length(),
213            VectorGeometry::MultiLineString(g) => g.line_length(),
214            VectorGeometry::Polygon(g) => g.line_length(),
215            VectorGeometry::MultiPolygon(g) => g.line_length(),
216        }
217    }
218}
219impl<M: Clone + Default> LengthOfLines for VectorPointGeometry<M> {
220    fn line_length(&self) -> f64 {
221        self.coordinates.line_length()
222    }
223}
224impl<M: Clone + Default> LengthOfLines for VectorMultiPointGeometry<M> {
225    fn line_length(&self) -> f64 {
226        self.coordinates.line_length()
227    }
228}
229impl<M: Clone + Default> LengthOfLines for VectorMultiLineStringGeometry<M> {
230    fn line_length(&self) -> f64 {
231        self.coordinates.line_length()
232    }
233}
234impl<M: Clone + Default> LengthOfLines for VectorMultiPolygonGeometry<M> {
235    fn line_length(&self) -> f64 {
236        self.coordinates.line_length()
237    }
238}
239
240// Vector Point Types
241
242impl<M: Clone + Default> LengthOfLines for VectorPoint<M> {
243    fn line_length(&self) -> f64 {
244        0.
245    }
246}
247impl<M: Clone + Default> LengthOfLines for VectorMultiPoint<M> {
248    fn line_length(&self) -> f64 {
249        let mut res = 0.;
250        let mut prev: Option<&VectorPoint<M>> = None;
251        for p in self {
252            if let Some(prev) = prev {
253                res += prev.distance(p);
254            }
255            prev = Some(p);
256        }
257        res
258    }
259}
260impl<M: Clone + Default> LengthOfLines for VectorMultiLineString<M> {
261    fn line_length(&self) -> f64 {
262        let mut res = 0.;
263        for p in self {
264            res += p.line_length();
265        }
266        res
267    }
268}
269impl<M: Clone + Default> LengthOfLines for VectorMultiPolygon<M> {
270    fn line_length(&self) -> f64 {
271        let mut res = 0.;
272        for p in self {
273            res += p.line_length();
274        }
275        res
276    }
277}
278
279/// Assuming two points are on the surface of the earth as Lon-Lat coordinates, find the distance
280/// between them using the haversine formula.  Returns the distance in degrees.
281///
282/// # Parameters
283/// Both points require the [`GetXY`] trait to be implemented
284/// - `a`: The first point
285/// - `b`: The second point
286///
287/// # Returns
288/// - `f64`: The distance between the two points
289pub fn haversine_distance<P: GetXY, Q: GetXY>(a: &P, b: &Q) -> f64 {
290    let lat1 = a.y().to_radians();
291    let lat2 = b.y().to_radians();
292    let lon1 = a.x().to_radians();
293    let lon2 = b.x().to_radians();
294    let dlat = sin(0.5 * (lat2 - lat1));
295    let dlon = sin(0.5 * (lon2 - lon1));
296    let x = dlat * dlat + dlon * dlon * cos(lat1) * cos(lat2);
297    (2. * asin(sqrt(fmin(1., x)))).to_degrees()
298}
299
300/// Find the euclidean distance between two points
301///
302/// # Parameters
303/// Both points require the [`GetXY`] trait to be implemented
304/// - `a`: The first point
305/// - `b`: The second point
306///
307/// # Returns
308/// - `f64`: The distance between the two points
309pub fn euclidean_distance<P: GetXY + GetZ, Q: GetXY + GetZ>(a: &P, b: &Q) -> f64 {
310    if let Some(z1) = a.z()
311        && let Some(z2) = b.z()
312    {
313        return sqrt(pow(b.x() - a.x(), 2.) + pow(b.y() - a.y(), 2.) + pow(z2 - z1, 2.));
314    }
315    sqrt(pow(b.x() - a.x(), 2.) + pow(b.y() - a.y(), 2.))
316}