gistools/geometry/tools/polys/
area.rs

1use crate::space::EARTH_RADIUS;
2use alloc::vec::Vec;
3use libm::sin;
4use s2json::{
5    Feature, Features, Geometry, GetXY, MultiLineString, MultiLineString3D,
6    MultiLineString3DGeometry, MultiLineStringGeometry, MultiPoint, MultiPoint3D,
7    MultiPoint3DGeometry, MultiPointGeometry, MultiPolygon, MultiPolygon3D, MultiPolygon3DGeometry,
8    MultiPolygonGeometry, Point, Point3D, Point3DGeometry, PointGeometry, VectorFeature,
9    VectorGeometry, VectorMultiLineString, VectorMultiLineStringGeometry, VectorMultiPoint,
10    VectorMultiPointGeometry, VectorMultiPolygon, VectorMultiPolygonGeometry, VectorPoint,
11    VectorPointGeometry,
12};
13
14/// Get the area of the polygon. Lines return 0 if not closed. Other geometries return 0.
15///
16/// Assumes geometry is in Lon-Lat space
17///
18/// If no radius is provided, the Earth's radius is used
19///
20/// This trait is implemented for:
21/// - [`Feature`]
22/// - [`Geometry`]
23/// - [`PointGeometry`]
24/// - [`MultiPointGeometry`]
25/// - [`s2json::LineStringGeometry`]
26/// - [`MultiLineStringGeometry`]
27/// - [`MultiPolygonGeometry`]
28/// - [`Point3DGeometry`]
29/// - [`MultiPoint3DGeometry`]
30/// - [`s2json::LineString3DGeometry`]
31/// - [`MultiLineString3DGeometry`]
32/// - [`MultiPolygon3DGeometry`]
33/// - [`VectorFeature`]
34/// - [`VectorGeometry`]
35/// - [`VectorPointGeometry`]
36/// - [`VectorMultiPointGeometry`]
37/// - [`s2json::VectorLineStringGeometry`]
38/// - [`VectorMultiLineStringGeometry`]
39/// - [`VectorMultiPolygonGeometry`]
40/// - [`VectorMultiPoint`]
41/// - [`VectorMultiLineString`]
42/// - [`VectorMultiPolygon`]
43/// - [`Features`]
44/// - `&Vec<P>` or `&[P]` where P implements [`GetXY`]
45///
46/// And all specific geometries of the above enums
47pub trait Area {
48    /// Get the area of the polygon. Lines return 0 if not closed. Other geometries return 0.
49    ///
50    /// Assumes geometry is in Lon-Lat space
51    ///
52    /// If no radius is provided, the Earth's radius is used
53    fn area(&self, radius: Option<f64>) -> f64;
54}
55
56// Feature and below
57
58impl<P: GetXY> Area for &Vec<P> {
59    fn area(&self, radius: Option<f64>) -> f64 {
60        ring_area(self, radius.unwrap_or(EARTH_RADIUS))
61    }
62}
63impl<P: GetXY> Area for &mut Vec<P> {
64    fn area(&self, radius: Option<f64>) -> f64 {
65        ring_area(self, radius.unwrap_or(EARTH_RADIUS))
66    }
67}
68impl<P: GetXY> Area for &[P] {
69    fn area(&self, radius: Option<f64>) -> f64 {
70        ring_area(self, radius.unwrap_or(EARTH_RADIUS))
71    }
72}
73impl<P: GetXY> Area for &mut [P] {
74    fn area(&self, radius: Option<f64>) -> f64 {
75        ring_area(self, radius.unwrap_or(EARTH_RADIUS))
76    }
77}
78
79impl<M, P: Clone + Default, D: Clone + Default> Area for Feature<M, P, D> {
80    fn area(&self, radius: Option<f64>) -> f64 {
81        self.geometry.area(radius)
82    }
83}
84impl<M: Clone + Default> Area for Geometry<M> {
85    fn area(&self, radius: Option<f64>) -> f64 {
86        match self {
87            Geometry::Point(g) => g.area(radius),
88            Geometry::MultiPoint(g) => g.area(radius),
89            Geometry::LineString(g) => g.area(radius),
90            Geometry::MultiLineString(g) => g.area(radius),
91            Geometry::Polygon(g) => g.area(radius),
92            Geometry::MultiPolygon(g) => g.area(radius),
93            Geometry::Point3D(g) => g.area(radius),
94            Geometry::MultiPoint3D(g) => g.area(radius),
95            Geometry::LineString3D(g) => g.area(radius),
96            Geometry::MultiLineString3D(g) => g.area(radius),
97            Geometry::Polygon3D(g) => g.area(radius),
98            Geometry::MultiPolygon3D(g) => g.area(radius),
99        }
100    }
101}
102impl<M: Clone + Default> Area for PointGeometry<M> {
103    fn area(&self, radius: Option<f64>) -> f64 {
104        self.coordinates.area(radius)
105    }
106}
107impl<M: Clone + Default> Area for MultiPointGeometry<M> {
108    fn area(&self, radius: Option<f64>) -> f64 {
109        self.coordinates.area(radius)
110    }
111}
112impl<M: Clone + Default> Area for MultiLineStringGeometry<M> {
113    fn area(&self, radius: Option<f64>) -> f64 {
114        self.coordinates.area(radius)
115    }
116}
117impl<M: Clone + Default> Area for MultiPolygonGeometry<M> {
118    fn area(&self, radius: Option<f64>) -> f64 {
119        self.coordinates.area(radius)
120    }
121}
122impl<M: Clone + Default> Area for Point3DGeometry<M> {
123    fn area(&self, radius: Option<f64>) -> f64 {
124        self.coordinates.area(radius)
125    }
126}
127impl<M: Clone + Default> Area for MultiPoint3DGeometry<M> {
128    fn area(&self, radius: Option<f64>) -> f64 {
129        self.coordinates.area(radius)
130    }
131}
132impl<M: Clone + Default> Area for MultiLineString3DGeometry<M> {
133    fn area(&self, radius: Option<f64>) -> f64 {
134        self.coordinates.area(radius)
135    }
136}
137impl<M: Clone + Default> Area for MultiPolygon3DGeometry<M> {
138    fn area(&self, radius: Option<f64>) -> f64 {
139        self.coordinates.area(radius)
140    }
141}
142
143// Feature Point types
144
145impl Area for Point {
146    fn area(&self, _radius: Option<f64>) -> f64 {
147        0.
148    }
149}
150impl Area for MultiPoint {
151    fn area(&self, radius: Option<f64>) -> f64 {
152        if self.first() != self.last() {
153            0.
154        } else {
155            ring_area(self, radius.unwrap_or(EARTH_RADIUS))
156        }
157    }
158}
159impl Area for MultiLineString {
160    fn area(&self, radius: Option<f64>) -> f64 {
161        // first line adds, all others subtract
162        let mut total = 0.;
163        for (i, line) in self.iter().enumerate() {
164            if i == 0 {
165                total += line.area(radius);
166            } else {
167                total -= line.area(radius);
168            }
169        }
170        total
171    }
172}
173impl Area for MultiPolygon {
174    fn area(&self, radius: Option<f64>) -> f64 {
175        let mut total = 0.;
176        for poly in self {
177            total += poly.area(radius);
178        }
179        total
180    }
181}
182impl Area for Point3D {
183    fn area(&self, _radius: Option<f64>) -> f64 {
184        0.
185    }
186}
187impl Area for MultiPoint3D {
188    fn area(&self, radius: Option<f64>) -> f64 {
189        if self.first() != self.last() {
190            0.
191        } else {
192            ring_area(self, radius.unwrap_or(EARTH_RADIUS))
193        }
194    }
195}
196impl Area for MultiLineString3D {
197    fn area(&self, radius: Option<f64>) -> f64 {
198        // first line adds, all others subtract
199        let mut total = 0.;
200        for (i, line) in self.iter().enumerate() {
201            if i == 0 {
202                total += line.area(radius);
203            } else {
204                total -= line.area(radius);
205            }
206        }
207        total
208    }
209}
210impl Area for MultiPolygon3D {
211    fn area(&self, radius: Option<f64>) -> f64 {
212        let mut total = 0.;
213        for poly in self {
214            total += poly.area(radius);
215        }
216        total
217    }
218}
219
220// Vector Feature and below
221
222impl<M, P: Clone + Default, D: Clone + Default> Area for VectorFeature<M, P, D> {
223    fn area(&self, radius: Option<f64>) -> f64 {
224        self.geometry.area(radius)
225    }
226}
227impl<M: Clone + Default> Area for VectorGeometry<M> {
228    fn area(&self, radius: Option<f64>) -> f64 {
229        match self {
230            VectorGeometry::Point(g) => g.area(radius),
231            VectorGeometry::MultiPoint(g) => g.area(radius),
232            VectorGeometry::LineString(g) => g.area(radius),
233            VectorGeometry::MultiLineString(g) => g.area(radius),
234            VectorGeometry::Polygon(g) => g.area(radius),
235            VectorGeometry::MultiPolygon(g) => g.area(radius),
236        }
237    }
238}
239impl<M: Clone + Default> Area for VectorPointGeometry<M> {
240    fn area(&self, radius: Option<f64>) -> f64 {
241        self.coordinates.area(radius)
242    }
243}
244impl<M: Clone + Default> Area for VectorMultiPointGeometry<M> {
245    fn area(&self, radius: Option<f64>) -> f64 {
246        self.coordinates.area(radius)
247    }
248}
249impl<M: Clone + Default> Area for VectorMultiLineStringGeometry<M> {
250    fn area(&self, radius: Option<f64>) -> f64 {
251        self.coordinates.area(radius)
252    }
253}
254impl<M: Clone + Default> Area for VectorMultiPolygonGeometry<M> {
255    fn area(&self, radius: Option<f64>) -> f64 {
256        self.coordinates.area(radius)
257    }
258}
259
260// Vector Point Types
261
262impl<M: Clone + Default> Area for VectorPoint<M> {
263    fn area(&self, _radius: Option<f64>) -> f64 {
264        0.
265    }
266}
267impl<M: Clone + Default> Area for VectorMultiPoint<M> {
268    fn area(&self, radius: Option<f64>) -> f64 {
269        if self.first() != self.last() {
270            0.
271        } else {
272            ring_area(self, radius.unwrap_or(EARTH_RADIUS))
273        }
274    }
275}
276impl<M: Clone + Default> Area for VectorMultiLineString<M> {
277    fn area(&self, radius: Option<f64>) -> f64 {
278        // first line adds, all others subtract
279        let mut total = 0.;
280        for (i, line) in self.iter().enumerate() {
281            if i == 0 {
282                total += line.area(radius);
283            } else {
284                total -= line.area(radius);
285            }
286        }
287        total
288    }
289}
290impl<M: Clone + Default> Area for VectorMultiPolygon<M> {
291    fn area(&self, radius: Option<f64>) -> f64 {
292        let mut total = 0.;
293        for poly in self {
294            total += poly.area(radius);
295        }
296        total
297    }
298}
299
300// Features
301
302impl<M, P: Clone + Default, D: Clone + Default> Area for Features<M, P, D> {
303    fn area(&self, radius: Option<f64>) -> f64 {
304        match self {
305            Features::Feature(f) => f.area(radius),
306            Features::VectorFeature(f) => f.area(radius),
307        }
308    }
309}
310
311/// Calculate the approximate area of the polygon were it projected onto the planet.
312/// Note that this area will be positive if ring is oriented counter-clockwise,
313/// otherwise it will be negative.
314///
315/// Reference:
316/// Robert. G. Chamberlain and William H. Duquette, "Some Algorithms for Polygons on a Sphere",
317/// JPL Publication 07-03, Jet Propulsion
318/// Laboratory, Pasadena, CA, June 2007 `https://trs.jpl.nasa.gov/handle/2014/40409`
319///
320/// ## Parameters
321/// - `coords`: ring Coordinates in lon-lat space
322/// - `planetRadius`: the radius of the planet (Earth by default)
323///
324/// ## Returns
325/// The approximate signed geodesic area of the polygon in square meters.
326pub fn ring_area<P: GetXY>(coords: &[P], radius: f64) -> f64 {
327    let coords_length = coords.len() - 1;
328    let factor = (radius * radius) / 2.;
329
330    if coords_length <= 2 {
331        return 0.;
332    }
333
334    let mut total = 0.;
335    let mut i = 0;
336    while i < coords_length {
337        let lower = &coords[i];
338        let middle = &coords[if i + 1 == coords_length { 0 } else { i + 1 }];
339        let upper = &coords[if i + 2 >= coords_length { (i + 2) % coords_length } else { i + 2 }];
340
341        let lower_x = lower.x().to_radians();
342        let middle_y = middle.y().to_radians();
343        let upper_x = upper.x().to_radians();
344
345        total += (upper_x - lower_x) * sin(middle_y);
346
347        i += 1;
348    }
349
350    -(total * factor)
351}
gistools/geometry/tools/polys/area.rs

gistools/geometry/tools/polys/
area.rs
