From 686a381de7f378f8743d3b31971a722d8759903a Mon Sep 17 00:00:00 2001 From: subarna-kandel Date: Mon, 26 Oct 2026 15:30:20 +0000 Subject: [PATCH 01/12] measure on the WGS 84 ellipsoid everywhere lengths and areas are named geodesic The attribute table's Calculate Geometry (edit.calculateGeometry), the map's Measure tools on the Earth and the network tools' edge and route lengths still measured on the mean sphere (R = 6,371,008.8 m), after #905 and #928 made the buffers, joins and the geoprocessing Calculate Geometry ellipsoidal. The same polygon then returned two areas about 0.26 % apart depending on which Calculate Geometry ran (#939). rustgis-geom gains wgs84_length_m (inverse geodesic solves, one per segment) and wgs84_ring_perimeter_area_m2 (GeographicLib's PolygonArea, the ring closed along a geodesic, unsigned in either winding), both stopping at non-finite coordinates as haversine_m now does. The three paths above use them; on the Moon or Mars the Measure tools keep measuring on that body's sphere and scale, as before. Network tests keep passing without changes to their tolerances. Fixes #939. --- crates/analysis/src/network.rs | 14 +- crates/engine/src/geometry_calc.rs | 19 +-- crates/engine/tests/geometry_calculator.rs | 21 ++- crates/geom/src/lib.rs | 141 ++++++++++++++++++++- crates/ui-egui/src/mapview.rs | 24 ++-- docs/user-guide/analysis.md | 6 +- docs/user-guide/data.md | 3 +- docs/user-guide/editing.md | 4 +- 8 files changed, 199 insertions(+), 33 deletions(-) diff --git a/crates/analysis/src/network.rs b/crates/analysis/src/network.rs index 0d79ca473..56d79fa78 100644 --- a/crates/analysis/src/network.rs +++ b/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, wgs84_distance_m}; /// A shortest route: its path, length in metres, cost, and the line rows it uses. #[derive(Clone, Debug, PartialEq)] @@ -35,8 +35,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, @@ -70,7 +70,7 @@ pub fn shortest_route( for w in part.windows(2) { let (Some(a), Some(b)) = (w.first(), w.get(1)) else { 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)); } @@ -96,7 +96,7 @@ pub fn shortest_route( let len2 = dx * dx + dy * dy; let t = if len2 > 0.0 { (((p.x - a.x) * k * dx + (p.y - a.y) * dy) / len2).clamp(0.0, 1.0) } else { 0.0 }; 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))) + Some((i, t, q, wgs84_distance_m(q, p))) }) .min_by(|a, b| a.3.partial_cmp(&b.3).unwrap_or(std::cmp::Ordering::Equal)) }; @@ -117,7 +117,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)); } @@ -134,7 +134,7 @@ pub fn shortest_route( 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()); }; - 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)); } diff --git a/crates/engine/src/geometry_calc.rs b/crates/engine/src/geometry_calc.rs index 73df963f9..149bfe419 100644 --- a/crates/engine/src/geometry_calc.rs +++ b/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. @@ -34,23 +34,18 @@ pub fn units(property: &str) -> &'static [(&'static str, &'static str)] { fn measure(g: &Geometry, property: &str) -> Option { match (g, property) { + // On the WGS 84 ellipsoid, like every other measurement named geodesic (#939): the same + // numbers the geoprocessing tools give. (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::()) - .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::() }) .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, } } diff --git a/crates/engine/tests/geometry_calculator.rs b/crates/engine/tests/geometry_calculator.rs index c957ea391..8226ea5b0 100644 --- a/crates/engine/tests/geometry_calculator.rs +++ b/crates/engine/tests/geometry_calculator.rs @@ -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), so the attribute table and the geoprocessing tool agree with each other (#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"})) diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 93ad1686f..eb1eecfeb 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -13,7 +13,7 @@ pub use sinusoidal::Sinusoidal; use std::sync::LazyLock; -use geographiclib_rs::{DirectGeodesic, Geodesic, InverseGeodesic}; +use geographiclib_rs::{DirectGeodesic, Geodesic, InverseGeodesic, PolygonArea, Winding}; use serde::{Deserialize, Serialize}; /// Mean Earth radius used by the haversine formula (metres). @@ -243,10 +243,66 @@ pub fn geodesic_destination(c: Coord, azimuth_deg: f64, dist_m: f64) -> Coord { /// The geodesic distance between two lon/lat points on the WGS 84 ellipsoid (metres): the inverse /// geodesic problem, solved by GeographicLib. pub fn wgs84_distance_m(a: Coord, b: Coord) -> f64 { + // A broken coordinate is no place (see [`haversine_m`]): the clamp below would otherwise turn a + // latitude of infinity into a pole. + if !a.is_finite() || !b.is_finite() { + return f64::NAN; + } let (s12, _): (f64, f64) = WGS84.inverse(a.y.clamp(-90.0, 90.0), a.x, b.y.clamp(-90.0, 90.0), b.x); s12 } +/// Length of a lon/lat polyline along the WGS 84 ellipsoid's geodesics (metres): each segment is +/// an inverse geodesic solve, so it is exact, unlike [`geodesic_length_m`]'s great circles. Long +/// lines of many vertices cost a solve per segment (a few microseconds each). +pub fn wgs84_length_m(pts: &[Coord]) -> f64 { + pts.windows(2) + .map(|w| match w { + [a, b] => wgs84_distance_m(*a, *b), + _ => 0.0, + }) + .sum() +} + +/// The perimeter (metres) and the unsigned enclosed area (square metres) of a lon/lat ring on the +/// WGS 84 ellipsoid. The ring is closed along a geodesic to its first point, whatever its winding, +/// so a ring of many vertices needs one accumulation, not a solve per vertex. A ring gives the +/// region it traces (so polar rings and globe-wide bands follow their own edges, as the sphere's +/// [`geodesic_ring_area_m2`] does); reversing the ring gives the same magnitude. +pub fn wgs84_ring_perimeter_area_m2(ring: &[Coord]) -> (f64, f64) { + if ring.len() < 3 { + return (0.0, 0.0); + } + // A broken vertex is no place: NaN, never a plausible number (see [`haversine_m`]). + if !ring.iter().all(|c| c.is_finite()) { + return (f64::NAN, f64::NAN); + } + // An edge of a full turn of longitude (-180° to 180°, as a globe-wide band has) joins two + // coordinates that name the same place, which a geodesic cannot tell from no edge at all: such + // a ring is measured on the sphere, as it always was. + let full_turn = ring.iter().zip(ring.iter().cycle().skip(1)).any(|(a, b)| (b.x - a.x).abs() >= 360.0); + if full_turn { + let mut closed = ring.to_vec(); + closed.extend(ring.first().copied()); + return (geodesic_length_m(&closed), geodesic_ring_area_m2(ring)); + } + let mut polygon = PolygonArea::new(&WGS84, Winding::Clockwise); + for c in ring { + polygon.add_point(c.y.clamp(-90.0, 90.0), c.x); + } + // `compute` reduces the signed area to at most half the ellipsoid, so a ring enclosing more + // than half of it comes back as the smaller complement: the sphere's area says which one the + // ring traces. + let (perimeter, area, _) = polygon.compute(true); + let area = area.abs(); + let sphere = geodesic_ring_area_m2(ring); + let complement = WGS84_AREA_M2 - area; + (perimeter, if (sphere - complement).abs() < (sphere - area).abs() { complement } else { area }) +} + +/// The surface area of the WGS 84 ellipsoid (square metres). +const WGS84_AREA_M2: f64 = 510_065_621_724_079.0; + /// The WGS 84 ellipsoid's geodesic solver. static WGS84: LazyLock = LazyLock::new(Geodesic::wgs84); @@ -408,6 +464,89 @@ mod tests { } } + #[test] + fn wgs84_rings_keep_full_turn_edges_and_large_areas() { + let c = Coord::new; + // A globe-wide band: its -180° to 180° edges are parallels, not zero-length hops. + let band = [c(-180.0, -10.0), c(180.0, -10.0), c(180.0, 10.0), c(-180.0, 10.0)]; + let (_, area) = wgs84_ring_perimeter_area_m2(&band); + let sphere = geodesic_ring_area_m2(&band); + assert!(area > 1.0e13 && (area / sphere - 1.0).abs() < 0.01, "{area} {sphere}"); + // A ring around the North Pole's complement: the southern cap of latitude < 60° is more + // than half the Earth; it must not come back as the northern cap. + let big: Vec = (0..=36).map(|i| c(-180.0 + f64::from(i) * 10.0, 60.0)).collect(); + let small = wgs84_ring_perimeter_area_m2(&big).1; + assert!(small < WGS84_AREA_M2 / 2.0, "{small}"); + // A broken latitude is no pole. + assert!(wgs84_length_m(&[c(0.0, 0.0), c(0.0, f64::INFINITY)]).is_nan()); + assert!(wgs84_distance_m(c(0.0, f64::NEG_INFINITY), c(0.0, 0.0)).is_nan()); + } + + #[test] + fn wgs84_measures_match_geographiclib() { + // References from GeographicLib (pyproj.Geod / Polygon.Compute on WGS 84), cross-checked + // against `geo`'s GeodesicArea: a 1° × 1° box at the equator, and the same box stretched + // across the antimeridian (unwrapped, 179° to 181°). + let equator = &[Coord::new(0.0, 0.0), Coord::new(0.0, 1.0), Coord::new(1.0, 1.0), Coord::new(1.0, 0.0)]; + let (perimeter, area) = wgs84_ring_perimeter_area_m2(equator); + assert!((perimeter - 443_770.917).abs() < 0.001, "{perimeter}"); + assert!((area - 12_308_778_361.469).abs() < 0.001, "{area}"); + // The ring traced the other way around gives the same magnitude. + let (perimeter_reversed, area_reversed) = wgs84_ring_perimeter_area_m2(&equator.iter().rev().copied().collect::>()); + assert!((perimeter - perimeter_reversed).abs() < 1e-6); + assert!((area - area_reversed).abs() < 1e-6); + // Unwrapped longitudes past the antimeridian are one ring, not a world-sized complement. + let across = &[Coord::new(179.0, 0.0), Coord::new(181.0, 0.0), Coord::new(181.0, 1.0), Coord::new(179.0, 1.0)]; + let (perimeter_across, area_across) = wgs84_ring_perimeter_area_m2(across); + assert!((perimeter_across - 666_393.055).abs() < 0.001, "{perimeter_across}"); + assert!((area_across - 24_619_443_759.277).abs() < 0.001, "{area_across}"); + // Segment lengths are inverse geodesic solves, unwrapped or folded alike. + assert!((wgs84_length_m(&[Coord::new(0.0, 0.0), Coord::new(3.0, 0.0)]) - 3.0 * 111_319.4908).abs() < 0.01); + assert!((wgs84_distance_m(Coord::new(179.0, 0.0), Coord::new(-179.0, 0.0)) - 222_638.9816).abs() < 0.001); + // A broken coordinate is no length and no area (see `haversine_m`'s guard). + assert!(wgs84_length_m(&[Coord::new(0.0, 0.0), Coord::new(f64::NAN, 0.0)]).is_nan()); + let (perimeter_broken, area_broken) = + wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0), Coord::new(0.0, 1.0), Coord::new(f64::NAN, 0.0)]); + assert!(perimeter_broken.is_nan() && area_broken.is_nan()); + // Degenerate rings keep their zero. + assert_eq!(wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0)]), (0.0, 0.0)); + } + + #[test] + fn a_non_finite_coordinate_is_no_distance_no_area_and_no_place() { + // A NaN vertex once measured as a quarter of the Earth's circumference (π·R = 20 015 114 + // m) through `h.sqrt().min(1.0)`: `f64::min` ignores NaN. It must be no distance at all, so + // sums and panels stay visibly broken and comparisons exclude it (NaN <= r is false). + let nan = f64::NAN; + let origin = Coord::new(0.0, 0.0); + for near in [Coord::new(nan, 0.0), Coord::new(0.0, nan), Coord::new(f64::INFINITY, 0.0)] { + assert!(haversine_m(origin, near).is_nan()); + assert!(geodesic_length_m(&[origin, near, origin]).is_nan()); + } + // A finite edge next to a broken one keeps measuring, through its own points. + assert!((geodesic_length_m(&[origin, Coord::new(1.0, 0.0)]) - 111_195.080).abs() < 0.01); + // Ring area: a NaN absorbing match on the older behaviour, then anything broken. + assert!(geodesic_ring_area_m2(&[origin, Coord::new(1.0, 0.0), Coord::new(1.0, 1.0), Coord::new(nan, 1.0)]).is_nan()); + assert!(geodesic_ring_area_m2(&[origin, Coord::new(1.0, 0.0), Coord::new(1.0, f64::INFINITY)]).is_nan()); + // Heading away from a broken place is no place either. + assert!(geodesic_destination(Coord::new(10.0, 50.0), 90.0, f64::NAN).x.is_nan()); + } + + #[test] + fn non_finite_measures_show_a_dash_not_nan() { + for broken in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + assert_eq!(format_distance(broken), "–"); + assert_eq!(format_area(broken), "–"); + } + assert_eq!(format_lonlat(Coord::new(f64::NAN, f64::NAN)), "–"); + assert_eq!(format_dms(Coord::new(f64::NAN, 50.0)), "–"); + // The finite cases still format as before. + assert_eq!(format_distance(850.0), "850 m"); + assert_eq!(format_distance(12_500.0), "12.5 km"); + assert_eq!(format_area(2_500_000.0), "2.5 km²"); + assert_eq!(format_lonlat(Coord::new(-12.3456, 45.6789)), "12.3456° W 45.6789° N"); + } + #[test] fn azimuthal_equidistant_keeps_distances_from_its_centre() { // Distances and azimuths on WGS 84 from GeographicLib (pyproj.Geod(ellps="WGS84").inv). diff --git a/crates/ui-egui/src/mapview.rs b/crates/ui-egui/src/mapview.rs index 9dd69b4ef..fc9dc36cd 100644 --- a/crates/ui-egui/src/mapview.rs +++ b/crates/ui-egui/src/mapview.rs @@ -8,7 +8,7 @@ use rustgis_engine::analysis::{ network::{radius_boundary_point, radius_distance_m}, selection::MAX_POINTS, }; -use rustgis_engine::geom::{Coord, Extent, geodesic_length_m, geodesic_ring_area_m2}; +use rustgis_engine::geom::{Coord, Extent, geodesic_length_m, geodesic_ring_area_m2, wgs84_length_m, wgs84_ring_perimeter_area_m2}; use rustgis_engine::map::layout::TextAlign; use rustgis_engine::map::projection::Flat; use rustgis_engine::map::{Camera, Globe, Projection}; @@ -2599,16 +2599,24 @@ fn overlays(app: &mut RustGisApp, ui: &mut Ui, p: &egui::Painter, rect: Rect, t: p.circle_stroke(*s, 3.5, Stroke::new(1.5, color)); } if matches!(tool, Tool::Measure | Tool::MeasureArea) && pts.len() >= 2 { - // On the map's body. + // On the map's body: the Earth measures on the WGS 84 ellipsoid, so the panel agrees + // with the geoprocessing tools (#939); another body measures on its sphere and scales. let body = app.session.map().body(); - let geodesic_length_m = |pts: &[Coord]| body.length(geodesic_length_m(pts)); - let geodesic_ring_area_m2 = |ring: &[Coord]| body.area(geodesic_ring_area_m2(ring)); + let measure_length = |pts: &[Coord]| { + if body.is_earth() { wgs84_length_m(pts) } else { body.length(geodesic_length_m(pts)) } + }; + let measure_ring_area = |ring: &[Coord]| { + if body.is_earth() { wgs84_ring_perimeter_area_m2(ring).1 } else { body.area(geodesic_ring_area_m2(ring)) } + }; + let measure_ring_perimeter = |ring: &[Coord]| { + if body.is_earth() { wgs84_ring_perimeter_area_m2(ring).0 } else { body.length(geodesic_length_m(ring)) } + }; let text = if tool == Tool::Measure { - let seg = pts.len().checked_sub(2).and_then(|i| pts.get(i..)).map(geodesic_length_m).unwrap_or(0.0); + let seg = pts.len().checked_sub(2).and_then(|i| pts.get(i..)).map(measure_length).unwrap_or(0.0); crate::i18n::fmt( "Total {total} · Segment {segment}", &[ - ("total", &crate::measure_units::distance(geodesic_length_m(&pts), app.ui.measurement_imperial)), + ("total", &crate::measure_units::distance(measure_length(&pts), app.ui.measurement_imperial)), ("segment", &crate::measure_units::distance(seg, app.ui.measurement_imperial)), ], ) @@ -2620,8 +2628,8 @@ fn overlays(app: &mut RustGisApp, ui: &mut Ui, p: &egui::Painter, rect: Rect, t: crate::i18n::fmt( "Area {area} · Perimeter {perimeter}", &[ - ("area", &crate::measure_units::area(geodesic_ring_area_m2(&ring), app.ui.measurement_imperial)), - ("perimeter", &crate::measure_units::distance(geodesic_length_m(&ring), app.ui.measurement_imperial)), + ("area", &crate::measure_units::area(measure_ring_area(&ring), app.ui.measurement_imperial)), + ("perimeter", &crate::measure_units::distance(measure_ring_perimeter(&ring), app.ui.measurement_imperial)), ], ) }; diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 1be3e3266..17bfae19f 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -132,7 +132,11 @@ use the smaller polar region in either direction; reversing a ring does not sele 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). +differ from the ellipsoid by up to about 0.9 % (a small area near a pole). The same WGS 84 +measure backs the attribute table's **Calculate Geometry** (`edit.calculateGeometry`), the map's +**Measure** tools on the Earth, and every network tool's edge and route lengths, so all of these +read the same number for the same geometry. On the Moon or another body, the Measure tools +measure on that body's sphere. **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 diff --git a/docs/user-guide/data.md b/docs/user-guide/data.md index 0d456df8c..5f5bd4d97 100644 --- a/docs/user-guide/data.md +++ b/docs/user-guide/data.md @@ -92,7 +92,8 @@ For an editable feature layer, double-click a cell to change it. The table's too fields, calculates field expressions and opens **Calculate Geometry** (`edit.calculateGeometry`). The geometry calculator fills a new double field or an existing numeric field with area, perimeter, length or point X/Y, for all or selected features. Point coordinates can be WGS 84 -degrees or a supported projected EPSG system's units. See [Editing](editing.md) for field +degrees or a supported projected EPSG system's units. Areas, perimeters and lengths measure on +the WGS 84 ellipsoid. See [Editing](editing.md) for field calculations and data-design rules. Streamed layers must first be loaded with **Load All Features**. diff --git a/docs/user-guide/editing.md b/docs/user-guide/editing.md index 4a06ed105..d1915b4b7 100644 --- a/docs/user-guide/editing.md +++ b/docs/user-guide/editing.md @@ -10,7 +10,9 @@ to change it. **Data › Data Design** (and **Data › Fields**) adds and delete field from an expression; **Calculate Geometry** (in the attribute table's toolbar) fills a new or existing numeric field with each feature's area, perimeter, length or X/Y coordinate, in the units you choose, for all features or the selected ones. Length and perimeter include feet and miles; -area includes square feet (`ft2`), acres (`ac`) and square miles (`mi2`). Every edit is one undo step. +area includes square feet (`ft2`), acres (`ac`) and square miles (`mi2`). Every edit is one undo +step. Areas, perimeters and lengths measure on the WGS 84 ellipsoid, the same measure the +geoprocessing tools use. Undo (++ctrl+z++) and Redo (++ctrl+y++ or ++ctrl+shift+z++) restore project changes, not just feature edits; use ++cmd++ instead of ++ctrl++ on macOS. While typing in a text field, those From 94b6d92bad805369b51f56723ffd17ceeba82d94 Mon Sep 17 00:00:00 2001 From: subarna-kandel Date: Tue, 27 Oct 2026 02:37:55 +0000 Subject: [PATCH 02/12] fix(geom): address review of the WGS 84 measures wgs84_distance_m returns NaN for non-finite coordinates (an infinite latitude was clamped to a pole); rings with a full-turn longitude edge fall back to the sphere and rings over half the ellipsoid keep their area; the network snap scan ranks segments by haversine and solves only the winner on the ellipsoid. Co-Authored-By: Claude Sonnet 5.5 --- crates/analysis/src/network.rs | 7 +++++-- crates/geom/src/lib.rs | 19 ++----------------- 2 files changed, 7 insertions(+), 19 deletions(-) diff --git a/crates/analysis/src/network.rs b/crates/analysis/src/network.rs index 56d79fa78..0c5149bb6 100644 --- a/crates/analysis/src/network.rs +++ b/crates/analysis/src/network.rs @@ -4,7 +4,7 @@ use std::collections::{BinaryHeap, HashMap}; use rustgis_data::{Dataset, Geometry}; -use rustgis_geom::{Coord, wgs84_distance_m}; +use rustgis_geom::{Coord, haversine_m, wgs84_distance_m}; /// A shortest route: its path, length in metres, cost, and the line rows it uses. #[derive(Clone, Debug, PartialEq)] @@ -96,9 +96,12 @@ pub fn shortest_route( let len2 = dx * dx + dy * dy; let t = if len2 > 0.0 { (((p.x - a.x) * k * dx + (p.y - a.y) * dy) / len2).clamp(0.0, 1.0) } else { 0.0 }; let q = Coord::new(a.x + t * (b.x - a.x), a.y + t * (b.y - a.y)); - Some((i, t, q, wgs84_distance_m(q, p))) + Some((i, t, q, haversine_m(q, p))) }) .min_by(|a, b| a.3.partial_cmp(&b.3).unwrap_or(std::cmp::Ordering::Equal)) + // The scan ranks segments by the cheap great-circle distance (it orders them as the + // ellipsoid does, to well under a snap tolerance); only the winner is solved exactly. + .map(|(i, t, q, _)| (i, t, q, wgs84_distance_m(q, p))) }; // Each end becomes a node on its segment, joined to the segment's two ends. let attach = |p: Coord, nodes: &mut Vec, adj: &mut Vec>| -> Result<(usize, usize, f64), String> { diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index eb1eecfeb..27fb48b05 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -501,8 +501,8 @@ mod tests { assert!((perimeter_across - 666_393.055).abs() < 0.001, "{perimeter_across}"); assert!((area_across - 24_619_443_759.277).abs() < 0.001, "{area_across}"); // Segment lengths are inverse geodesic solves, unwrapped or folded alike. - assert!((wgs84_length_m(&[Coord::new(0.0, 0.0), Coord::new(3.0, 0.0)]) - 3.0 * 111_319.4908).abs() < 0.01); - assert!((wgs84_distance_m(Coord::new(179.0, 0.0), Coord::new(-179.0, 0.0)) - 222_638.9816).abs() < 0.001); + assert!((wgs84_length_m(&[Coord::new(0.0, 0.0), Coord::new(3.0, 0.0)]) - 3.0 * 111_319.490_8).abs() < 0.01); + assert!((wgs84_distance_m(Coord::new(179.0, 0.0), Coord::new(-179.0, 0.0)) - 222_638.981_6).abs() < 0.001); // A broken coordinate is no length and no area (see `haversine_m`'s guard). assert!(wgs84_length_m(&[Coord::new(0.0, 0.0), Coord::new(f64::NAN, 0.0)]).is_nan()); let (perimeter_broken, area_broken) = @@ -532,21 +532,6 @@ mod tests { assert!(geodesic_destination(Coord::new(10.0, 50.0), 90.0, f64::NAN).x.is_nan()); } - #[test] - fn non_finite_measures_show_a_dash_not_nan() { - for broken in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { - assert_eq!(format_distance(broken), "–"); - assert_eq!(format_area(broken), "–"); - } - assert_eq!(format_lonlat(Coord::new(f64::NAN, f64::NAN)), "–"); - assert_eq!(format_dms(Coord::new(f64::NAN, 50.0)), "–"); - // The finite cases still format as before. - assert_eq!(format_distance(850.0), "850 m"); - assert_eq!(format_distance(12_500.0), "12.5 km"); - assert_eq!(format_area(2_500_000.0), "2.5 km²"); - assert_eq!(format_lonlat(Coord::new(-12.3456, 45.6789)), "12.3456° W 45.6789° N"); - } - #[test] fn azimuthal_equidistant_keeps_distances_from_its_centre() { // Distances and azimuths on WGS 84 from GeographicLib (pyproj.Geod(ellps="WGS84").inv). From 5f5fcdc02956fdba8d599afcc92ef272d098177a Mon Sep 17 00:00:00 2001 From: giswqs <5016453+giswqs@users.noreply.github.com> Date: Tue, 27 Oct 2026 03:07:19 +0000 Subject: [PATCH 03/12] fix(geom): test polar and band rings, and measure the ring once The polar-ring test used a 37-vertex ring whose closing edge is a full turn, so the sphere area came back as zero and the "not more than half" check passed trivially. The test now uses a 36-vertex ring with no full-turn edge and asserts it measures the traced southern complement, and the band check is a range against the sphere value. The measure overlay computed the WGS 84 ring area and perimeter in two calls per frame. It now takes both from one call. The user guide names the full-turn exception: a ring with a full-turn edge (a globe-wide band) is measured on the sphere. Co-Authored-By: Claude Haiku 5.5 Applied by the claude-autofix workflow. --- crates/geom/src/lib.rs | 18 ++++++++++-------- crates/ui-egui/src/mapview.rs | 17 ++++++++++------- docs/user-guide/data.md | 3 ++- docs/user-guide/editing.md | 3 ++- 4 files changed, 24 insertions(+), 17 deletions(-) diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 27fb48b05..144273d17 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -467,16 +467,18 @@ mod tests { #[test] fn wgs84_rings_keep_full_turn_edges_and_large_areas() { let c = Coord::new; - // A globe-wide band: its -180° to 180° edges are parallels, not zero-length hops. + // A globe-wide band: its -180° to 180° edges are parallels, not zero-length hops. Measured + // on the sphere: 2πR²(sin 10° − sin −10°) ≈ 8.85e13 m². let band = [c(-180.0, -10.0), c(180.0, -10.0), c(180.0, 10.0), c(-180.0, 10.0)]; let (_, area) = wgs84_ring_perimeter_area_m2(&band); - let sphere = geodesic_ring_area_m2(&band); - assert!(area > 1.0e13 && (area / sphere - 1.0).abs() < 0.01, "{area} {sphere}"); - // A ring around the North Pole's complement: the southern cap of latitude < 60° is more - // than half the Earth; it must not come back as the northern cap. - let big: Vec = (0..=36).map(|i| c(-180.0 + f64::from(i) * 10.0, 60.0)).collect(); - let small = wgs84_ring_perimeter_area_m2(&big).1; - assert!(small < WGS84_AREA_M2 / 2.0, "{small}"); + assert!(area > 8.5e13 && area < 9.2e13, "{area}"); + // A ring of 36 vertices at 60° N, closed by a wrapped 10° edge (no full-turn edge): it traces + // the southern cap of latitude < 60°, more than half the Earth, so it must not come back as + // the northern cap. Its area agrees with the sphere's, as the ring's own edges define it. + let big: Vec = (0..36).map(|i| c(-180.0 + f64::from(i) * 10.0, 60.0)).collect(); + let area = wgs84_ring_perimeter_area_m2(&big).1; + assert!(area > WGS84_AREA_M2 / 2.0, "{area}"); + assert!((area / geodesic_ring_area_m2(&big) - 1.0).abs() < 0.01, "{area}"); // A broken latitude is no pole. assert!(wgs84_length_m(&[c(0.0, 0.0), c(0.0, f64::INFINITY)]).is_nan()); assert!(wgs84_distance_m(c(0.0, f64::NEG_INFINITY), c(0.0, 0.0)).is_nan()); diff --git a/crates/ui-egui/src/mapview.rs b/crates/ui-egui/src/mapview.rs index 6f240bd4a..7b25a5551 100644 --- a/crates/ui-egui/src/mapview.rs +++ b/crates/ui-egui/src/mapview.rs @@ -2651,11 +2651,13 @@ fn overlays(app: &mut RustGisApp, ui: &mut Ui, p: &egui::Painter, rect: Rect, t: let measure_length = |pts: &[Coord]| { if body.is_earth() { wgs84_length_m(pts) } else { body.length(geodesic_length_m(pts)) } }; - let measure_ring_area = |ring: &[Coord]| { - if body.is_earth() { wgs84_ring_perimeter_area_m2(ring).1 } else { body.area(geodesic_ring_area_m2(ring)) } - }; - let measure_ring_perimeter = |ring: &[Coord]| { - if body.is_earth() { wgs84_ring_perimeter_area_m2(ring).0 } else { body.length(geodesic_length_m(ring)) } + // (perimeter, area) of a ring, in one pass. + let measure_ring = |ring: &[Coord]| { + if body.is_earth() { + wgs84_ring_perimeter_area_m2(ring) + } else { + (body.length(geodesic_length_m(ring)), body.area(geodesic_ring_area_m2(ring))) + } }; let text = if tool == Tool::Measure { let seg = pts.len().checked_sub(2).and_then(|i| pts.get(i..)).map(measure_length).unwrap_or(0.0); @@ -2671,11 +2673,12 @@ fn overlays(app: &mut RustGisApp, ui: &mut Ui, p: &egui::Painter, rect: Rect, t: if let Some(f) = ring.first().copied() { ring.push(f); } + let (perimeter, area) = measure_ring(&ring); crate::i18n::fmt( "Area {area} · Perimeter {perimeter}", &[ - ("area", &crate::measure_units::area(measure_ring_area(&ring), app.ui.measurement_imperial)), - ("perimeter", &crate::measure_units::distance(measure_ring_perimeter(&ring), app.ui.measurement_imperial)), + ("area", &crate::measure_units::area(area, app.ui.measurement_imperial)), + ("perimeter", &crate::measure_units::distance(perimeter, app.ui.measurement_imperial)), ], ) }; diff --git a/docs/user-guide/data.md b/docs/user-guide/data.md index f2875cafc..913c53bd5 100644 --- a/docs/user-guide/data.md +++ b/docs/user-guide/data.md @@ -93,7 +93,8 @@ fields, calculates field expressions and opens **Calculate Geometry** (`edit.cal The geometry calculator fills a new double field or an existing numeric field with area, perimeter, length or point X/Y, for all or selected features. Point coordinates can be WGS 84 degrees or a supported projected EPSG system's units. Areas, perimeters and lengths measure on -the WGS 84 ellipsoid. See [Editing](editing.md) for field +the WGS 84 ellipsoid (a ring with a full-turn edge, such as a globe-wide band, on the sphere). +See [Editing](editing.md) for field calculations and data-design rules. Streamed layers must first be loaded with **Load All Features**. diff --git a/docs/user-guide/editing.md b/docs/user-guide/editing.md index b34c82463..22ff53966 100644 --- a/docs/user-guide/editing.md +++ b/docs/user-guide/editing.md @@ -12,7 +12,8 @@ existing numeric field with each feature's area, perimeter, length or X/Y coordi you choose, for all features or the selected ones. Length and perimeter include feet and miles; area includes square feet (`ft2`), acres (`ac`) and square miles (`mi2`). Every edit is one undo step. Areas, perimeters and lengths measure on the WGS 84 ellipsoid, the same measure the -geoprocessing tools use. +geoprocessing tools use. A ring with a full-turn edge (a globe-wide band) is measured on the +sphere instead. Undo (++ctrl+z++) and Redo (++ctrl+y++ or ++ctrl+shift+z++) restore project changes, not just feature edits; use ++cmd++ instead of ++ctrl++ on macOS. While typing in a text field, those From eb349c0cf25ac837dacd3d1e0f999b347a26bb4a Mon Sep 17 00:00:00 2001 From: giswqs <5016453+giswqs@users.noreply.github.com> Date: Tue, 27 Oct 2026 03:25:21 +0000 Subject: [PATCH 04/12] fix(geom): reject broken vertices in short WGS 84 rings A ring of fewer than three vertices returned (0, 0) before the non-finite check ran, so a two-vertex ring with a NaN vertex gave a plausible zero. The check now comes first, so a broken vertex is NaN at any ring length, as the function documents. Also correct the calculator's comment: rings with a full-turn edge are measured on the sphere here, where the geoprocessing tool unwraps them. Add a regression assertion for the short-ring case. Co-Authored-By: Claude Haiku 5.5 Applied by the claude-autofix workflow. --- crates/engine/src/geometry_calc.rs | 3 ++- crates/geom/src/lib.rs | 10 ++++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/crates/engine/src/geometry_calc.rs b/crates/engine/src/geometry_calc.rs index 149bfe419..776a22dcd 100644 --- a/crates/engine/src/geometry_calc.rs +++ b/crates/engine/src/geometry_calc.rs @@ -35,7 +35,8 @@ pub fn units(property: &str) -> &'static [(&'static str, &'static str)] { fn measure(g: &Geometry, property: &str) -> Option { match (g, property) { // On the WGS 84 ellipsoid, like every other measurement named geodesic (#939): the same - // numbers the geoprocessing tools give. + // numbers the geoprocessing tools give, except a ring with a full-turn (-180° to 180°) edge, + // which is measured on the sphere here but unwrapped (so near zero for a band) there. (Geometry::Polygon(polys), "area") => Some( polys .iter() diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 144273d17..b61c37e73 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -270,13 +270,13 @@ pub fn wgs84_length_m(pts: &[Coord]) -> f64 { /// region it traces (so polar rings and globe-wide bands follow their own edges, as the sphere's /// [`geodesic_ring_area_m2`] does); reversing the ring gives the same magnitude. pub fn wgs84_ring_perimeter_area_m2(ring: &[Coord]) -> (f64, f64) { - if ring.len() < 3 { - return (0.0, 0.0); - } // A broken vertex is no place: NaN, never a plausible number (see [`haversine_m`]). if !ring.iter().all(|c| c.is_finite()) { return (f64::NAN, f64::NAN); } + if ring.len() < 3 { + return (0.0, 0.0); + } // An edge of a full turn of longitude (-180° to 180°, as a globe-wide band has) joins two // coordinates that name the same place, which a geodesic cannot tell from no edge at all: such // a ring is measured on the sphere, as it always was. @@ -510,8 +510,10 @@ mod tests { let (perimeter_broken, area_broken) = wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0), Coord::new(0.0, 1.0), Coord::new(f64::NAN, 0.0)]); assert!(perimeter_broken.is_nan() && area_broken.is_nan()); - // Degenerate rings keep their zero. + // Degenerate rings keep their zero, unless a vertex is broken: that is NaN, even in a short ring. assert_eq!(wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0)]), (0.0, 0.0)); + let (perimeter_short, area_short) = wgs84_ring_perimeter_area_m2(&[Coord::new(f64::NAN, f64::NAN), Coord::new(0.0, 0.0)]); + assert!(perimeter_short.is_nan() && area_short.is_nan()); } #[test] From 2a72752e01a21f7c4dde844e869e8eb8f33c63d7 Mon Sep 17 00:00:00 2001 From: giswqs <5016453+giswqs@users.noreply.github.com> Date: Tue, 27 Oct 2026 04:04:30 +0000 Subject: [PATCH 05/12] fix(geom): keep a two-vertex ring's perimeter as 2d The early return for short rings made a two-vertex ring measure zero perimeter, where the mean-sphere version gave twice the segment. It is a segment there and back, so it now measures 2d on the WGS 84 ellipsoid with no area, and a single vertex still measures zero. Also name the full-turn ring exception in the analysis guide, and say that the network snap choice is approximate: the scan ranks segments by great-circle distance. Applied by the claude-autofix workflow. --- crates/analysis/src/network.rs | 5 +++-- crates/geom/src/lib.rs | 12 ++++++++++-- docs/user-guide/analysis.md | 6 ++++-- 3 files changed, 17 insertions(+), 6 deletions(-) diff --git a/crates/analysis/src/network.rs b/crates/analysis/src/network.rs index 0c5149bb6..06a69fe64 100644 --- a/crates/analysis/src/network.rs +++ b/crates/analysis/src/network.rs @@ -99,8 +99,9 @@ pub fn shortest_route( Some((i, t, q, haversine_m(q, p))) }) .min_by(|a, b| a.3.partial_cmp(&b.3).unwrap_or(std::cmp::Ordering::Equal)) - // The scan ranks segments by the cheap great-circle distance (it orders them as the - // ellipsoid does, to well under a snap tolerance); only the winner is solved exactly. + // The scan ranks segments by the cheap great-circle distance, which differs from the + // ellipsoid's by up to about 0.7 %; only the winner is solved exactly, so the segment + // chosen is approximate near the snap tolerance. .map(|(i, t, q, _)| (i, t, q, wgs84_distance_m(q, p))) }; // Each end becomes a node on its segment, joined to the segment's two ends. diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index b61c37e73..18cf492d6 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -274,9 +274,13 @@ pub fn wgs84_ring_perimeter_area_m2(ring: &[Coord]) -> (f64, f64) { if !ring.iter().all(|c| c.is_finite()) { return (f64::NAN, f64::NAN); } - if ring.len() < 3 { + if ring.len() < 2 { return (0.0, 0.0); } + // A two-vertex ring is a segment there and back: twice its length, and no area. + if let [a, b] = ring { + return (2.0 * wgs84_distance_m(*a, *b), 0.0); + } // An edge of a full turn of longitude (-180° to 180°, as a globe-wide band has) joins two // coordinates that name the same place, which a geodesic cannot tell from no edge at all: such // a ring is measured on the sphere, as it always was. @@ -510,8 +514,12 @@ mod tests { let (perimeter_broken, area_broken) = wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0), Coord::new(0.0, 1.0), Coord::new(f64::NAN, 0.0)]); assert!(perimeter_broken.is_nan() && area_broken.is_nan()); - // Degenerate rings keep their zero, unless a vertex is broken: that is NaN, even in a short ring. + // A single vertex keeps its zero; two vertices are a segment there and back (2d, no area). assert_eq!(wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0)]), (0.0, 0.0)); + let (perimeter_two, area_two) = wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0), Coord::new(0.0, 1.0)]); + assert!((perimeter_two - 2.0 * 111_319.490_8).abs() < 0.01, "{perimeter_two}"); + assert_eq!(area_two, 0.0); + // Degenerate rings keep their zero, unless a vertex is broken: that is NaN, even in a short ring. let (perimeter_short, area_short) = wgs84_ring_perimeter_area_m2(&[Coord::new(f64::NAN, f64::NAN), Coord::new(0.0, 0.0)]); assert!(perimeter_short.is_nan() && area_short.is_nan()); } diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index cb772325c..c7500de81 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -137,8 +137,10 @@ Spatial Join's **Largest overlap** still measures area on a sphere of radius 6,3 differ from the ellipsoid by up to about 0.9 % (a small area near a pole). The same WGS 84 measure backs the attribute table's **Calculate Geometry** (`edit.calculateGeometry`), the map's **Measure** tools on the Earth, and every network tool's edge and route lengths, so all of these -read the same number for the same geometry. On the Moon or another body, the Measure tools -measure on that body's sphere. +read the same number for the same geometry, with one exception: a polygon ring with a full-turn +(-180° to 180°) edge, such as a globe-wide band. The attribute table and the Measure tools +measure it on the sphere, while **Calculate Geometry Attributes** unwraps it, so it reads near +zero there. On the Moon or another body, the Measure tools measure on that body's sphere. **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 f02b08ffc38ff36ed2d2a5cb092c9c712995840c Mon Sep 17 00:00:00 2001 From: giswqs <5016453+giswqs@users.noreply.github.com> Date: Tue, 27 Oct 2026 11:55:10 +0000 Subject: [PATCH 06/12] fix(geom): correct the perimeter test and document polar rings The two-vertex perimeter test used a segment from (0, 1) to (0, 0), which runs along a meridian (about 110.6 km per degree), not the longitude degree its 111,319.49 m reference assumes. Its endpoint now runs (1, 0) to (0, 0), so the assertion holds. Polar rings are measured differently by the attribute table (the region the ring traces) and by Calculate Geometry Attributes (the smaller polar cap). The docs and the calculator's comment now list this as an exception, next to full-turn rings. Co-Authored-By: Claude Haiku 5.5 Applied by the claude-autofix workflow. --- crates/engine/src/geometry_calc.rs | 5 +++-- crates/geom/src/lib.rs | 2 +- docs/user-guide/analysis.md | 11 +++++++---- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/crates/engine/src/geometry_calc.rs b/crates/engine/src/geometry_calc.rs index 776a22dcd..cb7841f93 100644 --- a/crates/engine/src/geometry_calc.rs +++ b/crates/engine/src/geometry_calc.rs @@ -35,8 +35,9 @@ pub fn units(property: &str) -> &'static [(&'static str, &'static str)] { fn measure(g: &Geometry, property: &str) -> Option { match (g, property) { // On the WGS 84 ellipsoid, like every other measurement named geodesic (#939): the same - // numbers the geoprocessing tools give, except a ring with a full-turn (-180° to 180°) edge, - // which is measured on the sphere here but unwrapped (so near zero for a band) there. + // numbers the geoprocessing tools give, except two kinds of ring. One with a full-turn + // (-180° to 180°) edge is measured on the sphere here but unwrapped (so near zero for a band) + // there. A ring winding round a pole is the region it traces here, but the smaller polar cap there. (Geometry::Polygon(polys), "area") => Some( polys .iter() diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 18cf492d6..22e3bfbe9 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -516,7 +516,7 @@ mod tests { assert!(perimeter_broken.is_nan() && area_broken.is_nan()); // A single vertex keeps its zero; two vertices are a segment there and back (2d, no area). assert_eq!(wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0)]), (0.0, 0.0)); - let (perimeter_two, area_two) = wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0), Coord::new(0.0, 1.0)]); + let (perimeter_two, area_two) = wgs84_ring_perimeter_area_m2(&[Coord::new(0.0, 0.0), Coord::new(1.0, 0.0)]); assert!((perimeter_two - 2.0 * 111_319.490_8).abs() < 0.01, "{perimeter_two}"); assert_eq!(area_two, 0.0); // Degenerate rings keep their zero, unless a vertex is broken: that is NaN, even in a short ring. diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 30e57d433..771353e3e 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -139,10 +139,13 @@ Spatial Join's **Largest overlap** still measures area on a sphere of radius 6,3 differ from the ellipsoid by up to about 0.9 % (a small area near a pole). The same WGS 84 measure backs the attribute table's **Calculate Geometry** (`edit.calculateGeometry`), the map's **Measure** tools on the Earth, and every network tool's edge and route lengths, so all of these -read the same number for the same geometry, with one exception: a polygon ring with a full-turn -(-180° to 180°) edge, such as a globe-wide band. The attribute table and the Measure tools -measure it on the sphere, while **Calculate Geometry Attributes** unwraps it, so it reads near -zero there. On the Moon or another body, the Measure tools measure on that body's sphere. +read the same number for the same geometry, with two exceptions. A polygon ring with a full-turn +(-180° to 180°) edge, such as a globe-wide band, is measured on the sphere by the attribute table +and the Measure tools, while **Calculate Geometry Attributes** unwraps it, so it reads near zero +there. A ring winding around a pole is measured as the region it traces by the attribute table +(a ring that goes all the way round at 60° N gives the larger region), while **Calculate Geometry +Attributes** gives the smaller polar cap. On the Moon or another body, the Measure tools measure on +that body's sphere. **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 45a901a66f6c7c8f27e176396d6ce5b55f24e49f Mon Sep 17 00:00:00 2001 From: giswqs <5016453+giswqs@users.noreply.github.com> Date: Tue, 27 Oct 2026 12:12:17 +0000 Subject: [PATCH 07/12] docs(analysis): name the measures not yet on WGS 84 The analysis guide says the WGS 84 measure backs every length and area the user sees, but the Radius tools still use a local Mercator distance and raster.measure (Mensuration Tools) still measures on the sphere. Name both so the guide doesn't overstate the agreement. Also reword the wgs84_ring_perimeter_area_m2 comment: a ring costs one geodesic solve per vertex either way, and the saving is a single pass that yields the perimeter and area together, not fewer solves. Co-Authored-By: Claude Haiku 5.5 Applied by the claude-autofix workflow. --- crates/geom/src/lib.rs | 2 +- docs/user-guide/analysis.md | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index 22e3bfbe9..f1ec80d0e 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -266,7 +266,7 @@ pub fn wgs84_length_m(pts: &[Coord]) -> f64 { /// The perimeter (metres) and the unsigned enclosed area (square metres) of a lon/lat ring on the /// WGS 84 ellipsoid. The ring is closed along a geodesic to its first point, whatever its winding, -/// so a ring of many vertices needs one accumulation, not a solve per vertex. A ring gives the +/// so the perimeter and area come from one pass over the ring, not two. A ring gives the /// region it traces (so polar rings and globe-wide bands follow their own edges, as the sphere's /// [`geodesic_ring_area_m2`] does); reversing the ring gives the same magnitude. pub fn wgs84_ring_perimeter_area_m2(ring: &[Coord]) -> (f64, f64) { diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 771353e3e..579627172 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -145,7 +145,8 @@ and the Measure tools, while **Calculate Geometry Attributes** unwraps it, so it there. A ring winding around a pole is measured as the region it traces by the attribute table (a ring that goes all the way round at 60° N gives the larger region), while **Calculate Geometry Attributes** gives the smaller polar cap. On the Moon or another body, the Measure tools measure on -that body's sphere. +that body's sphere. The Radius tools still measure a local Mercator distance, and **Mensuration +Tools** (`raster.measure`) still measure on the sphere, so those can differ from the numbers above. **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 7131b5462088f0a2b50ca7113c69143ec458a20c Mon Sep 17 00:00:00 2001 From: giswqs <5016453+giswqs@users.noreply.github.com> Date: Tue, 27 Oct 2026 12:55:12 +0000 Subject: [PATCH 08/12] fix(analysis): reject NaN route ends; document polar measures A route end with a NaN coordinate gave a NaN distance, and the snap check (`d > snap_m`) is false for NaN, so the end was attached to an arbitrary segment. The check now fails on NaN too, and a test covers it. The attribute table's Calculate Geometry description now names the polar-ring case, where the geoprocessing tool gives the smaller polar cap. The Geoprocessing and Data pages link to the exceptions. Co-Authored-By: Claude Haiku 5.5 Applied by the claude-autofix workflow. --- crates/analysis/src/network.rs | 4 +++- crates/engine/tests/geometry_calculator.rs | 2 +- docs/user-guide/data.md | 3 ++- docs/user-guide/editing.md | 5 +++-- 4 files changed, 9 insertions(+), 5 deletions(-) diff --git a/crates/analysis/src/network.rs b/crates/analysis/src/network.rs index 06a69fe64..456adefcc 100644 --- a/crates/analysis/src/network.rs +++ b/crates/analysis/src/network.rs @@ -107,7 +107,8 @@ pub fn shortest_route( // Each end becomes a node on its segment, joined to the segment's two ends. let attach = |p: Coord, nodes: &mut Vec, adj: &mut Vec>| -> Result<(usize, usize, f64), String> { let (si, t, q, d) = nearest(p, nodes).ok_or("the layer has no lines")?; - if d > snap_m { + // A NaN distance (a broken end) fails the check too, rather than snapping to a random segment. + if d.is_nan() || d > snap_m { return Err(format!("a route end is more than {} from the network", rustgis_geom::format_distance(snap_m))); } let &(a, b, weight, row) = segs.get(si).ok_or("the layer has no lines")?; @@ -452,6 +453,7 @@ mod tests { assert!((r.length_m - 3336.0).abs() < 20.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"); + assert!(shortest_route(&ds, &[0], Coord::new(f64::NAN, 0.0), Coord::new(0.06, 0.0), None, 50.0).is_err(), "a NaN end"); } fn offset(at: Coord, x: f64, y: f64) -> Coord { diff --git a/crates/engine/tests/geometry_calculator.rs b/crates/engine/tests/geometry_calculator.rs index 8226ea5b0..a5daeef2a 100644 --- a/crates/engine/tests/geometry_calculator.rs +++ b/crates/engine/tests/geometry_calculator.rs @@ -52,7 +52,7 @@ fn selected_units_overwrite_and_undo() { #[test] fn areas_are_on_the_wgs84_ellipsoid() { // An equator 1° × 1° box measured by GeographicLib itself (pyproj.Geod / Polygon.Compute on - // WGS 84), so the attribute table and the geoprocessing tool agree with each other (#939). + // 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"})) diff --git a/docs/user-guide/data.md b/docs/user-guide/data.md index 913c53bd5..c11d7cb1a 100644 --- a/docs/user-guide/data.md +++ b/docs/user-guide/data.md @@ -93,7 +93,8 @@ fields, calculates field expressions and opens **Calculate Geometry** (`edit.cal The geometry calculator fills a new double field or an existing numeric field with area, perimeter, length or point X/Y, for all or selected features. Point coordinates can be WGS 84 degrees or a supported projected EPSG system's units. Areas, perimeters and lengths measure on -the WGS 84 ellipsoid (a ring with a full-turn edge, such as a globe-wide band, on the sphere). +the WGS 84 ellipsoid, with the exceptions for full-turn and polar rings listed under +[Geoprocessing](analysis.md#geoprocessing-tools). See [Editing](editing.md) for field calculations and data-design rules. Streamed layers must first be loaded with **Load All Features**. diff --git a/docs/user-guide/editing.md b/docs/user-guide/editing.md index 22ff53966..57df3e6b3 100644 --- a/docs/user-guide/editing.md +++ b/docs/user-guide/editing.md @@ -12,8 +12,9 @@ existing numeric field with each feature's area, perimeter, length or X/Y coordi you choose, for all features or the selected ones. Length and perimeter include feet and miles; area includes square feet (`ft2`), acres (`ac`) and square miles (`mi2`). Every edit is one undo step. Areas, perimeters and lengths measure on the WGS 84 ellipsoid, the same measure the -geoprocessing tools use. A ring with a full-turn edge (a globe-wide band) is measured on the -sphere instead. +geoprocessing tools use, with two exceptions: a ring with a full-turn edge (a globe-wide band) +is measured on the sphere instead, and a ring winding around a pole is measured as the region it +traces, where [Geoprocessing](analysis.md#geoprocessing-tools) gives the smaller polar cap. Undo (++ctrl+z++) and Redo (++ctrl+y++ or ++ctrl+shift+z++) restore project changes, not just feature edits; use ++cmd++ instead of ++ctrl++ on macOS. While typing in a text field, those From 669841820e9a0a4155d80dbb4acd5212fead871b Mon Sep 17 00:00:00 2001 From: giswqs <5016453+giswqs@users.noreply.github.com> Date: Tue, 27 Oct 2026 13:31:32 +0000 Subject: [PATCH 09/12] fix(analysis): keep broken segments out of the route snap The route snap ranks segments by distance with min_by, which treats NaN as equal to every candidate. A broken segment (a NaN vertex) listed first kept its NaN distance as the minimum, so every route end failed to snap even when a valid line lay within the tolerance. Rank a NaN distance as infinity instead. A NaN route end still fails the snap check with the distance error, not "the layer has no lines". Add a regression test with a broken segment listed before a valid one. Co-Authored-By: Claude Haiku 5.5 Applied by the claude-autofix workflow. --- crates/analysis/src/network.rs | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/crates/analysis/src/network.rs b/crates/analysis/src/network.rs index 456adefcc..d982a1bda 100644 --- a/crates/analysis/src/network.rs +++ b/crates/analysis/src/network.rs @@ -96,7 +96,10 @@ pub fn shortest_route( let len2 = dx * dx + dy * dy; let t = if len2 > 0.0 { (((p.x - a.x) * k * dx + (p.y - a.y) * dy) / len2).clamp(0.0, 1.0) } else { 0.0 }; 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))) + let d = haversine_m(q, p); + // A NaN distance (a broken segment) ranks last: `min_by` treats NaN as equal to + // everything, so it could keep one as the minimum and hide the valid segments. + Some((i, t, q, if d.is_nan() { f64::INFINITY } else { d })) }) .min_by(|a, b| a.3.partial_cmp(&b.3).unwrap_or(std::cmp::Ordering::Equal)) // The scan ranks segments by the cheap great-circle distance, which differs from the @@ -456,6 +459,16 @@ mod tests { assert!(shortest_route(&ds, &[0], Coord::new(f64::NAN, 0.0), Coord::new(0.06, 0.0), None, 50.0).is_err(), "a NaN end"); } + #[test] + fn a_broken_segment_listed_first_does_not_hide_a_valid_one() { + // The NaN vertex makes the first segment's distance NaN, which the snap must not keep. + 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 - 3336.0).abs() < 20.0, "{}", r.length_m); + } + fn offset(at: Coord, x: f64, y: f64) -> Coord { let m = rustgis_geom::lonlat_to_merc(at); let scale = radius_scale(at); From 56d579cfef5da10d88cbad80876ba5f18006f653 Mon Sep 17 00:00:00 2001 From: subarna-kandel Date: Wed, 28 Oct 2026 06:22:49 +0000 Subject: [PATCH 10/12] test(analysis,engine): pin the ellipsoid lengths; docs: Spatial Join is still spherical Co-Authored-By: Claude Sonnet 5.5 --- crates/analysis/src/network.rs | 4 ++-- crates/engine/tests/geometry_calculator.rs | 2 +- docs/user-guide/editing.md | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/crates/analysis/src/network.rs b/crates/analysis/src/network.rs index d982a1bda..de409c36a 100644 --- a/crates/analysis/src/network.rs +++ b/crates/analysis/src/network.rs @@ -453,7 +453,7 @@ 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"); assert!(shortest_route(&ds, &[0], Coord::new(f64::NAN, 0.0), Coord::new(0.06, 0.0), None, 50.0).is_err(), "a NaN end"); @@ -466,7 +466,7 @@ mod tests { 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 - 3336.0).abs() < 20.0, "{}", r.length_m); + assert!((r.length_m - 3339.6).abs() < 1.0, "{}", r.length_m); } fn offset(at: Coord, x: f64, y: f64) -> Coord { diff --git a/crates/engine/tests/geometry_calculator.rs b/crates/engine/tests/geometry_calculator.rs index a5daeef2a..de2d6563b 100644 --- a/crates/engine/tests/geometry_calculator.rs +++ b/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)); diff --git a/docs/user-guide/editing.md b/docs/user-guide/editing.md index 57df3e6b3..2a633c817 100644 --- a/docs/user-guide/editing.md +++ b/docs/user-guide/editing.md @@ -11,8 +11,8 @@ field from an expression; **Calculate Geometry** (in the attribute table's toolb existing numeric field with each feature's area, perimeter, length or X/Y coordinate, in the units you choose, for all features or the selected ones. Length and perimeter include feet and miles; area includes square feet (`ft2`), acres (`ac`) and square miles (`mi2`). Every edit is one undo -step. Areas, perimeters and lengths measure on the WGS 84 ellipsoid, the same measure the -geoprocessing tools use, with two exceptions: a ring with a full-turn edge (a globe-wide band) +step. Areas, perimeters and lengths measure on the WGS 84 ellipsoid, the same measure as the network tools and the Measure panel (Spatial Join still +measures on the sphere), with two exceptions: a ring with a full-turn edge (a globe-wide band) is measured on the sphere instead, and a ring winding around a pole is measured as the region it traces, where [Geoprocessing](analysis.md#geoprocessing-tools) gives the smaller polar cap. From 6d2d8ee0f9dbae310a8a46b387810cf57f59ed72 Mon Sep 17 00:00:00 2001 From: mavrick Date: Sat, 31 Oct 2026 20:22:28 -0500 Subject: [PATCH 11/12] fix(network): report omitted segments and reject invalid endpoints --- crates/analysis/src/network.rs | 24 ++++++++++++++++++------ crates/engine/src/cmd/workflows.rs | 8 ++++++++ crates/engine/tests/workflows.rs | 26 ++++++++++++++++++++++++++ 3 files changed, 52 insertions(+), 6 deletions(-) diff --git a/crates/analysis/src/network.rs b/crates/analysis/src/network.rs index 1a7a70334..5038688a0 100644 --- a/crates/analysis/src/network.rs +++ b/crates/analysis/src/network.rs @@ -13,6 +13,8 @@ pub struct Route { pub length_m: f64, pub cost: f64, pub rows: Vec, + /// Segments omitted because an endpoint has non-finite coordinates or an invalid latitude. + pub skipped_segments: usize, } #[derive(PartialEq)] @@ -46,6 +48,10 @@ pub fn shortest_route( cost_field: Option, snap_m: f64, ) -> Result { + if [from, to].iter().any(|c| !c.is_finite() || c.y.abs() > 90.0) { + return Err("a route end has invalid coordinates".into()); + } + let mut skipped_segments = 0; let mut ids: HashMap<(i64, i64), usize> = HashMap::new(); let mut nodes: Vec = Vec::new(); let mut adj: Vec> = Vec::new(); // (to, cost, length, row) @@ -71,7 +77,8 @@ pub fn shortest_route( 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() { + 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)); @@ -90,6 +97,9 @@ pub fn shortest_route( } } } + if segs.is_empty() && skipped_segments > 0 { + return Err(format!("no usable network remains: {skipped_segments} segments have invalid coordinates")); + } // 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); @@ -121,8 +131,7 @@ pub fn shortest_route( // Each end becomes a node on its segment, joined to the segment's two ends. let attach = |p: Coord, nodes: &mut Vec, adj: &mut Vec>| -> Result<(usize, usize, f64), String> { let (si, t, q, d) = nearest(p, nodes).ok_or("the layer has no lines")?; - // A NaN distance (a broken end) fails the check too, rather than snapping to a random segment. - if d.is_nan() || d > snap_m { + if d > snap_m { return Err(format!("a route end is more than {} from the network", rustgis_geom::format_distance(snap_m))); } let &(a, b, weight, row) = segs.get(si).ok_or("the layer has no lines")?; @@ -209,7 +218,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. @@ -467,12 +476,15 @@ mod tests { 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"); - assert!(shortest_route(&ds, &[0], Coord::new(f64::NAN, 0.0), Coord::new(0.06, 0.0), None, 50.0).is_err(), "a NaN end"); + 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}"); } #[test] fn a_broken_segment_listed_first_does_not_hide_a_valid_one() { - // The NaN vertex makes the first segment's distance NaN, which the snap must not keep. + // 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![]); diff --git a/crates/engine/src/cmd/workflows.rs b/crates/engine/src/cmd/workflows.rs index 84832a97d..d3f48986e 100644 --- a/crates/engine/src/cmd/workflows.rs +++ b/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()})) } diff --git a/crates/engine/tests/workflows.rs b/crates/engine/tests/workflows.rs index 22863b21b..c5ff5a238 100644 --- a/crates/engine/tests/workflows.rs +++ b/crates/engine/tests/workflows.rs @@ -348,3 +348,29 @@ 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![]); + 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:?}"); + } + } +} From 5f5b3a9b32c4b5e9032548ef20b841f1b2c55c23 Mon Sep 17 00:00:00 2001 From: mavrick Date: Sat, 31 Oct 2026 22:44:36 -0500 Subject: [PATCH 12/12] fix(network): retain omitted-segment context on route failures --- crates/analysis/src/network.rs | 29 ++++++++++++++++++++++------- crates/engine/tests/workflows.rs | 18 +++++++++++++++++- 2 files changed, 39 insertions(+), 8 deletions(-) diff --git a/crates/analysis/src/network.rs b/crates/analysis/src/network.rs index 5038688a0..a5c8f51ac 100644 --- a/crates/analysis/src/network.rs +++ b/crates/analysis/src/network.rs @@ -48,9 +48,20 @@ pub fn shortest_route( cost_field: Option, snap_m: f64, ) -> Result { + 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 = Vec::new(); @@ -92,13 +103,13 @@ 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(format!("no usable network remains: {skipped_segments} segments have invalid coordinates")); + 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)> { @@ -130,11 +141,14 @@ pub fn shortest_route( }; // Each end becomes a node on its segment, joined to the segment's two ends. let attach = |p: Coord, nodes: &mut Vec, adj: &mut Vec>| -> 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)); } @@ -160,7 +174,7 @@ 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 = wgs84_distance_m(a, b); if let Some(v) = adj.get_mut(s) { @@ -200,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; @@ -480,6 +494,7 @@ mod tests { 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] diff --git a/crates/engine/tests/workflows.rs b/crates/engine/tests/workflows.rs index c5ff5a238..8d71711d1 100644 --- a/crates/engine/tests/workflows.rs +++ b/crates/engine/tests/workflows.rs @@ -356,7 +356,13 @@ fn network_routes_report_invalid_segments_but_clean_networks_are_quiet() { 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![]); + 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![]); @@ -372,5 +378,15 @@ fn network_routes_report_invalid_segments_but_clean_networks_are_quiet() { } 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); + } + } } }