geographiclib_rs/
polygon_area.rs

1use crate::geomath::ang_diff;
2use crate::geomath::ang_normalize;
3use crate::Geodesic;
4
5use crate::geodesic_capability as caps;
6
7const POLYGONAREA_MASK: u64 =
8    caps::LATITUDE | caps::LONGITUDE | caps::DISTANCE | caps::AREA | caps::LONG_UNROLL;
9
10#[cfg(feature = "accurate")]
11use accurate::traits::*;
12
13/// Clockwise or CounterClockwise winding
14///
15/// The standard winding of a Simple Feature polygon is counter-clockwise. However, if the polygon is a hole, then the winding is clockwise.
16/// ESRI Shapefile polygons are opposite, with the outer-ring being clockwise and holes being counter-clockwise.
17#[derive(Debug, Default, Copy, Clone)]
18pub enum Winding {
19    Clockwise,
20    #[default]
21    CounterClockwise,
22}
23
24/// Compute the perimeter and area of a polygon on a Geodesic.
25#[derive(Debug, Clone)]
26pub struct PolygonArea<'a> {
27    geoid: &'a Geodesic,
28    winding: Winding,
29    num: usize,
30
31    #[cfg(not(feature = "accurate"))]
32    areasum: f64,
33
34    #[cfg(feature = "accurate")]
35    areasum: accurate::sum::Sum2<f64>,
36
37    #[cfg(not(feature = "accurate"))]
38    perimetersum: f64,
39
40    #[cfg(feature = "accurate")]
41    perimetersum: accurate::sum::Sum2<f64>,
42
43    crossings: i64,
44    initial_lat: f64,
45    initial_lon: f64,
46    latest_lat: f64,
47    latest_lon: f64,
48}
49
50/// PolygonArea can be used to compute the perimeter and area of a polygon on a Geodesic.
51///
52/// # Example
53/// ```rust
54/// use geographiclib_rs::{Geodesic, PolygonArea, Winding};
55///
56/// let g = Geodesic::wgs84();
57/// let mut pa = PolygonArea::new(&g, Winding::CounterClockwise);
58///
59/// pa.add_point(0.0, 0.0);
60/// pa.add_point(0.0, 1.0);
61/// pa.add_point(1.0, 1.0);
62/// pa.add_point(1.0, 0.0);
63///
64/// let (perimeter, area, _num) = pa.compute(false);
65///
66/// use approx::assert_relative_eq;
67/// assert_relative_eq!(perimeter, 443770.917248302);
68/// assert_relative_eq!(area, 12308778361.469452);
69/// ```
70impl<'a> PolygonArea<'a> {
71    /// Create a new PolygonArea using a Geodesic.
72    pub fn new(geoid: &'a Geodesic, winding: Winding) -> Self {
73        PolygonArea {
74            geoid,
75            winding,
76
77            num: 0,
78
79            #[cfg(not(feature = "accurate"))]
80            areasum: 0.0,
81
82            #[cfg(feature = "accurate")]
83            areasum: accurate::sum::Sum2::zero(),
84
85            #[cfg(not(feature = "accurate"))]
86            perimetersum: 0.0,
87
88            #[cfg(feature = "accurate")]
89            perimetersum: accurate::sum::Sum2::zero(),
90
91            crossings: 0,
92            initial_lat: 0.0,
93            initial_lon: 0.0,
94            latest_lat: 0.0,
95            latest_lon: 0.0,
96        }
97    }
98
99    /// Add a point to the polygon
100    pub fn add_point(&mut self, lat: f64, lon: f64) {
101        if self.num == 0 {
102            self.initial_lat = lat;
103            self.initial_lon = lon;
104        } else {
105            #[allow(non_snake_case)]
106            let (_a12, s12, _salp1, _calp1, _salp2, _calp2, _m12, _M12, _M21, S12) = self
107                .geoid
108                ._gen_inverse(self.latest_lat, self.latest_lon, lat, lon, POLYGONAREA_MASK);
109            self.perimetersum += s12;
110            self.areasum += S12;
111            self.crossings += PolygonArea::transit(self.latest_lon, lon);
112        }
113        self.latest_lat = lat;
114        self.latest_lon = lon;
115        self.num += 1;
116    }
117
118    /// Add an edge to the polygon using an azimuth (in degrees) and a distance (in meters). This can only be called after at least one point has been added.
119    ///
120    /// # Panics
121    /// Panics if no points have been added yet.
122    pub fn add_edge(&mut self, azimuth: f64, distance: f64) {
123        if self.num == 0 {
124            panic!("PolygonArea::add_edge: No points added yet");
125        }
126
127        #[allow(non_snake_case)]
128        let (_a12, lat, lon, _azi2, _s12, _m12, _M12, _M21, S12) = self.geoid._gen_direct(
129            self.latest_lat,
130            self.latest_lon,
131            azimuth,
132            false,
133            distance,
134            POLYGONAREA_MASK,
135        );
136        self.perimetersum += distance;
137        self.areasum += S12;
138        self.crossings += PolygonArea::transitdirect(self.latest_lon, lon);
139        self.latest_lat = lat;
140        self.latest_lon = lon;
141        self.num += 1;
142    }
143
144    /// Consumes the PolygonArea and returns the following tuple:
145    ///  - 0: Perimeter in (meters) of the polygon.
146    ///  - 1: Area (meters²) of the polygon.
147    ///  - 2: Number of points added to the polygon.
148    ///
149    /// # Parameters
150    ///
151    /// - `sign`: Whether to allow negative values for the area.
152    ///   - `true`: Interpret an inversely wound polygon to be a "negative" area. This will produce incorrect results if the polygon covers over half the geodesic. See "Interpreting negative area values" below.
153    ///   - `false`: Always return a positive area. Inversely wound polygons are treated as if they are always wound the same way as the winding specified during creation of the PolygonArea. This is useful if you are dealing with very large polygons that might cover over half the geodesic, or if you are certain of your winding.
154    ///
155    /// # Interpreting negative area values
156    ///
157    /// A negative value can mean one of two things:
158    /// 1. The winding of the polygon is opposite the winding specified during creation of the PolygonArea.
159    /// 2. The polygon is larger than half the planet. In this case, to get the final area of the polygon, add the area of the planet to the result. If you expect to be dealing with polygons of this size pass `signed = false` to `compute()` to get the correct result.
160    ///
161    /// # Large polygon example
162    /// ```rust
163    /// use geographiclib_rs::{Geodesic, PolygonArea, Winding};
164    /// let g = Geodesic::wgs84();
165    ///
166    /// // Describe a polygon that covers all of the earth EXCEPT this small square.
167    /// // The outside of the polygon is in this square, the inside of the polygon is the rest of the earth.
168    /// let mut pa = PolygonArea::new(&g, Winding::CounterClockwise);
169    /// pa.add_point(0.0, 0.0);
170    /// pa.add_point(1.0, 0.0);
171    /// pa.add_point(1.0, 1.0);
172    /// pa.add_point(0.0, 1.0);
173    ///
174    /// let (_perimeter, area, _count) = pa.compute(false);
175    ///
176    /// // Over 5 trillion square meters!
177    /// assert_eq!(area, 510053312945726.94);
178    /// ```
179    pub fn compute(mut self, sign: bool) -> (f64, f64, usize) {
180        #[allow(non_snake_case)]
181        let (_a12, s12, _salp1, _calp1, _salp2, _calp2, _m12, _M12, _M21, S12) =
182            self.geoid._gen_inverse(
183                self.latest_lat,
184                self.latest_lon,
185                self.initial_lat,
186                self.initial_lon,
187                POLYGONAREA_MASK,
188            );
189        self.perimetersum += s12;
190        self.areasum += S12;
191
192        let perimetersum;
193        let areasum;
194
195        #[cfg(not(feature = "accurate"))]
196        {
197            perimetersum = self.perimetersum;
198            areasum = self.areasum;
199        }
200        #[cfg(feature = "accurate")]
201        {
202            perimetersum = self.perimetersum.sum();
203            areasum = self.areasum.sum();
204        }
205
206        self.crossings += PolygonArea::transit(self.latest_lon, self.initial_lon);
207
208        // Properly take into account crossings when calculating area.
209        let areasum = self.reduce_area(areasum, sign);
210
211        (perimetersum, areasum, self.num)
212    }
213
214    /// Check what the perimeter and area would be if this point was added to the polygon without actually adding it
215    pub fn test_point(&self, lat: f64, lon: f64, sign: bool) -> (f64, f64, usize) {
216        let mut pa = self.clone();
217        pa.add_point(lat, lon);
218        pa.compute(sign)
219    }
220
221    /// Check what the perimeter and area would be if this edge was added to the polygon without actually adding it
222    pub fn test_edge(&self, azimuth: f64, distance: f64, sign: bool) -> (f64, f64, usize) {
223        let mut pa = self.clone();
224        pa.add_edge(azimuth, distance);
225        pa.compute(sign)
226    }
227
228    // Return 1 or -1 if crossing prime meridian in east or west direction.
229    // Otherwise return zero.  longitude = +/-0 considered to be positive.
230    fn transit(lon1: f64, lon2: f64) -> i64 {
231        let (lon12, _lon12s) = ang_diff(lon1, lon2);
232        let lon1 = ang_normalize(lon1);
233        let lon2 = ang_normalize(lon2);
234
235        // Translation from the following cpp code:
236        //  https://github.com/geographiclib/geographiclib/blob/8bc13eb53acdd8bc4fbe4de212d42dbb29779476/src/PolygonArea.cpp#L22
237        //
238        //  lon12 > 0 && ((lon1 < 0 && lon2 >= 0) ||
239        //  (lon1 > 0 && lon2 == 0)) ? 1 :
240        //  (lon12 < 0 && lon1 >= 0 && lon2 < 0 ? -1 : 0);
241
242        if lon12 > 0.0 && ((lon1 < 0.0 && lon2 >= 0.0) || (lon1 > 0.0 && lon2 == 0.0)) {
243            1
244        } else if lon12 < 0.0 && lon1 >= 0.0 && lon2 < 0.0 {
245            -1
246        } else {
247            0
248        }
249    }
250
251    fn transitdirect(lon1: f64, lon2: f64) -> i64 {
252        // We want to compute exactly: floor(lon2 / 360) - floor(lon1 / 360)
253        let lon1 = lon1 % 720.0;
254        let lon2 = lon2 % 720.0;
255
256        let a = if 0.0 <= lon2 && lon2 < 360.0 { 0 } else { 1 };
257
258        let b = if 0.0 <= lon1 && lon1 < 360.0 { 0 } else { 1 };
259
260        a - b
261    }
262
263    fn reduce_area(&self, area: f64, signed: bool) -> f64 {
264        let geoid_area = self.geoid.area(); // Area of the planet
265        let mut area = area % geoid_area;
266
267        // Translation of the following cpp code:
268        // if (crossings & 1) area += (area < 0 ? 1 : -1) * _area0/2;
269        if self.crossings % 2 != 0 {
270            if area < 0.0 {
271                area += geoid_area / 2.0;
272            } else {
273                area -= geoid_area / 2.0;
274            }
275        }
276
277        // Area is with the clockwise sense. If needed convert to counter-clockwise convention.
278        area = match self.winding {
279            Winding::Clockwise => area,
280            Winding::CounterClockwise => -area,
281        };
282
283        if signed {
284            // Put area in (-geoid_area/2, geoid_area/2]
285            if area > geoid_area / 2.0 {
286                area -= geoid_area;
287            } else if area <= -geoid_area / 2.0 {
288                area += geoid_area;
289            }
290        } else {
291            // Negative values mean the polygon is larger than half the planet. Correct for this.
292            if area < 0.0 {
293                area += geoid_area;
294            }
295        }
296
297        area
298    }
299}
300
301#[cfg(test)]
302mod tests {
303    use super::*;
304    use crate::Geodesic;
305    use approx::assert_relative_eq;
306
307    #[test]
308    fn test_simple_polygonarea() {
309        let geoid = Geodesic::wgs84();
310        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
311
312        pa.add_point(0.0, 0.0);
313        pa.add_point(0.0, 1.0);
314        pa.add_point(1.0, 1.0);
315        pa.add_point(1.0, 0.0);
316
317        let (perimeter, area, _count) = pa.compute(true);
318
319        assert_relative_eq!(perimeter, 443770.917, epsilon = 1.0e-3);
320        assert_relative_eq!(area, 12308778361.469, epsilon = 1.0e-3);
321
322        let mut pa = PolygonArea::new(&geoid, Winding::Clockwise);
323
324        pa.add_point(0.0, 0.0);
325        pa.add_point(0.0, 1.0);
326        pa.add_point(1.0, 1.0);
327        pa.add_point(1.0, 0.0);
328
329        let (perimeter, area, _count) = pa.compute(true);
330
331        assert_relative_eq!(perimeter, 443770.917, epsilon = 1.0e-3);
332        assert_relative_eq!(area, -12308778361.469, epsilon = 1.0e-3);
333    }
334
335    #[test]
336    fn test_planimeter0() {
337        // Copied from https://github.com/geographiclib/geographiclib-octave/blob/0662e05a432a040a60ab27c779fa09b554177ba9/inst/geographiclib_test.m#L644
338
339        let geoid = Geodesic::wgs84();
340
341        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
342        pa.add_point(89.0, 0.0);
343        pa.add_point(89.0, 90.0);
344        pa.add_point(89.0, 180.0);
345        pa.add_point(89.0, 270.0);
346        let (perimeter, area, _count) = pa.compute(true);
347        assert_relative_eq!(perimeter, 631819.8745, epsilon = 1.0e-4);
348        assert_relative_eq!(area, 24952305678.0, epsilon = 1.0);
349
350        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
351        pa.add_point(-89.0, 0.0);
352        pa.add_point(-89.0, 90.0);
353        pa.add_point(-89.0, 180.0);
354        pa.add_point(-89.0, 270.0);
355        let (perimeter, area, _count) = pa.compute(true);
356        assert_relative_eq!(perimeter, 631819.8745, epsilon = 1.0e-4);
357        assert_relative_eq!(area, -24952305678.0, epsilon = 1.0);
358
359        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
360        pa.add_point(0.0, -1.0);
361        pa.add_point(-1.0, 0.0);
362        pa.add_point(0.0, 1.0);
363        pa.add_point(1.0, 0.0);
364        let (perimeter, area, _count) = pa.compute(true);
365        assert_relative_eq!(perimeter, 627598.2731, epsilon = 1.0e-4);
366        assert_relative_eq!(area, 24619419146.0, epsilon = 1.0);
367
368        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
369        pa.add_point(90.0, 0.0);
370        pa.add_point(0.0, 0.0);
371        pa.add_point(0.0, 90.0);
372        let (perimeter, area, _count) = pa.compute(true);
373        assert_relative_eq!(perimeter, 30022685.0, epsilon = 1.0);
374        assert_relative_eq!(area, 63758202715511.0, epsilon = 1.0);
375    }
376
377    #[test]
378    fn test_planimeter5() {
379        // Copied from https://github.com/geographiclib/geographiclib-octave/blob/0662e05a432a040a60ab27c779fa09b554177ba9/inst/geographiclib_test.m#L670
380
381        let geoid = Geodesic::wgs84();
382        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
383        pa.add_point(89.0, 0.1);
384        pa.add_point(89.0, 90.1);
385        pa.add_point(89.0, -179.9);
386        let (perimeter, area, _) = pa.compute(true);
387        assert_relative_eq!(perimeter, 539297.0, epsilon = 1.0);
388        assert_relative_eq!(area, 12476152838.5, epsilon = 1.0);
389    }
390
391    #[test]
392    fn test_planimeter6() {
393        // Copied from https://github.com/geographiclib/geographiclib-octave/blob/0662e05a432a040a60ab27c779fa09b554177ba9/inst/geographiclib_test.m#L679
394
395        let geoid = Geodesic::wgs84();
396
397        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
398        pa.add_point(9.0, -0.00000000000001);
399        pa.add_point(9.0, 180.0);
400        pa.add_point(9.0, 0.0);
401        let (perimeter, area, _) = pa.compute(true);
402        assert_relative_eq!(perimeter, 36026861.0, epsilon = 1.0);
403        assert_relative_eq!(area, 0.0, epsilon = 1.0);
404
405        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
406        pa.add_point(9.0, 0.00000000000001);
407        pa.add_point(9.0, 0.0);
408        pa.add_point(9.0, 180.0);
409        let (perimeter, area, _) = pa.compute(true);
410        assert_relative_eq!(perimeter, 36026861.0, epsilon = 1.0);
411        assert_relative_eq!(area, 0.0, epsilon = 1.0);
412
413        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
414        pa.add_point(9.0, 0.00000000000001);
415        pa.add_point(9.0, 180.0);
416        pa.add_point(9.0, 0.0);
417        let (perimeter, area, _) = pa.compute(true);
418        assert_relative_eq!(perimeter, 36026861.0, epsilon = 1.0);
419        assert_relative_eq!(area, 0.0, epsilon = 1.0);
420
421        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
422        pa.add_point(9.0, -0.00000000000001);
423        pa.add_point(9.0, 0.0);
424        pa.add_point(9.0, 180.0);
425        let (perimeter, area, _) = pa.compute(true);
426        assert_relative_eq!(perimeter, 36026861.0, epsilon = 1.0);
427        assert_relative_eq!(area, 0.0, epsilon = 1.0);
428    }
429
430    #[test]
431    fn test_planimeter12() {
432        // Copied from https://github.com/geographiclib/geographiclib-octave/blob/0662e05a432a040a60ab27c779fa09b554177ba9/inst/geographiclib_test.m#L701
433        let geoid = Geodesic::wgs84();
434
435        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
436        pa.add_point(66.562222222, 0.0);
437        pa.add_point(66.562222222, 180.0);
438        pa.add_point(66.562222222, 360.0);
439        let (perimeter, area, _) = pa.compute(true);
440        assert_relative_eq!(perimeter, 10465729.0, epsilon = 1.0);
441        assert_relative_eq!(area, 0.0);
442    }
443
444    #[test]
445    fn test_planimeter12r() {
446        // Copied from https://github.com/geographiclib/geographiclib-octave/blob/0662e05a432a040a60ab27c779fa09b554177ba9/inst/geographiclib_test.m#L710
447        let geoid = Geodesic::wgs84();
448
449        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
450
451        pa.add_point(66.562222222, -0.0);
452        pa.add_point(66.562222222, -180.0);
453        pa.add_point(66.562222222, -360.0);
454
455        let (perimeter, area, _) = pa.compute(true);
456        assert_relative_eq!(perimeter, 10465729.0, epsilon = 1.0);
457        assert_relative_eq!(area, 0.0);
458    }
459
460    #[test]
461    fn test_planimeter13() {
462        // Copied from https://github.com/geographiclib/geographiclib-octave/blob/0662e05a432a040a60ab27c779fa09b554177ba9/inst/geographiclib_test.m#L719
463
464        let geoid = Geodesic::wgs84();
465
466        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
467        pa.add_point(89.0, -360.0);
468        pa.add_point(89.0, -240.0);
469        pa.add_point(89.0, -120.0);
470        pa.add_point(89.0, 0.0);
471        pa.add_point(89.0, 120.0);
472        pa.add_point(89.0, 240.0);
473        let (perimeter, area, _) = pa.compute(true);
474        assert_relative_eq!(perimeter, 1160741.0, epsilon = 1.0);
475        assert_relative_eq!(area, 32415230256.0, epsilon = 1.0);
476    }
477
478    #[test]
479    fn test_planimeter15() {
480        // Copied from https://github.com/geographiclib/geographiclib-octave/blob/0662e05a432a040a60ab27c779fa09b554177ba9/inst/geographiclib_test.m#LL728C14-L728C26
481
482        let geoid = Geodesic::wgs84();
483
484        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
485        pa.add_point(2.0, 1.0);
486        pa.add_point(1.0, 2.0);
487        pa.add_point(3.0, 3.0);
488        let (_, area, _) = pa.compute(true);
489        assert_relative_eq!(area, 18454562325.45119, epsilon = 1.0e-4);
490
491        // Switching the winding
492        let mut pa = PolygonArea::new(&geoid, Winding::Clockwise);
493        pa.add_point(2.0, 1.0);
494        pa.add_point(1.0, 2.0);
495        pa.add_point(3.0, 3.0);
496        let (_, area, _) = pa.compute(true);
497        assert_relative_eq!(area, -18454562325.45119, epsilon = 1.0e-4);
498
499        // Swaping lat and lon
500        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
501        pa.add_point(1.0, 2.0);
502        pa.add_point(2.0, 1.0);
503        pa.add_point(3.0, 3.0);
504        let (_, area, _) = pa.compute(true);
505        assert_relative_eq!(area, -18454562325.45119, epsilon = 1.0e-4);
506    }
507
508    #[test]
509    fn test_planimeter19() {
510        // Test degenerate polygons
511
512        // Copied from https://github.com/geographiclib/geographiclib-js/blob/57137fdcf4ba56718c64b909b00331754b6efceb/geodesic/test/geodesictest.js#L801
513        // Testing degenrate polygons.
514        let geoid = Geodesic::wgs84();
515
516        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
517        let (perimeter, area, _) = pa.clone().compute(true);
518        assert_relative_eq!(perimeter, 0.0);
519        assert_relative_eq!(area, 0.0);
520
521        let (perimeter, area, _) = pa.test_point(1.0, 1.0, true);
522        assert_relative_eq!(perimeter, 0.0);
523        assert_relative_eq!(area, 0.0);
524
525        let result = std::panic::catch_unwind(|| {
526            let (_, _, _) = pa.test_edge(90.0, 1000.0, true);
527        });
528        assert!(result.is_err());
529
530        pa.add_point(1.0, 1.0);
531        let (perimeter, area, _) = pa.clone().compute(true);
532        assert_relative_eq!(perimeter, 0.0);
533        assert_relative_eq!(area, 0.0);
534    }
535
536    #[test]
537    fn test_planimeter21() {
538        // Copied From: https://github.com/geographiclib/geographiclib-js/blob/57137fdcf4ba56718c64b909b00331754b6efceb/geodesic/test/geodesictest.js#LL836C13-L836C13
539
540        let g = Geodesic::wgs84();
541
542        let lat = 45.0;
543        let azi = 39.2144607176828184218;
544        let s = 8420705.40957178156285;
545        let r = 39433884866571.4277; // Area for one circuit
546        let a0 = g.area(); // Ellipsoid area
547
548        let mut pa_clockwise = PolygonArea::new(&g, Winding::Clockwise);
549        let mut pa_counter = PolygonArea::new(&g, Winding::CounterClockwise);
550
551        pa_clockwise.add_point(lat, 60.0);
552        pa_clockwise.add_point(lat, 180.0);
553        pa_clockwise.add_point(lat, -60.0);
554        pa_clockwise.add_point(lat, 60.0);
555        pa_clockwise.add_point(lat, 180.0);
556        pa_clockwise.add_point(lat, -60.0);
557
558        pa_counter.add_point(lat, 60.0);
559        pa_counter.add_point(lat, 180.0);
560        pa_counter.add_point(lat, -60.0);
561        pa_counter.add_point(lat, 60.0);
562        pa_counter.add_point(lat, 180.0);
563        pa_counter.add_point(lat, -60.0);
564
565        for i in 3..=4 {
566            pa_clockwise.add_point(lat, 60.0);
567            pa_clockwise.add_point(lat, 180.0);
568
569            pa_counter.add_point(lat, 60.0);
570            pa_counter.add_point(lat, 180.0);
571
572            let (_, area, _) = pa_counter.test_point(lat, -60.0, true);
573            assert_relative_eq!(area, i as f64 * r, epsilon = 0.5);
574
575            let (_, area, _) = pa_counter.test_point(lat, -60.0, false);
576            assert_relative_eq!(area, i as f64 * r, epsilon = 0.5);
577
578            let (_, area, _) = pa_clockwise.test_point(lat, -60.0, true);
579            assert_relative_eq!(area, -i as f64 * r, epsilon = 0.5);
580
581            let (_, area, _) = pa_clockwise.test_point(lat, -60.0, false);
582            assert_relative_eq!(area, (-i as f64 * r) + a0, epsilon = 0.5);
583
584            let (_, area, _) = pa_counter.test_edge(azi, s, true);
585            assert_relative_eq!(area, i as f64 * r, epsilon = 0.5);
586
587            let (_, area, _) = pa_counter.test_edge(azi, s, false);
588            assert_relative_eq!(area, i as f64 * r, epsilon = 0.5);
589
590            let (_, area, _) = pa_clockwise.test_edge(azi, s, true);
591            assert_relative_eq!(area, -i as f64 * r, epsilon = 0.5);
592
593            let (_, area, _) = pa_clockwise.test_edge(azi, s, false);
594            assert_relative_eq!(area, (-i as f64 * r) + a0, epsilon = 0.5);
595
596            pa_clockwise.add_point(lat, -60.0);
597            pa_counter.add_point(lat, -60.0);
598
599            let (_, area, _) = pa_counter.clone().compute(true);
600            assert_relative_eq!(area, r * i as f64, epsilon = 0.5);
601
602            let (_, area, _) = pa_counter.clone().compute(false);
603            assert_relative_eq!(area, r * i as f64, epsilon = 0.5);
604
605            let (_, area, _) = pa_clockwise.clone().compute(true);
606            assert_relative_eq!(area, -(r * i as f64), epsilon = 0.5);
607
608            let (_, area, _) = pa_clockwise.clone().compute(false);
609            assert_relative_eq!(area, -(r * i as f64) + a0, epsilon = 0.5);
610        }
611    }
612
613    #[test]
614    fn test_planimeter29() {
615        // Check transitdirect vs transit zero handling consistency
616        // Copied from: https://github.com/geographiclib/geographiclib-js/blob/57137fdcf4ba56718c64b909b00331754b6efceb/geodesic/test/geodesictest.js#L883
617
618        let geoid = Geodesic::wgs84();
619        let mut pa = PolygonArea::new(&geoid, Winding::CounterClockwise);
620        pa.add_point(0.0, 0.0);
621        pa.add_edge(90.0, 1000.0);
622        pa.add_edge(0.0, 1000.0);
623        pa.add_edge(-90.0, 1000.0);
624        let (_, area, _) = pa.compute(true);
625        assert_relative_eq!(area, 1000000.0, epsilon = 0.01);
626    }
627}
geographiclib_rs/polygon_area.rs

geographiclib_rs/
polygon_area.rs
