From 6069ab1968baeb527214e4467d40d7231f9f1f87 Mon Sep 17 00:00:00 2001 From: choupeanut Date: Mon, 26 Oct 2026 15:13:23 +0800 Subject: [PATCH 1/2] fix(analysis): calculate geometry on the WGS 84 ellipsoid --- crates/analysis/src/tools.rs | 81 +++++++++++++++++++++++++++++++++--- crates/engine/src/tests.rs | 23 ++++++++++ docs/user-guide/analysis.md | 9 ++-- 3 files changed, 105 insertions(+), 8 deletions(-) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index b749c6724..0a142d4bd 100644 --- a/crates/analysis/src/tools.rs +++ b/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 { } fn calculate_geometry(inp: &dyn Inputs, p: &Json) -> Result { + 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,9 @@ fn calculate_geometry(inp: &dyn Inputs, p: &Json) -> Result 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); - a += if i == 0 { ra } else { -ra }; - per += geodesic_length_m(r); + let (rp, ra) = geo::Polygon::new(line(r), vec![]).geodesic_perimeter_area_signed(); + a += if i == 0 { ra.abs() } else { -ra.abs() }; + per += rp; } } } @@ -1265,7 +1267,7 @@ fn calculate_geometry(inp: &dyn Inputs, p: &Json) -> Result 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 +2245,75 @@ 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_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] diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index def24f2b7..907bb9dc4 100644 --- a/crates/engine/src/tests.rs +++ b/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); +} diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 0604740cc..6646e01eb 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -123,9 +123,12 @@ 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. 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 From 242886f01f15ccc8245130c8e420ab51bbe8fd8c Mon Sep 17 00:00:00 2001 From: choupeanut Date: Mon, 26 Oct 2026 17:42:46 +0800 Subject: [PATCH 2/2] fix(analysis): preserve large geodesic polygon areas --- crates/analysis/src/tools.rs | 101 ++++++++++++++++++++++++++++++++++- docs/user-guide/analysis.md | 6 ++- 2 files changed, 104 insertions(+), 3 deletions(-) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index 0a142d4bd..956361a52 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -1249,8 +1249,25 @@ fn calculate_geometry(inp: &dyn Inputs, p: &Json) -> Result if let Some(Geometry::Polygon(polys)) = &f.geometry { for rings in polys { for (i, r) in rings.iter().enumerate() { - let (rp, ra) = geo::Polygon::new(line(r), vec![]).geodesic_perimeter_area_signed(); - a += if i == 0 { ra.abs() } else { -ra.abs() }; + 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 += rp; } } @@ -2295,6 +2312,86 @@ mod tests { 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::>(); + 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![]); diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 6646e01eb..9e24ca0c9 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -125,7 +125,11 @@ buffered in stretches of a few hundred kilometres to limit projection error. Neg distances shrink polygons. Buffer and Select Layer By Location limit distances to 5,000 km. **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. The output fields remain +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).