oxigdal_algorithms/vector/
area.rs

1//! Area calculation for polygonal geometries
2//!
3//! This module provides functions for computing areas of polygons using
4//! different methods appropriate for planar and geodetic coordinates.
5//!
6//! # Area Methods
7//!
8//! - **Planar**: Fast area calculation assuming Cartesian coordinates (projected data)
9//! - **Geodetic**: Accurate area on Earth's surface using WGS84 ellipsoid
10//! - **Signed**: Preserves orientation information (counter-clockwise = positive)
11//!
12//! # Examples
13//!
14//! ```
15//! use oxigdal_algorithms::vector::{Coordinate, Polygon, LineString, area_polygon, AreaMethod};
16//!
17//! # fn main() -> Result<(), Box<dyn std::error::Error>> {
18//! let coords = vec![
19//!     Coordinate::new_2d(0.0, 0.0),
20//!     Coordinate::new_2d(10.0, 0.0),
21//!     Coordinate::new_2d(10.0, 10.0),
22//!     Coordinate::new_2d(0.0, 10.0),
23//!     Coordinate::new_2d(0.0, 0.0),
24//! ];
25//! let exterior = LineString::new(coords)?;
26//! let polygon = Polygon::new(exterior, vec![])?;
27//! let area_planar = area_polygon(&polygon, AreaMethod::Planar)?;
28//! // Planar area = 100.0 square units
29//! # Ok(())
30//! # }
31//! ```
32
33use crate::error::{AlgorithmError, Result};
34use oxigdal_core::vector::{Coordinate, Geometry, MultiPolygon, Polygon};
35
36#[cfg(feature = "std")]
37use std::f64::consts::PI;
38
39#[cfg(not(feature = "std"))]
40use core::f64::consts::PI;
41
42/// WGS84 ellipsoid semi-major axis (meters)
43const WGS84_A: f64 = 6_378_137.0;
44
45/// WGS84 ellipsoid semi-minor axis (meters)
46const WGS84_B: f64 = 6_356_752.314_245;
47
48/// WGS84 ellipsoid flattening
49const WGS84_F: f64 = (WGS84_A - WGS84_B) / WGS84_A;
50
51/// Area calculation method
52#[derive(Debug, Clone, Copy, PartialEq, Eq)]
53pub enum AreaMethod {
54    /// Planar area using Cartesian coordinates (fast, for projected data)
55    Planar,
56    /// Geodetic area on WGS84 ellipsoid (accurate for lat/lon coordinates)
57    Geodetic,
58    /// Signed planar area (preserves orientation)
59    SignedPlanar,
60}
61
62/// Computes the area of a geometry
63///
64/// # Arguments
65///
66/// * `geometry` - Input geometry
67/// * `method` - Area calculation method
68///
69/// # Returns
70///
71/// Area in appropriate units (square meters for geodetic, coordinate units for planar)
72///
73/// # Errors
74///
75/// Returns error if geometry is not polygonal or invalid
76pub fn area(geometry: &Geometry, method: AreaMethod) -> Result<f64> {
77    match geometry {
78        Geometry::Polygon(p) => area_polygon(p, method),
79        Geometry::MultiPolygon(mp) => area_multipolygon(mp, method),
80        _ => Err(AlgorithmError::GeometryError {
81            message: "area calculation requires Polygon or MultiPolygon geometry".to_string(),
82        }),
83    }
84}
85
86/// Computes the area of a polygon
87///
88/// # Arguments
89///
90/// * `polygon` - Input polygon
91/// * `method` - Area calculation method
92///
93/// # Returns
94///
95/// Area value (always non-negative unless using SignedPlanar method)
96///
97/// # Errors
98///
99/// Returns error if polygon is invalid
100pub fn area_polygon(polygon: &Polygon, method: AreaMethod) -> Result<f64> {
101    if polygon.exterior.coords.len() < 3 {
102        return Err(AlgorithmError::InsufficientData {
103            operation: "area_polygon",
104            message: "polygon must have at least 3 coordinates".to_string(),
105        });
106    }
107
108    match method {
109        AreaMethod::Planar => Ok(area_polygon_planar(polygon)),
110        AreaMethod::Geodetic => area_polygon_geodetic(polygon),
111        AreaMethod::SignedPlanar => Ok(area_polygon_signed(polygon)),
112    }
113}
114
115/// Computes the area of a multipolygon
116///
117/// # Arguments
118///
119/// * `multipolygon` - Input multipolygon
120/// * `method` - Area calculation method
121///
122/// # Returns
123///
124/// Total area of all polygons
125///
126/// # Errors
127///
128/// Returns error if any polygon is invalid
129pub fn area_multipolygon(multipolygon: &MultiPolygon, method: AreaMethod) -> Result<f64> {
130    if multipolygon.polygons.is_empty() {
131        return Ok(0.0);
132    }
133
134    let mut total = 0.0;
135    for polygon in &multipolygon.polygons {
136        total += area_polygon(polygon, method)?;
137    }
138
139    Ok(total)
140}
141
142/// Computes planar area of a polygon (fast, Cartesian)
143///
144/// Uses the shoelace formula (Gauss's area formula).
145fn area_polygon_planar(polygon: &Polygon) -> f64 {
146    let mut area = ring_area_planar(&polygon.exterior.coords).abs();
147
148    // Subtract holes
149    for hole in &polygon.interiors {
150        area -= ring_area_planar(&hole.coords).abs();
151    }
152
153    area
154}
155
156/// Computes signed planar area of a polygon
157///
158/// Positive for counter-clockwise orientation, negative for clockwise.
159fn area_polygon_signed(polygon: &Polygon) -> f64 {
160    let mut area = ring_area_planar(&polygon.exterior.coords);
161
162    // Subtract holes (with their signs)
163    for hole in &polygon.interiors {
164        area -= ring_area_planar(&hole.coords);
165    }
166
167    area
168}
169
170/// Computes planar area of a ring using shoelace formula
171fn ring_area_planar(coords: &[Coordinate]) -> f64 {
172    if coords.len() < 3 {
173        return 0.0;
174    }
175
176    let mut area = 0.0;
177    let n = coords.len();
178
179    for i in 0..n {
180        let j = (i + 1) % n;
181        area += coords[i].x * coords[j].y;
182        area -= coords[j].x * coords[i].y;
183    }
184
185    area / 2.0
186}
187
188/// Computes geodetic area of a polygon on WGS84 ellipsoid
189///
190/// Uses the method from Bessel (1825) for geodetic polygons.
191/// Coordinates must be in longitude/latitude (degrees).
192fn area_polygon_geodetic(polygon: &Polygon) -> Result<f64> {
193    let mut area = ring_area_geodetic(&polygon.exterior.coords)?;
194
195    // Subtract holes
196    for hole in &polygon.interiors {
197        area -= ring_area_geodetic(&hole.coords)?;
198    }
199
200    Ok(area.abs())
201}
202
203/// Computes geodetic area of a ring
204///
205/// Implementation based on:
206/// Chamberlain, R. G., and W. H. Duquette. "Some algorithms for polygons on a sphere."
207/// JPL Publication 07-03 (2007).
208fn ring_area_geodetic(coords: &[Coordinate]) -> Result<f64> {
209    if coords.len() < 3 {
210        return Ok(0.0);
211    }
212
213    let n = coords.len();
214    let mut area = 0.0;
215
216    for i in 0..n - 1 {
217        let p1 = &coords[i];
218        let p2 = &coords[i + 1];
219
220        // Convert to radians
221        let lat1 = p1.y * PI / 180.0;
222        let lon1 = p1.x * PI / 180.0;
223        let lat2 = p2.y * PI / 180.0;
224        let lon2 = p2.x * PI / 180.0;
225
226        // Validate latitude range
227        if lat1.abs() > PI / 2.0 || lat2.abs() > PI / 2.0 {
228            return Err(AlgorithmError::InvalidParameter {
229                parameter: "latitude",
230                message: "latitude must be between -90 and 90 degrees".to_string(),
231            });
232        }
233
234        // Calculate area contribution of this edge
235        let dlon = lon2 - lon1;
236
237        // Handle antimeridian crossing
238        let dlon = if dlon > PI {
239            dlon - 2.0 * PI
240        } else if dlon < -PI {
241            dlon + 2.0 * PI
242        } else {
243            dlon
244        };
245
246        // Accumulate area using trapezoid rule on the sphere
247        area += dlon * (lat1.sin() + lat2.sin());
248    }
249
250    // Convert to square meters using WGS84 parameters
251    let authalic_radius =
252        ((WGS84_A * WGS84_A + WGS84_B * WGS84_B / (1.0 - WGS84_F).sqrt()) / 2.0).sqrt();
253    let area_m2 = area.abs() * authalic_radius * authalic_radius / 2.0;
254
255    Ok(area_m2)
256}
257
258/// Checks if a polygon is oriented counter-clockwise
259///
260/// # Arguments
261///
262/// * `polygon` - Input polygon
263///
264/// # Returns
265///
266/// True if exterior ring is counter-clockwise, false otherwise
267pub fn is_counter_clockwise(polygon: &Polygon) -> bool {
268    ring_area_planar(&polygon.exterior.coords) > 0.0
269}
270
271/// Checks if a polygon is oriented clockwise
272///
273/// # Arguments
274///
275/// * `polygon` - Input polygon
276///
277/// # Returns
278///
279/// True if exterior ring is clockwise, false otherwise
280pub fn is_clockwise(polygon: &Polygon) -> bool {
281    ring_area_planar(&polygon.exterior.coords) < 0.0
282}
283
284#[cfg(test)]
285mod tests {
286    use super::*;
287    use oxigdal_core::vector::LineString;
288
289    fn create_square(size: f64) -> Result<Polygon> {
290        let coords = vec![
291            Coordinate::new_2d(0.0, 0.0),
292            Coordinate::new_2d(size, 0.0),
293            Coordinate::new_2d(size, size),
294            Coordinate::new_2d(0.0, size),
295            Coordinate::new_2d(0.0, 0.0),
296        ];
297        let exterior = LineString::new(coords).map_err(|e| AlgorithmError::Core(e))?;
298        Polygon::new(exterior, vec![]).map_err(|e| AlgorithmError::Core(e))
299    }
300
301    fn create_square_with_hole() -> Result<Polygon> {
302        let exterior_coords = vec![
303            Coordinate::new_2d(0.0, 0.0),
304            Coordinate::new_2d(10.0, 0.0),
305            Coordinate::new_2d(10.0, 10.0),
306            Coordinate::new_2d(0.0, 10.0),
307            Coordinate::new_2d(0.0, 0.0),
308        ];
309        let hole_coords = vec![
310            Coordinate::new_2d(2.0, 2.0),
311            Coordinate::new_2d(8.0, 2.0),
312            Coordinate::new_2d(8.0, 8.0),
313            Coordinate::new_2d(2.0, 8.0),
314            Coordinate::new_2d(2.0, 2.0),
315        ];
316
317        let exterior = LineString::new(exterior_coords).map_err(|e| AlgorithmError::Core(e))?;
318        let hole = LineString::new(hole_coords).map_err(|e| AlgorithmError::Core(e))?;
319
320        Polygon::new(exterior, vec![hole]).map_err(|e| AlgorithmError::Core(e))
321    }
322
323    #[test]
324    fn test_area_polygon_planar() {
325        let poly = create_square(10.0);
326        assert!(poly.is_ok());
327
328        if let Ok(p) = poly {
329            let result = area_polygon(&p, AreaMethod::Planar);
330            assert!(result.is_ok());
331
332            if let Ok(area) = result {
333                assert!((area - 100.0).abs() < 1e-10);
334            }
335        }
336    }
337
338    #[test]
339    fn test_area_polygon_with_hole() {
340        let poly = create_square_with_hole();
341        assert!(poly.is_ok());
342
343        if let Ok(p) = poly {
344            let result = area_polygon(&p, AreaMethod::Planar);
345            assert!(result.is_ok());
346
347            if let Ok(area) = result {
348                // Outer area = 100, hole area = 36, effective = 64
349                assert!((area - 64.0).abs() < 1e-10);
350            }
351        }
352    }
353
354    #[test]
355    fn test_area_polygon_signed() {
356        let poly = create_square(10.0);
357        assert!(poly.is_ok());
358
359        if let Ok(p) = poly {
360            let result = area_polygon(&p, AreaMethod::SignedPlanar);
361            assert!(result.is_ok());
362
363            if let Ok(area) = result {
364                // CCW orientation = positive area
365                assert!(area > 0.0);
366                assert!((area.abs() - 100.0).abs() < 1e-10);
367            }
368        }
369    }
370
371    #[test]
372    fn test_area_multipolygon() {
373        let poly1 = create_square(5.0);
374        let poly2 = create_square(3.0);
375
376        assert!(poly1.is_ok() && poly2.is_ok());
377
378        if let (Ok(p1), Ok(p2)) = (poly1, poly2) {
379            let mp = MultiPolygon::new(vec![p1, p2]);
380            let result = area_multipolygon(&mp, AreaMethod::Planar);
381            assert!(result.is_ok());
382
383            if let Ok(area) = result {
384                // 5*5 + 3*3 = 25 + 9 = 34
385                assert!((area - 34.0).abs() < 1e-10);
386            }
387        }
388    }
389
390    #[test]
391    fn test_area_geodetic_small_polygon() {
392        // Small polygon in degrees (approximately 1 degree square near equator)
393        let coords = vec![
394            Coordinate::new_2d(0.0, 0.0),
395            Coordinate::new_2d(1.0, 0.0),
396            Coordinate::new_2d(1.0, 1.0),
397            Coordinate::new_2d(0.0, 1.0),
398            Coordinate::new_2d(0.0, 0.0),
399        ];
400        let exterior = LineString::new(coords);
401        assert!(exterior.is_ok());
402
403        if let Ok(ext) = exterior {
404            let poly = Polygon::new(ext, vec![]);
405            assert!(poly.is_ok());
406
407            if let Ok(p) = poly {
408                let result = area_polygon(&p, AreaMethod::Geodetic);
409                assert!(result.is_ok());
410
411                if let Ok(area) = result {
412                    // 1 degree at equator is approximately 111 km
413                    // So 1 degree square is approximately 111 * 111 = 12,321 km^2
414                    // = 12,321,000,000 m^2
415                    assert!(area > 1e10); // Should be in billions of square meters
416                    assert!(area < 2e10); // Rough upper bound
417                }
418            }
419        }
420    }
421
422    #[test]
423    fn test_is_counter_clockwise() {
424        let poly = create_square(10.0);
425        assert!(poly.is_ok());
426
427        if let Ok(p) = poly {
428            assert!(is_counter_clockwise(&p));
429            assert!(!is_clockwise(&p));
430        }
431    }
432
433    #[test]
434    fn test_ring_area_planar() {
435        let coords = vec![
436            Coordinate::new_2d(0.0, 0.0),
437            Coordinate::new_2d(10.0, 0.0),
438            Coordinate::new_2d(10.0, 10.0),
439            Coordinate::new_2d(0.0, 10.0),
440            Coordinate::new_2d(0.0, 0.0),
441        ];
442
443        let area = ring_area_planar(&coords);
444        assert!((area.abs() - 100.0).abs() < 1e-10);
445    }
446
447    #[test]
448    fn test_area_empty_multipolygon() {
449        let mp = MultiPolygon::empty();
450        let result = area_multipolygon(&mp, AreaMethod::Planar);
451        assert!(result.is_ok());
452
453        if let Ok(area) = result {
454            assert_eq!(area, 0.0);
455        }
456    }
457
458    #[test]
459    fn test_area_invalid_geometry() {
460        let point = Geometry::Point(oxigdal_core::vector::Point::new(0.0, 0.0));
461        let result = area(&point, AreaMethod::Planar);
462        assert!(result.is_err());
463    }
464
465    #[test]
466    fn test_geodetic_area_invalid_latitude() {
467        // Latitude > 90 degrees
468        let coords = vec![
469            Coordinate::new_2d(0.0, 0.0),
470            Coordinate::new_2d(1.0, 0.0),
471            Coordinate::new_2d(1.0, 100.0), // Invalid
472            Coordinate::new_2d(0.0, 1.0),
473            Coordinate::new_2d(0.0, 0.0),
474        ];
475        let exterior = LineString::new(coords);
476        assert!(exterior.is_ok());
477
478        if let Ok(ext) = exterior {
479            let poly = Polygon::new(ext, vec![]);
480            assert!(poly.is_ok());
481
482            if let Ok(p) = poly {
483                let result = area_polygon(&p, AreaMethod::Geodetic);
484                assert!(result.is_err());
485            }
486        }
487    }
488}
oxigdal_algorithms/vector/area.rs

oxigdal_algorithms/vector/
area.rs
