geometry_strategy/spherical/
area_chamberlain_duquette.rs

1//! Chamberlain–Duquette spherical polygon area.
2//!
3//! This is an opt-in, lower-cost alternative to Boost.Geometry's default
4//! spherical-excess area strategy. Boost has no corresponding strategy; the
5//! formula comes from Chamberlain & Duquette, *Some Algorithms for Polygons on
6//! a Sphere* (JPL Publication 07-03, 2007).
7
8use alloc::vec::Vec;
9
10#[cfg(not(feature = "std"))]
11use geometry_coords::math::Float;
12use geometry_cs::{CoordinateSystem, SphericalFamily};
13use geometry_tag::SameAs;
14use geometry_trait::{Point, PointOrder, Polygon, Ring};
15
16use crate::area::AreaStrategy;
17use crate::normalise::{HasAngularUnits, lonlat_radians};
18
19/// Chamberlain–Duquette area for a spherical polygon.
20///
21/// For every ring vertex `i`, the strategy accumulates
22/// `(lon[i+1] - lon[i-1]) * sin(lat[i])`, then multiplies by `R² / 2`.
23/// The result follows this crate's signed-area convention: a ring matching its
24/// declared order is positive, while an oppositely wound ring is negative.
25#[derive(Debug, Clone, Copy)]
26pub struct ChamberlainDuquetteArea {
27    /// Sphere radius in output units.
28    pub radius: f64,
29}
30
31impl ChamberlainDuquetteArea {
32    /// Mean Earth radius in metres.
33    pub const EARTH: Self = Self {
34        radius: 6_371_008.8,
35    };
36
37    /// Unit sphere; output is a solid angle in steradians.
38    pub const UNIT: Self = Self { radius: 1.0 };
39}
40
41impl Default for ChamberlainDuquetteArea {
42    fn default() -> Self {
43        Self::EARTH
44    }
45}
46
47impl<Pg> AreaStrategy<Pg> for ChamberlainDuquetteArea
48where
49    Pg: Polygon,
50    Pg::Point: Point<Scalar = f64>,
51    <Pg::Point as Point>::Cs: HasAngularUnits + CoordinateSystem,
52    <<Pg::Point as Point>::Cs as CoordinateSystem>::Family: SameAs<SphericalFamily>,
53{
54    type Out = f64;
55
56    fn area(&self, polygon: &Pg) -> Self::Out {
57        let mut total = ring_area(polygon.exterior(), self.radius);
58        for interior in polygon.interiors() {
59            total += ring_area(interior, self.radius);
60        }
61        total
62    }
63}
64
65fn ring_area<R>(ring: &R, radius: f64) -> f64
66where
67    R: Ring,
68    R::Point: Point<Scalar = f64>,
69    <R::Point as Point>::Cs: HasAngularUnits,
70{
71    let mut coordinates: Vec<(f64, f64)> = ring.points().map(lonlat_radians).collect();
72    if coordinates.len() > 1 && coordinates.first() == coordinates.last() {
73        coordinates.pop();
74    }
75    if coordinates.len() < 3 {
76        return 0.0;
77    }
78
79    let mut sum = 0.0;
80    for index in 0..coordinates.len() {
81        let previous = coordinates[(index + coordinates.len() - 1) % coordinates.len()];
82        let current = coordinates[index];
83        let next = coordinates[(index + 1) % coordinates.len()];
84        sum += longitude_delta(next.0 - previous.0) * current.1.sin();
85    }
86    let signed = sum * radius * radius / 2.0;
87    match ring.point_order() {
88        PointOrder::Clockwise => signed,
89        PointOrder::CounterClockwise => -signed,
90    }
91}
92
93fn longitude_delta(delta: f64) -> f64 {
94    (delta + core::f64::consts::PI).rem_euclid(core::f64::consts::TAU) - core::f64::consts::PI
95}
96
97#[cfg(test)]
98mod tests {
99    use geometry_cs::{Degree, Spherical};
100    use geometry_model::{Point2D, Polygon, Ring};
101
102    use super::*;
103
104    #[test]
105    fn one_degree_square_on_unit_sphere() {
106        type P = Point2D<f64, Spherical<Degree>>;
107        let polygon: Polygon<P> = Polygon::new(Ring::from_vec(alloc::vec![
108            P::new(0.0, 0.0),
109            P::new(0.0, 1.0),
110            P::new(1.0, 1.0),
111            P::new(1.0, 0.0),
112            P::new(0.0, 0.0),
113        ]));
114        let area = ChamberlainDuquetteArea::UNIT.area(&polygon);
115        assert!((area - 0.000_304_601_954_726_850_5).abs() < 1e-12);
116    }
117}
geometry_strategy/spherical/area_chamberlain_duquette.rs

geometry_strategy/spherical/
area_chamberlain_duquette.rs
