Measure on the WGS 84 ellipsoid in the attribute table, measure tool and network tools (#939)

crates/analysis/src/network.rs

@@ -4,7 +4,7 @@
 use std::collections::{BinaryHeap, HashMap};
 
 use rustgis_data::{Dataset, Geometry};
-use rustgis_geom::{Coord, haversine_m};
+use rustgis_geom::{Coord, EARTH_RADIUS_M, WGS84_MIN_RADIUS_M, haversine_m, wgs84_distance_m};
 
 /// A shortest route: its path, length in metres, cost, and the line rows it uses.
 #[derive(Clone, Debug, PartialEq)]
@@ -13,6 +13,8 @@ pub struct Route {
     pub length_m: f64,
     pub cost: f64,
     pub rows: Vec<usize>,
+    /// Segments omitted because an endpoint has non-finite coordinates or an invalid latitude.
+    pub skipped_segments: usize,
 }
 
 #[derive(PartialEq)]
@@ -35,8 +37,8 @@ fn key(c: Coord) -> (i64, i64) {
 }
 
 /// The shortest route between two points over the lines of `rows`, joining where vertices
-/// coincide. Each edge costs its length times the row's `cost_field` value (1 when empty; lines
-/// costing zero or less are closed). Each
+/// coincide. Each edge costs its length (on the WGS 84 ellipsoid) times the row's `cost_field`
+/// value (1 when empty; lines costing zero or less are closed). Each
 /// end snaps to the nearest point on a line (within `snap_m`), where the route starts or stops.
 pub fn shortest_route(
     ds: &Dataset,
@@ -46,6 +48,21 @@ pub fn shortest_route(
     cost_field: Option<usize>,
     snap_m: f64,
 ) -> Result<Route, String> {
+    if !snap_m.is_finite() {
+        return Err("the snap distance must be finite".into());
+    }
+    if [from, to].iter().any(|c| !c.is_finite() || c.y.abs() > 90.0) {
+        return Err("a route end has invalid coordinates".into());
+    }
+    // Preserve clean-network errors; failed partial routes must still disclose omissions.
+    let route_error = |reason: String, skipped_segments: usize| {
+        if skipped_segments == 0 {
+            reason
+        } else {
+            format!("{reason}; {skipped_segments} segments with invalid coordinates were omitted from '{}'", ds.name)
+        }
+    };
+    let mut skipped_segments = 0;
     let mut ids: HashMap<(i64, i64), usize> = HashMap::new();
     let mut nodes: Vec<Coord> = Vec::new();
     let mut adj: Vec<Vec<(usize, f64, f64, usize)>> = Vec::new(); // (to, cost, length, row)
@@ -69,8 +86,14 @@ pub fn shortest_route(
         for part in parts {
             for w in part.windows(2) {
                 let (Some(a), Some(b)) = (w.first(), w.get(1)) else { continue };
+                // Non-finite coordinates cast to ordinary integer keys (NaN becomes zero).
+                // Keep unroutable segments out before they can corrupt a valid shared node.
+                if !a.is_finite() || !b.is_finite() || a.y.abs() > 90.0 || b.y.abs() > 90.0 {
+                    skipped_segments += 1;
+                    continue;
+                }
                 let (ia, ib) = (node(*a, &mut ids, &mut nodes, &mut adj), node(*b, &mut ids, &mut nodes, &mut adj));
-                let len = haversine_m(*a, *b);
+                let len = wgs84_distance_m(*a, *b);
                 if let Some(v) = adj.get_mut(ia) {
                     v.push((ib, len * weight, len, *r));
                 }
@@ -80,11 +103,14 @@ pub fn shortest_route(
                 segs.push((ia, ib, weight, *r));
                 // Checked per segment, so one huge feature can't grow the network unbounded.
                 if segs.len() > 10_000_000 || nodes.len() > 5_000_000 {
-                    return Err("the network is too large".into());
+                    return Err(route_error("the network is too large".into(), skipped_segments));
                 }
             }
         }
     }
+    if segs.is_empty() && skipped_segments > 0 {
+        return Err(route_error("no usable network remains".into(), skipped_segments));
+    }
     // The closest point on any segment: (segment, fraction along it, the point, metres away).
     let nearest = |p: Coord, nodes: &[Coord]| -> Option<(usize, f64, Coord, f64)> {
         let k = p.y.to_radians().cos().max(0.01);
@@ -98,15 +124,31 @@ pub fn shortest_route(
                 let q = Coord::new(a.x + t * (b.x - a.x), a.y + t * (b.y - a.y));
                 Some((i, t, q, haversine_m(q, p)))
             })
-            .min_by(|a, b| a.3.partial_cmp(&b.3).unwrap_or(std::cmp::Ordering::Equal))
+            .fold(None, |best: Option<(usize, f64, Coord, f64)>, (i, t, q, sphere)| {
+                // The minimum WGS 84 radius gives a lower bound, so a farther spherical
+                // candidate cannot be discarded merely because the two metrics rank it differently.
+                let lower = if (-90.0..=90.0).contains(&p.y) && (-90.0..=90.0).contains(&q.y) {
+                    sphere * WGS84_MIN_RADIUS_M / EARTH_RADIUS_M
+                } else {
+                    0.0 // The inverse solver clamps latitudes; the spherical bound would not apply.
+                };
+                if best.as_ref().is_some_and(|b| lower > b.3) {
+                    return best;
+                }
+                let distance = wgs84_distance_m(q, p);
+                if distance.is_finite() && best.as_ref().is_none_or(|b| distance < b.3) { Some((i, t, q, distance)) } else { best }
+            })
     };
     // Each end becomes a node on its segment, joined to the segment's two ends.
     let attach = |p: Coord, nodes: &mut Vec<Coord>, adj: &mut Vec<Vec<(usize, f64, f64, usize)>>| -> Result<(usize, usize, f64), String> {
-        let (si, t, q, d) = nearest(p, nodes).ok_or("the layer has no lines")?;
+        let (si, t, q, d) = nearest(p, nodes).ok_or_else(|| route_error("the layer has no lines".into(), skipped_segments))?;
         if d > snap_m {
-            return Err(format!("a route end is more than {} from the network", rustgis_geom::format_distance(snap_m)));
+            return Err(route_error(
+                format!("a route end is more than {} from the network", rustgis_geom::format_distance(snap_m)),
+                skipped_segments,
+            ));
         }
-        let &(a, b, weight, row) = segs.get(si).ok_or("the layer has no lines")?;
+        let &(a, b, weight, row) = segs.get(si).ok_or_else(|| route_error("the layer has no lines".into(), skipped_segments))?;
         if t <= 0.0 {
             return Ok((a, si, t));
         }
@@ -117,7 +159,7 @@ pub fn shortest_route(
         adj.push(Vec::new());
         let n = nodes.len() - 1;
         for end in [a, b] {
-            let len = haversine_m(q, nodes.get(end).copied().unwrap_or(q));
+            let len = wgs84_distance_m(q, nodes.get(end).copied().unwrap_or(q));
             if let Some(v) = adj.get_mut(n) {
                 v.push((end, len * weight, len, row));
             }
@@ -132,9 +174,9 @@ pub fn shortest_route(
     // Both ends on one segment: they're joined directly along it.
     if s_seg == t_seg && s != t && (0.0..=1.0).contains(&s_t) && (0.0..=1.0).contains(&t_t) {
         let (Some(&(_, _, weight, row)), Some(a), Some(b)) = (segs.get(s_seg), nodes.get(s).copied(), nodes.get(t).copied()) else {
-            return Err("the layer has no lines".into());
+            return Err(route_error("the layer has no lines".into(), skipped_segments));
         };
-        let len = haversine_m(a, b);
+        let len = wgs84_distance_m(a, b);
         if let Some(v) = adj.get_mut(s) {
             v.push((t, len * weight, len, row));
         }
@@ -172,7 +214,7 @@ pub fn shortest_route(
     }
     let cost = dist.get(t).copied().unwrap_or(f64::INFINITY);
     if !cost.is_finite() {
-        return Err("the two points aren't connected by the network".into());
+        return Err(route_error("the two points aren't connected by the network".into(), skipped_segments));
     }
     let mut path = Vec::new();
     let mut length = 0.0;
@@ -190,7 +232,7 @@ pub fn shortest_route(
     path.reverse();
     used.reverse();
     used.dedup();
-    Ok(Route { path, length_m: length, cost, rows: used })
+    Ok(Route { path, length_m: length, cost, rows: used, skipped_segments })
 }
 
 /// Most random points made in one call, and point-in-polygon tries spent making them.
@@ -445,9 +487,49 @@ mod tests {
         ds.push(Some(Geometry::Line(vec![vec![Coord::new(0.0, 0.0), Coord::new(0.09, 0.0)]])), vec![]);
         // Clicks beside its middle, far from both vertices but close to the road.
         let r = shortest_route(&ds, &[0], Coord::new(0.03, 0.0001), Coord::new(0.06, -0.0001), None, 50.0).unwrap();
-        assert!((r.length_m - 3336.0).abs() < 20.0, "{}", r.length_m);
+        assert!((r.length_m - 3339.6).abs() < 1.0, "{}", r.length_m);
         assert!((r.path[0].x - 0.03).abs() < 1e-9 && (r.path.last().unwrap().x - 0.06).abs() < 1e-9);
         assert!(shortest_route(&ds, &[0], Coord::new(0.03, 0.01), Coord::new(0.06, 0.0), None, 50.0).is_err(), "1 km away");
+        let error = shortest_route(&ds, &[0], Coord::new(f64::NAN, 0.0), Coord::new(0.06, 0.0), None, 50.0).unwrap_err();
+        assert!(error.contains("route end"), "{error}");
+        let error = shortest_route(&ds, &[0], Coord::new(0.0, 95.0), Coord::new(0.06, 0.0), None, 50.0).unwrap_err();
+        assert!(error.contains("route end"), "{error}");
+        assert!(shortest_route(&ds, &[0], Coord::new(0.0, 1.0), Coord::new(0.06, 1.0), None, f64::NAN).is_err());
+    }
+
+    #[test]
+    fn a_broken_segment_listed_first_does_not_hide_a_valid_one() {
+        // The invalid segment is excluded before graph construction and cannot hide the road.
+        let mut ds = Dataset::new("roads", GeometryType::Polyline, vec![]);
+        ds.push(Some(Geometry::Line(vec![vec![Coord::new(0.5, 0.5), Coord::new(f64::NAN, 0.5)]])), vec![]);
+        ds.push(Some(Geometry::Line(vec![vec![Coord::new(0.0, 0.0), Coord::new(0.09, 0.0)]])), vec![]);
+        let r = shortest_route(&ds, &[0, 1], Coord::new(0.03, 0.0001), Coord::new(0.06, -0.0001), None, 50.0).unwrap();
+        assert!((r.length_m - 3339.6).abs() < 1.0, "{}", r.length_m);
+    }
+
+    #[test]
+    fn nonfinite_vertices_do_not_alias_valid_network_nodes() {
+        let mut ds = Dataset::new("roads", GeometryType::Polyline, vec![]);
+        ds.push(Some(Geometry::Line(vec![vec![Coord::new(f64::NAN, 0.5), Coord::new(0.5, 0.5)]])), vec![]);
+        ds.push(Some(Geometry::Line(vec![vec![Coord::new(0.0, 0.5), Coord::new(0.09, 0.5)]])), vec![]);
+        let r = shortest_route(&ds, &[0, 1], Coord::new(0.0, 0.5), Coord::new(0.09, 0.5), None, 1.0).unwrap();
+        assert!(r.path.iter().all(|c| c.is_finite()), "{:?}", r.path);
+        assert_eq!(r.rows, vec![1]);
+        assert!((r.length_m - 10_018.375).abs() < 1.0, "{}", r.length_m);
+    }
+
+    #[test]
+    fn snap_acceptance_uses_the_nearest_ellipsoidal_candidate() {
+        let mut ds = Dataset::new("roads", GeometryType::Polyline, vec![]);
+        // At the equator, 0.001 degrees east is nearer on the sphere than 0.001003
+        // degrees north. On WGS 84 only the northern segment is within 111 metres.
+        ds.push(Some(Geometry::Line(vec![vec![Coord::new(0.001, 0.0), Coord::new(0.002, 0.0)]])), vec![]);
+        ds.push(Some(Geometry::Line(vec![vec![Coord::new(0.0, 0.001003), Coord::new(0.0, 0.003)]])), vec![]);
+        let from = Coord::new(0.0, 0.0);
+        assert!(haversine_m(from, Coord::new(0.001, 0.0)) < haversine_m(from, Coord::new(0.0, 0.001003)));
+        let r = shortest_route(&ds, &[0, 1], from, Coord::new(0.0, 0.003), None, 111.0).unwrap();
+        assert_eq!(r.rows, vec![1]);
+        assert_eq!(r.path.first(), Some(&Coord::new(0.0, 0.001003)));
     }
 
     fn offset(at: Coord, x: f64, y: f64) -> Coord {

crates/engine/src/cmd/workflows.rs

@@ -678,6 +678,7 @@ fn interpolate(s: &mut Session, p: &Value) -> CmdResult {
     add_raster(s, r)
 }
 
+/// Creates a route layer and reports invalid source segments omitted from its network.
 fn route(s: &mut Session, p: &Value) -> CmdResult {
     let id = s.layer_id(p)?;
     let l = s.map().layer(id).ok_or("no such layer")?;
@@ -694,11 +695,18 @@ fn route(s: &mut Session, p: &Value) -> CmdResult {
         cost,
         f(p, "snapDistance").unwrap_or(1000.0).clamp(0.0, 1e6),
     )?;
+    let source_name = l.name.clone();
     let name = out_name(p, format!("Route on {}", l.name));
     let mut ds =
         Dataset::new(name, GeometryType::Polyline, vec![Field::new("length_km", FieldType::Double), Field::new("cost", FieldType::Double)]);
     ds.push(Some(Geometry::Line(vec![r.path.clone()])), vec![Attr::Float(r.length_m / 1000.0), Attr::Float(r.cost)]);
     let nid = super::layer::add_dataset(s, ds);
+    if r.skipped_segments > 0 {
+        s.notify(
+            format!("{source_name}: {} segments with invalid coordinates were omitted from the route network", r.skipped_segments),
+            false,
+        );
+    }
     Ok(json!({"layer": nid, "lengthKm": r.length_m / 1000.0, "cost": r.cost, "segments": r.rows.len()}))
 }
 

crates/engine/src/geometry_calc.rs

@@ -2,7 +2,7 @@
 //! X/Y coordinate, calculated from its WGS84 geometry, in a field.
 use crate::{CmdResult, Session};
 use rustgis_data::{FieldType, Geometry, GeometryType, Value, raster::Crs};
-use rustgis_geom::{geodesic_length_m, geodesic_ring_area_m2, spherical_centre};
+use rustgis_geom::{spherical_centre, wgs84_length_m, wgs84_ring_perimeter_area_m2};
 use serde_json::{Value as Json, json};
 
 /// The properties a layer of `kind` can calculate, as `(key, label)` pairs.
@@ -32,25 +32,21 @@ pub fn units(property: &str) -> &'static [(&'static str, &'static str)] {
     }
 }
 
+/// Measures WGS 84 geometry in metres or square metres; holes subtract from polygon area.
 fn measure(g: &Geometry, property: &str) -> Option<f64> {
     match (g, property) {
+        // Full-turn and polar rings retain the documented exceptions to geoprocessing parity:
+        // docs/user-guide/analysis.md (Calculate Geometry Attributes).
         (Geometry::Polygon(polys), "area") => Some(
             polys
                 .iter()
-                .map(|rings| rings.iter().enumerate().map(|(i, r)| geodesic_ring_area_m2(r) * if i == 0 { 1.0 } else { -1.0 }).sum::<f64>())
-                .sum(),
-        ),
-        (Geometry::Polygon(polys), "perimeter") => Some(
-            polys
-                .iter()
-                .flatten()
-                .map(|r| {
-                    let closing = r.first().zip(r.last()).filter(|(a, b)| a != b).map(|(a, b)| geodesic_length_m(&[*b, *a])).unwrap_or(0.0);
-                    geodesic_length_m(r) + closing
+                .map(|rings| {
+                    rings.iter().enumerate().map(|(i, r)| wgs84_ring_perimeter_area_m2(r).1 * if i == 0 { 1.0 } else { -1.0 }).sum::<f64>()
                 })
                 .sum(),
         ),
-        (Geometry::Line(parts), "length") => Some(parts.iter().map(|r| geodesic_length_m(r)).sum()),
+        (Geometry::Polygon(polys), "perimeter") => Some(polys.iter().flatten().map(|r| wgs84_ring_perimeter_area_m2(r).0).sum()),
+        (Geometry::Line(parts), "length") => Some(parts.iter().map(|r| wgs84_length_m(r)).sum()),
         _ => None,
     }
 }

crates/engine/tests/geometry_calculator.rs

@@ -28,7 +28,7 @@ fn selected_units_overwrite_and_undo() {
     p["overwrite"] = json!(true);
     assert_eq!(s.run("edit.calculateGeometry", &p).unwrap()["rows"], 1);
     let value = s.map().layer(id).unwrap().dataset.value(1, 0).as_f64().unwrap();
-    assert!((222.0..223.0).contains(&value));
+    assert!((222.6..222.7).contains(&value), "{value}"); // the ellipsoid: the sphere gives 222.39
     assert_eq!(s.map().layer(id).unwrap().dataset.value(0, 0).as_f64(), Some(99.0));
     s.run("edit.undo", &json!({})).unwrap();
     assert_eq!(s.map().layer(id).unwrap().dataset.value(1, 0).as_f64(), Some(99.0));
@@ -49,6 +49,22 @@ fn selected_units_overwrite_and_undo() {
     assert!(s.map().layer(id).unwrap().dataset.field_index("metres").is_none());
 }
 
+#[test]
+fn areas_are_on_the_wgs84_ellipsoid() {
+    // An equator 1° × 1° box measured by GeographicLib itself (pyproj.Geod / Polygon.Compute on
+    // WGS 84); the attribute table's Calculate Geometry values are checked against those (#939).
+    let ring = vec![Coord::new(0.0, 0.0), Coord::new(1.0, 0.0), Coord::new(1.0, 1.0), Coord::new(0.0, 1.0)];
+    let (mut s, id) = session(GeometryType::Polygon, vec![Some(Geometry::Polygon(vec![vec![ring]]))]);
+    s.run("edit.calculateGeometry", &json!({"layer":id,"field":"result","overwrite":true,"property":"area","unit":"m2","scope":"all"}))
+        .unwrap();
+    let actual = s.map().layer(id).unwrap().dataset.value(0, 0).as_f64().unwrap();
+    assert!((actual - 12_308_778_361.469).abs() < 0.001, "{actual}");
+    s.run("edit.calculateGeometry", &json!({"layer":id,"field":"result","overwrite":true,"property":"perimeter","unit":"m","scope":"all"}))
+        .unwrap();
+    let actual = s.map().layer(id).unwrap().dataset.value(0, 0).as_f64().unwrap();
+    assert!((actual - 443_770.917).abs() < 0.001, "{actual}");
+}
+
 #[test]
 fn polygon_holes_multipart_and_area_units() {
     let ring = |x: f64, size: f64| {
@@ -57,8 +73,9 @@ fn polygon_holes_multipart_and_area_units() {
     let outer = ring(0.0, 0.01);
     let hole = ring(0.002, 0.002);
     let part = ring(0.02, 0.005);
-    let expected = rustgis_engine::geom::geodesic_ring_area_m2(&outer) - rustgis_engine::geom::geodesic_ring_area_m2(&hole)
-        + rustgis_engine::geom::geodesic_ring_area_m2(&part);
+    let expected = rustgis_engine::geom::wgs84_ring_perimeter_area_m2(&outer).1
+        - rustgis_engine::geom::wgs84_ring_perimeter_area_m2(&hole).1
+        + rustgis_engine::geom::wgs84_ring_perimeter_area_m2(&part).1;
     let (mut s, id) = session(GeometryType::Polygon, vec![Some(Geometry::Polygon(vec![vec![outer, hole], vec![part]]))]);
     for (unit, divisor) in [("m2", 1.0), ("km2", 1e6), ("ha", 1e4), ("ac", 4046.8564224)] {
         s.run("edit.calculateGeometry", &json!({"layer":id,"field":"result","overwrite":true,"property":"area","unit":unit,"scope":"all"}))

crates/engine/tests/workflows.rs

@@ -348,3 +348,45 @@ fn shelby_county_tract_areas_match_the_census() {
     reopened.run("share.exportLayer", &json!({"layer": "shelby", "path": gpkg})).unwrap();
     check(&io::read_path(&gpkg).unwrap().0);
 }
+
+#[test]
+fn network_routes_report_invalid_segments_but_clean_networks_are_quiet() {
+    use rustgis_engine::data::{Geometry, GeometryType};
+    use rustgis_engine::geom::Coord;
+    for invalid in [false, true] {
+        let mut s = Session::new(Project::default());
+        let mut ds = Dataset::new("roads", GeometryType::Polyline, vec![]);
+        ds.push(
+            Some(Geometry::Line(vec![
+                vec![Coord::new(0.0, 0.0), Coord::new(0.01, 0.0)],
+                vec![Coord::new(0.02, 0.0), Coord::new(0.03, 0.0)],
+            ])),
+            vec![],
+        );
+        if invalid {
+            ds.push(Some(Geometry::Line(vec![vec![Coord::new(f64::NAN, 0.0), Coord::new(0.02, 0.0)]])), vec![]);
+            ds.push(Some(Geometry::Line(vec![vec![Coord::new(0.0, 95.0), Coord::new(0.02, 95.0)]])), vec![]);
+        }
+        let id = s.map_mut().add_dataset(ds);
+        let mark = s.notice_mark();
+        s.run("analysis.route", &json!({"layer":id, "from":[0.0,0.0], "to":[0.01,0.0], "snapDistance":1})).unwrap();
+        let notices = s.notices_since(mark);
+        if invalid {
+            assert_eq!(notices.len(), 1, "{notices:?}");
+            assert!(notices[0].text.contains("2 segments") && notices[0].text.contains("roads"), "{notices:?}");
+            assert!(!notices[0].error);
+        } else {
+            assert!(notices.is_empty(), "{notices:?}");
+        }
+        for (to, expected) in
+            [(0.03, "the two points aren't connected by the network"), (0.1, "a route end is more than 1 m from the network")]
+        {
+            let error = s.run("analysis.route", &json!({"layer":id, "from":[0.0,0.0], "to":[to,0.0], "snapDistance":1})).unwrap_err();
+            if invalid {
+                assert!(error.contains("2 segments") && error.contains("roads"), "{error}");
+            } else {
+                assert_eq!(error, expected);
+            }
+        }
+    }
+}