fix(analysis): calculate geometry attributes on the WGS 84 ellipsoid

crates/analysis/src/tools.rs

@@ -3,7 +3,7 @@
 use geo::{BooleanOps, Contains, ConvexHull, InteriorPoint, Intersects, MultiPolygon, Orient, Simplify, orient::Direction};
 use rustgis_data::raster::Crs;
 use rustgis_data::{Coord, Dataset, Field, FieldType, Geometry, GeometryType, Value};
-use rustgis_geom::{Extent, geodesic_length_m, geodesic_ring_area_m2};
+use rustgis_geom::Extent;
 use serde_json::Value as Json;
 
 use crate::convert::{from_line, from_multipolygon, from_point, to_geo, to_multiline, to_multipolygon};
@@ -1236,6 +1236,8 @@ fn count_points(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
 }
 
 fn calculate_geometry(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
+    use geo::{Geodesic, GeodesicArea, Length};
+    let line = |coords: &[Coord]| geo::LineString::new(coords.iter().map(|c| geo::Coord { x: c.x, y: c.y }).collect());
     let ds = input(inp, p, "input")?;
     let mut out = ds.clone();
     match ds.geometry_type {
@@ -1247,9 +1249,26 @@ fn calculate_geometry(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String>
                 if let Some(Geometry::Polygon(polys)) = &f.geometry {
                     for rings in polys {
                         for (i, r) in rings.iter().enumerate() {
-                            let ra = geodesic_ring_area_m2(r);
+                            let mut ring = line(r);
+                            ring.close();
+                            // Orient in continuous longitude space, including small dateline crossings.
+                            let mut longitude = ring.0.first().map_or(0.0, |c| c.x);
+                            for c in &mut ring.0 {
+                                longitude += (c.x - longitude + 180.0).rem_euclid(360.0) - 180.0;
+                                c.x = longitude;
+                            }
+                            let polar = ring.0.first().zip(ring.0.last()).is_some_and(|(a, b)| (a.x - b.x).abs() > 180.0);
+                            let orientable = !polar && ring.0.iter().all(|c| c.x.is_finite() && (-90.0..=90.0).contains(&c.y));
+                            let polygon = geo::Polygon::new(ring, vec![]);
+                            let (rp, ra) = if orientable {
+                                polygon.orient(Direction::Default).geodesic_perimeter_area_unsigned()
+                            } else {
+                                // Keep smaller polar caps; invalid coordinates bypass planar orientation.
+                                let (perimeter, area) = polygon.geodesic_perimeter_area_signed();
+                                (perimeter, area.abs())
+                            };
                             a += if i == 0 { ra } else { -ra };
-                            per += geodesic_length_m(r);
+                            per += rp;
                         }
                     }
                 }
@@ -1265,7 +1284,7 @@ fn calculate_geometry(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String>
             let fl = out.ensure_field("length_km", FieldType::Double);
             for f in &mut out.features {
                 let len: f64 = match &f.geometry {
-                    Some(Geometry::Line(parts)) => parts.iter().map(|p| geodesic_length_m(p)).sum(),
+                    Some(Geometry::Line(parts)) => parts.iter().map(|p| Geodesic.length(&line(p))).sum(),
                     _ => 0.0,
                 };
                 if let Some(v) = f.attrs.get_mut(fl) {
@@ -2243,6 +2262,155 @@ mod tests {
         }
     }
 
+    // Reference values from PROJ / pyproj 3.7.2 Geod(ellps="WGS84"), in km² and km.
+    #[test]
+    fn calculate_geometry_uses_wgs84_for_polygons() {
+        let rectangle = |x: f64, y: f64, w: f64, h: f64| {
+            vec![Coord::new(x, y), Coord::new(x + w, y), Coord::new(x + w, y + h), Coord::new(x, y + h), Coord::new(x, y)]
+        };
+        for (lat, area, perimeter) in [
+            (0.0, 132.93794889680612, 68.77959742996498),
+            (46.5, 92.07715389612198, 58.58311003354862),
+            (60.0, 67.04587460470582, 52.371257858418176),
+        ] {
+            let parts = vec![
+                vec![rectangle(10.5, lat, 0.1, 0.1), rectangle(10.52, lat + 0.02, 0.02, 0.02)],
+                vec![rectangle(11.0, lat, 0.03, 0.04)],
+            ];
+            for reverse in [false, true] {
+                let mut parts = parts.clone();
+                if reverse {
+                    for ring in parts.iter_mut().flatten() {
+                        ring.reverse();
+                    }
+                }
+                let mut ds = Dataset::new("zones", GeometryType::Polygon, vec![Field::new("code", FieldType::Text)]);
+                ds.push(Some(Geometry::Polygon(parts)), vec![Value::Text("Z1".into())]);
+                let inputs = Map(HashMap::from([("zones".into(), ds.clone())]));
+                let out = layer(crate::run("calculateGeometry", &inputs, &serde_json::json!({"input":"zones"})).unwrap());
+                assert_eq!(out.features[0].geometry, ds.features[0].geometry);
+                assert_eq!(out.value(0, 0), &Value::Text("Z1".into()));
+                let measured_area = out.value(0, out.field_index("area_km2").unwrap()).as_f64().unwrap();
+                let measured_perimeter = out.value(0, out.field_index("perimeter_km").unwrap()).as_f64().unwrap();
+                assert!((measured_area - area).abs() < 1e-6, "{lat}, reversed={reverse}: {measured_area} vs {area}");
+                assert!((measured_perimeter - perimeter).abs() < 1e-6, "{lat}: {measured_perimeter} vs {perimeter}");
+                assert_eq!(inputs.0["zones"], ds);
+            }
+        }
+        let ring = vec![
+            Coord::new(179.9, 10.0),
+            Coord::new(-179.9, 10.0),
+            Coord::new(-179.9, 10.1),
+            Coord::new(179.9, 10.1),
+            Coord::new(179.9, 10.0),
+        ];
+        let mut ds = Dataset::new("island", GeometryType::Polygon, vec![]);
+        ds.push(Some(Geometry::Polygon(vec![vec![ring]])), vec![]);
+        let inputs = Map(HashMap::from([("island".into(), ds)]));
+        let out = layer(crate::run("calculateGeometry", &inputs, &serde_json::json!({"input":"island"})).unwrap());
+        assert!((out.value(0, out.field_index("area_km2").unwrap()).as_f64().unwrap() - 242.50305444813918).abs() < 1e-6);
+        assert!((out.value(0, out.field_index("perimeter_km").unwrap()).as_f64().unwrap() - 65.97062684618793).abs() < 1e-6);
+    }
+
+    #[test]
+    fn calculate_geometry_keeps_large_regions_and_reversed_small_rings() {
+        let ring = |points: &[(f64, f64)]| points.iter().map(|&(x, y)| Coord::new(x, y)).collect::<Vec<_>>();
+        let band = ring(&[
+            (-170.0, -60.0),
+            (-85.0, -60.0),
+            (0.0, -60.0),
+            (85.0, -60.0),
+            (170.0, -60.0),
+            (170.0, 60.0),
+            (85.0, 60.0),
+            (0.0, 60.0),
+            (-85.0, 60.0),
+            (-170.0, 60.0),
+            (-170.0, -60.0),
+        ]);
+        let hole = ring(&[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]);
+        // GeographicLib 2.1 WGS84 Polygon.Compute: unsigned for the band, signed absolute for caps.
+        let cases = [
+            ("large", vec![band, hole], 435312377.1122553 - 12308.778361469453, 61868.48818267019 + 443.77091724830194),
+            (
+                "north",
+                vec![ring(&[(-135.0, 80.0), (-45.0, 80.0), (45.0, 80.0), (135.0, 80.0), (-135.0, 80.0)])],
+                2507270.031169875,
+                6301.599963614222,
+            ),
+            (
+                "south",
+                vec![ring(&[(-135.0, -80.0), (-45.0, -80.0), (45.0, -80.0), (135.0, -80.0), (-135.0, -80.0)])],
+                2507270.031169875,
+                6301.599963614222,
+            ),
+            (
+                "dateline",
+                vec![ring(&[(179.9, 10.0), (-179.9, 10.0), (-179.9, 10.1), (179.9, 10.1), (179.9, 10.0)])],
+                242.50305444813918,
+                65.97062684618793,
+            ),
+        ];
+        for (name, rings, area, perimeter) in cases {
+            for reverse in [false, true] {
+                let mut rings = rings.clone();
+                if reverse {
+                    for ring in &mut rings {
+                        ring.reverse();
+                    }
+                }
+                let geometry = Geometry::Polygon(vec![rings]);
+                let mut ds = Dataset::new(name, GeometryType::Polygon, vec![]);
+                ds.push(Some(geometry.clone()), vec![]);
+                let inputs = Map(HashMap::from([(name.into(), ds)]));
+                let out = layer(crate::run("calculateGeometry", &inputs, &serde_json::json!({"input":name})).unwrap());
+                let measured = out.value(0, out.field_index("area_km2").unwrap()).as_f64().unwrap();
+                let per = out.value(0, out.field_index("perimeter_km").unwrap()).as_f64().unwrap();
+                assert!((measured - area).abs() < 1e-6, "{name}, reversed={reverse}: {measured} vs {area}");
+                assert!((per - perimeter).abs() < 1e-6, "{name}, reversed={reverse}: {per} vs {perimeter}");
+                assert_eq!(out.features[0].geometry.as_ref(), Some(&geometry));
+            }
+        }
+    }
+
+    #[test]
+    fn calculate_geometry_does_not_orient_invalid_coordinates() {
+        for bad in [Coord::new(f64::NAN, 0.0), Coord::new(f64::INFINITY, 0.0), Coord::new(0.0, f64::NAN), Coord::new(0.0, 100.0)] {
+            let mut ds = Dataset::new("invalid", GeometryType::Polygon, vec![]);
+            ds.push(
+                Some(Geometry::Polygon(vec![vec![vec![
+                    Coord::new(0.0, 0.0),
+                    Coord::new(1.0, 0.0),
+                    bad,
+                    Coord::new(0.0, 1.0),
+                    Coord::new(0.0, 0.0),
+                ]]])),
+                vec![],
+            );
+            let inputs = Map(HashMap::from([("invalid".into(), ds)]));
+            assert!(crate::run("calculateGeometry", &inputs, &serde_json::json!({"input":"invalid"})).is_ok());
+        }
+    }
+
+    #[test]
+    fn calculate_geometry_uses_wgs84_for_lines() {
+        let mut ds = Dataset::new("routes", GeometryType::Polyline, vec![]);
+        let parts = vec![
+            vec![Coord::new(0.0, 0.0), Coord::new(1.0, 0.0)],
+            vec![Coord::new(10.0, 60.0), Coord::new(10.0, 61.0)],
+            vec![Coord::new(179.9, 10.0), Coord::new(-179.9, 10.0)],
+        ];
+        ds.push(Some(Geometry::Line(parts)), vec![]);
+        ds.push(None, vec![]);
+        let inputs = Map(HashMap::from([("routes".into(), ds.clone())]));
+        let out = layer(crate::run("calculateGeometry", &inputs, &serde_json::json!({"input":"routes"})).unwrap());
+        let length = out.field_index("length_km").unwrap();
+        let expected = 111.31949079327356 + 111.42072786973769 + 21.92787247793737;
+        assert!((out.value(0, length).as_f64().unwrap() - expected).abs() < 1e-6);
+        assert_eq!(out.value(1, length).as_f64(), Some(0.0));
+        assert_eq!(out.features[0].geometry, ds.features[0].geometry);
+    }
+
     /// Calculate Geometry stores full-precision values (#513): rounding to 0.01 km² gave every
     /// polygon under half a hectare an area of 0, and lines under 5 m a length of 0.
     #[test]

crates/engine/src/tests.rs

@@ -17395,3 +17395,26 @@ fn publication_locator_keeps_small_islands_across_the_antimeridian() {
         assert_eq!(geometry.as_ref(), Some(&original), "island at {lon}, {lat}");
     }
 }
+
+#[test]
+fn calculate_geometry_updates_ellipsoidal_lengths_in_one_undo_step() {
+    use rustgis_data::{Dataset, Field, FieldType, Geometry, GeometryType, Value};
+    use rustgis_geom::Coord;
+    let mut s = world();
+    let mut routes = Dataset::new("routes", GeometryType::Polyline, vec![Field::new("code", FieldType::Text)]);
+    routes.push(Some(Geometry::Line(vec![vec![Coord::new(0.0, 0.0), Coord::new(1.0, 0.0)]])), vec![Value::Text("R1".into())]);
+    s.map_mut().add_layer(rustgis_map::Layer::new(0, routes.clone()));
+    let depth = s.undo_depth();
+    s.run("gp.run", &json!({"tool":"calculateGeometry", "params":{"input":"routes"}})).unwrap();
+    let out = &s.map().layer_by_name("routes").unwrap().dataset;
+    let length = out.field_index("length_km").unwrap();
+    assert!((out.value(0, length).as_f64().unwrap() - 111.31949079327356).abs() < 1e-6);
+    assert_eq!(out.features[0].geometry, routes.features[0].geometry);
+    assert_eq!(out.value(0, 0), &Value::Text("R1".into()));
+    assert_eq!(s.undo_depth(), depth + 1);
+    assert!(s.undo());
+    assert_eq!(s.map().layer_by_name("routes").unwrap().dataset, routes);
+    assert!(s.redo());
+    let out = &s.map().layer_by_name("routes").unwrap().dataset;
+    assert!((out.value(0, length).as_f64().unwrap() - 111.31949079327356).abs() < 1e-6);
+}

docs/user-guide/analysis.md

@@ -123,9 +123,16 @@ projection. Each feature is buffered in an azimuthal equidistant projection on t
 centred on it, so a point's buffer is exact; long lines and large polygon outlines are
 buffered in stretches of a few hundred kilometres to limit projection error. Negative Buffer
 distances shrink polygons. Buffer and Select Layer By Location limit distances to 5,000 km.
-Lengths and areas (Calculate Geometry Attributes, and Spatial Join's **Largest overlap**) are
-still measured on a sphere of radius 6,371 km: on the ellipsoid, lengths differ by up to about
-0.6 % and areas by up to about 0.9 % (a small area near a pole).
+**Calculate Geometry Attributes** measures polygon areas and perimeters, and line lengths,
+on the WGS 84 ellipsoid, independently of the map's display projection. Holes are subtracted
+regardless of ring winding, and multipart measurements are summed. Ring directions are
+normalized in continuous longitude space, so non-polar regions larger than half the Earth
+keep their full area and small antimeridian crossings stay small. Rings winding around a pole
+use the smaller polar region in either direction; reversing a ring does not select its global
+complement. Edges follow the shortest geodesics between vertices. The output fields remain
+`area_km2`, `perimeter_km` and `length_km`; point coordinates are unchanged.
+Spatial Join's **Largest overlap** still measures area on a sphere of radius 6,371 km; it can
+differ from the ellipsoid by up to about 0.9 % (a small area near a pole).
 
 **IDW Interpolation** works like the QGIS tool of the same name: every point weighs on every
 cell by 1/distanceᴾ, where P is the distance coefficient (2 by default). Lines and polygons