fix(analysis): measure radius searches on the ground near a pole, and meet at the poles (#1033)

crates/analysis/src/convert.rs

@@ -1,18 +1,26 @@
-//! Conversions between RustGIS geometry and `geo` types (in Web Mercator metres).
+//! Conversions between RustGIS geometry and `geo` types, in the analysis plane
+//! ([`rustgis_geom::lonlat_to_plane`]): Web Mercator meters within ±85.05°, carried on to the poles
+//! beyond, so polar geometry isn't clamped onto 85.05° (#1033).
 
 use geo::{Coord as GC, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon};
 use rustgis_data::Geometry;
-use rustgis_geom::{Coord, lonlat_to_merc, merc_to_lonlat};
+use rustgis_geom::{Coord, lonlat_to_plane, plane_to_lonlat};
 
 /// Projects without folding the longitude: a ring kept unwrapped across the antimeridian (as a
 /// geodesic buffer is) stays one piece, rather than folding its far vertex across the world.
 fn gc(c: Coord) -> GC<f64> {
-    let m = lonlat_to_merc(c);
+    let m = lonlat_to_plane(c);
     GC { x: m.x, y: m.y }
 }
 
 fn back(c: &GC<f64>) -> Coord {
-    merc_to_lonlat(Coord::new(c.x, c.y))
+    plane_to_lonlat(Coord::new(c.x, c.y))
+}
+
+/// The poles something spanning latitudes `south` to `north` reaches (to within
+/// [`rustgis_geom::POLE_TOLERANCE`]): the North Pole, then the South Pole.
+pub(crate) fn reaches_poles(south: f64, north: f64) -> (bool, bool) {
+    (north >= 90.0 - rustgis_geom::POLE_TOLERANCE, south <= -90.0 + rustgis_geom::POLE_TOLERANCE)
 }
 
 pub fn ring(r: &[Coord]) -> LineString<f64> {
@@ -96,7 +104,7 @@ mod tests {
         let xs: Vec<f64> = mp.0.iter().flat_map(|p| p.exterior().0.iter().map(|c| c.x)).collect();
         let width = xs.iter().copied().fold(f64::NEG_INFINITY, f64::max) - xs.iter().copied().fold(f64::INFINITY, f64::min);
         assert!(width > 20_000.0 && width < 25_000.0, "{width}");
-        let back: Vec<Coord> = mp.0.iter().flat_map(|p| p.exterior().0.iter().map(|c| merc_to_lonlat(Coord::new(c.x, c.y)))).collect();
+        let back: Vec<Coord> = mp.0.iter().flat_map(|p| p.exterior().0.iter().map(|c| plane_to_lonlat(Coord::new(c.x, c.y)))).collect();
         for (a, b) in ring.iter().zip(back.iter()) {
             assert!((a.x - b.x).abs() < 1e-9 && (a.y - b.y).abs() < 1e-9, "{a:?} → {b:?}");
         }

crates/analysis/src/lib.rs

@@ -1,9 +1,13 @@
 //! L2: geoprocessing tools.
 //!
 //! Each tool is a [`ToolSpec`]: an id, a toolbox, parameters for the Tools pane, and a
-//! pure `run` function over datasets. Inputs refer to layers by name through [`Inputs`]. Work is
-//! done in Web Mercator metres, except buffers, which are geodesic: drawn in a projection centred
-//! on each feature (see [`geodesic`]) so a 100 km buffer is 100 km on the ground.
+//! pure `run` function over datasets. Inputs refer to layers by name through [`Inputs`].
+//! The geometry tools and topology checks work in the analysis plane
+//! ([`rustgis_geom::lonlat_to_plane`]: Web Mercator meters within ±85.05°, continued beyond to
+//! the poles), radius searches in Web Mercator meters (on the ground when the center or the
+//! feature's box reaches beyond 85.05°, see [`network::within`]), and buffers are geodesic: drawn
+//! in a projection centered on each feature (see [`geodesic`]) so a 100 km buffer is 100 km on the
+//! ground.
 //!
 //! Tools from external suites (the GeoLibre/Whitebox tools in `rustgis-geolibre`, and plugins
 //! in `rustgis-plugins`) join the built-in ones through [`register_external`]: they're listed by [`tools`], found by

crates/analysis/tests/reference.rs

@@ -104,7 +104,7 @@ fn polygons(gs: &[Geometry]) -> MultiPolygon<f64> {
     MultiPolygon::new(gs.iter().filter_map(to_multipolygon).flat_map(|m| m.0).collect())
 }
 
-/// The lines among `gs`, together, in Web Mercator metres.
+/// The lines among `gs`, together, in the analysis plane.
 fn lines(gs: &[Geometry]) -> MultiLineString<f64> {
     MultiLineString::new(gs.iter().filter_map(to_multiline).flat_map(|m| m.0).collect())
 }

crates/analysis/tests/tabulate_intersection.rs

@@ -255,10 +255,10 @@ fn overlapping_zones_still_say_what_lies_outside() {
     assert_eq!(got.iter().map(|(_, n)| n[2]).collect::<Vec<_>>(), [50.0, 50.0]);
 }
 
-/// A class feature reaching beyond the map's latitude limit, and features lacking values in two sum
-/// fields: one notice for each kind.
+/// Features lacking values in two sum fields: one notice for each field, and none for a zone and a
+/// class feature reaching past 85.05°, which are measured there like anywhere else (#1033).
 #[test]
-fn one_notice_per_kind_for_missing_values_and_the_latitude_limit() {
+fn one_notice_per_field_for_missing_values() {
     let polar = |x0: f64, x1: f64, y0: f64, y1: f64| {
         let ring = [(x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0)].map(|(x, y)| Coord::new(x, y)).to_vec();
         Some(Geometry::Polygon(vec![vec![ring]]))
@@ -285,12 +285,55 @@ fn one_notice_per_kind_for_missing_values_and_the_latitude_limit() {
             "Tabulate Intersection: 1 row",
             "1 zone overlaps no class feature, so it has no rows",
             "class features with no value in a sum field are left out of its shares: 2 class features in 'pop', 1 class feature in 'income'",
-            "1 zone feature and 1 class feature reach beyond 85.05° latitude, where the map ends: parts there are taken as lying along it, with no area",
+            "1 zone feature and 1 class feature have an edge beyond 85.05° running east or west: the analysis plane takes it straight in longitude and latitude beyond 85.05°, and their percentages and shares follow the plane, while area_km2 takes the shortest way on the ground, which near a pole can differ widely",
             "1 class feature lies partly or wholly outside the zones: its value outside isn't shared out",
         ]
     );
 }
 
+/// Polygons reaching past 85.05°, where the overlay goes on to the poles (#1033): a class feature
+/// from 84°N to 89°N shares out 91.42% of its value to a zone from 80°N to 88°N (on the sphere,
+/// the part from 84° to 88° over the whole), not all of it as when its part beyond 85.05° had no
+/// area; and one wholly beyond, from 86°N to 90°N, has its row in a zone from 80°N to the pole,
+/// 16.03% of it, rather than none.
+#[test]
+fn polygons_past_85_degrees_are_measured_to_the_pole() {
+    let polar = |x0: f64, x1: f64, y0: f64, y1: f64| {
+        let ring = [(x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0)].map(|(x, y)| Coord::new(x, y)).to_vec();
+        Some(Geometry::Polygon(vec![vec![ring]]))
+    };
+    let mut z = Dataset::new("zones", GeometryType::Polygon, vec![Field::new("zone", FieldType::Text)]);
+    z.push_next(polar(0.0, 10.0, 80.0, 88.0), vec![Value::Text("A".into())]);
+    z.push_next(polar(20.0, 30.0, 80.0, 90.0), vec![Value::Text("B".into())]);
+    let mut c = Dataset::new("classes", GeometryType::Polygon, vec![Field::new("pop", FieldType::Double)]);
+    c.push_next(polar(0.0, 10.0, 84.0, 89.0), vec![Value::Float(1000.0)]);
+    c.push_next(polar(20.0, 30.0, 86.0, 90.0), vec![Value::Float(500.0)]);
+    let l = layers(z, c);
+    let params = json!({"zones": "zones", "zoneFields": "zone", "classes": "classes", "sumFields": ["pop"]});
+    let res = rustgis_analysis::run("tabulateIntersection", &l, &params).unwrap();
+    // Each has an edge beyond 85.05° running east or west, which the map and the ground take
+    // differently: the tool says so.
+    let said = "2 zone features and 2 class features have an edge beyond 85.05° running east or west";
+    assert!(res.messages.iter().any(|m| m.starts_with(said)), "{:?}", res.messages);
+    let got = rows(&table(res));
+    // zone, then area_km2, percentage, pop.
+    assert_eq!(got.len(), 2, "{got:?}");
+    let (a, b) = (&got[0], &got[1]);
+    assert_eq!((a.0[0].as_str(), b.0[0].as_str()), ("A", "B"));
+    // The parts' areas on WGS 84, as polygons of geodesic edges (GeographicLib 2.1 gives 34,625.271
+    // and 17,323.559 km²), not 0 for the stretch beyond 85.05°.
+    assert!((a.1[0] - 34_625.271).abs() < 0.01 && (b.1[0] - 17_323.559).abs() < 0.01, "{a:?} {b:?}");
+    assert!((a.1[1] - 33.388).abs() < 0.01 && (a.1[2] - 914.216).abs() < 0.1, "{a:?}");
+    assert!((b.1[1] - 16.034).abs() < 0.01 && (b.1[2] - 500.0).abs() < 1e-6, "{b:?}");
+    // Within ±85.05°, or beyond it only along meridians, nothing of the kind.
+    let mut z = Dataset::new("zones", GeometryType::Polygon, vec![Field::new("zone", FieldType::Text)]);
+    z.push_next(polar(0.0, 10.0, 70.0, 85.0), vec![Value::Text("A".into())]);
+    let mut c = Dataset::new("classes", GeometryType::Polygon, vec![Field::new("pop", FieldType::Double)]);
+    c.push_next(polar(0.0, 10.0, 75.0, 85.0), vec![Value::Float(1000.0)]);
+    let res = rustgis_analysis::run("tabulateIntersection", &layers(z, c), &params).unwrap();
+    assert!(res.messages.iter().all(|m| !m.contains("beyond 85.05°")), "{:?}", res.messages);
+}
+
 /// Features of about 10 m in zones of half a degree: the overlay rounds their parts to a grid
 /// scaled to the zones, so a feature wholly inside a zone still counts whole, and one across the
 /// border between two zones shares its value out in full within that rounding, with no notice.

crates/engine/src/cmd/design.rs

@@ -1289,8 +1289,11 @@ fn topology_validate(s: &mut Session, p: &Json) -> CmdResult {
     use rustgis_analysis::topology as t;
     let tol = f(p, "toleranceM").unwrap_or(1.0).clamp(0.0, 10_000.0);
     let mut errors = Vec::new();
-    // Rules whose check stopped early (an error or work limit): their counts are partial.
+    // Rules whose check stopped early (an error or work limit), or measured some part more coarsely
+    // than it's meant to (a line winding round a pole many times): their counts may be off.
     let mut incomplete = Vec::new();
+    // Polygons a check had to leave out, by rule.
+    let mut skipped = Vec::new();
     for rule in s.map().editing.topology.clone() {
         let Some(l) = s.map().layer_by_name(rule.layer()) else { continue };
         let rows: Vec<usize> = l.visible_rows().iter().enumerate().filter(|(_, v)| **v).map(|(i, _)| i).collect();
@@ -1307,6 +1310,14 @@ fn topology_validate(s: &mut Session, p: &Json) -> CmdResult {
         if !found.complete {
             incomplete.push(rule.label().to_string());
         }
+        if found.skipped > 0 {
+            // The features left out are the covering layer's.
+            let of = match &rule {
+                TopologyRule::MustBeCoveredBy { by, .. } => by.as_str(),
+                _ => l.name.as_str(),
+            };
+            skipped.push(format!("{} on '{}': {} feature(s) of '{}'", rule.label(), l.name, found.skipped, of));
+        }
         errors.extend(found.errors.into_iter().map(|e| crate::TopologyIssue {
             rule: rule.label().to_string(),
             layer: l.name.clone(),
@@ -1315,6 +1326,18 @@ fn topology_validate(s: &mut Session, p: &Json) -> CmdResult {
             lat: e.at.y,
         }));
     }
+    if !incomplete.is_empty() {
+        s.warn(format!(
+            "Checks that stopped early (an error or work limit) or measured a line round a pole coarsely, so their counts may be off: {}",
+            incomplete.join(", ")
+        ));
+    }
+    if !skipped.is_empty() {
+        s.warn(format!(
+            "Polygons left out of the check at a pole (a coordinate missing or infinite, or a million turns around), so its result there may be wrong: {}",
+            skipped.join("; ")
+        ));
+    }
     let n = errors.len();
     let out: Vec<Json> =
         errors.iter().take(1000).map(|e| json!({"rule": e.rule, "layer": e.layer, "rows": e.rows, "lon": e.lon, "lat": e.lat})).collect();

crates/engine/src/cmd/select.rs

@@ -239,17 +239,38 @@ fn by_radius(s: &mut Session, p: &Value) -> CmdResult {
     flat.validate_selection_path(&boundary)?;
     let scale = s.scale();
     let mut hits = Vec::new();
+    let mut bounded = Vec::new();
     for layer in s.map_mut().selectable_layers(scale) {
         let drawn = layer.drawn_rows(scale);
         let rows: Vec<usize> = drawn.iter().enumerate().filter_map(|(row, shown)| shown.then_some(row)).collect();
-        let rows: Vec<usize> = network::within(&layer.dataset, &rows, center, radius).into_iter().map(|(row, _)| row).collect();
+        let (near, n) = network::within_counted(&layer.dataset, &rows, center, radius);
+        if n > 0 {
+            bounded.push((layer.name.clone(), n));
+        }
+        let rows: Vec<usize> = near.into_iter().map(|(row, _)| row).collect();
         if !rows.is_empty() {
             hits.push((layer.id, rows));
         }
     }
+    warn_bounded(s, &bounded);
     apply_hits(s, hits, m)
 }
 
+/// One warning for the features of each layer a radius search ([`network::within_counted`])
+/// counted at a bound on their distance near a pole, with their counts (#1033); none when there are
+/// none.
+pub(crate) fn warn_bounded(s: &mut Session, bounded: &[(String, usize)]) {
+    if bounded.is_empty() {
+        return;
+    }
+    let total: usize = bounded.iter().map(|(_, n)| n).sum();
+    let layers: Vec<String> = bounded.iter().map(|(layer, n)| format!("{layer} {n}")).collect();
+    s.warn(format!(
+        "{total} feature(s) near a pole ({}) couldn't be measured in full (past a search's 250,000 steps a layer, or along an edge over 512° long) and were counted at the least distance they could be: they may lie farther, even beyond the radius",
+        layers.join(", ")
+    ));
+}
+
 /// The hit count describes the query, regardless of whether its rows are added or removed.
 fn apply_hits(s: &mut Session, hits: Vec<(u32, Vec<usize>)>, m: Mode) -> CmdResult {
     if m == Mode::New {

crates/engine/src/cmd/workflows.rs

@@ -723,17 +723,24 @@ fn route(s: &mut Session, p: &Value) -> CmdResult {
 
 fn neighborhood(s: &mut Session, p: &Value) -> CmdResult {
     let at = Coord::new(f(p, "lon").ok_or("lon is required")?, f(p, "lat").ok_or("lat is required")?);
+    if !(at.x.is_finite() && at.y.abs() <= 90.0) {
+        return Err("lon must be a finite number and lat within -90 to 90".into());
+    }
     let radius = f(p, "radius").filter(|r| *r > 0.0).ok_or("radius is required (metres)")?.min(5e6);
     let only: Option<Vec<String>> =
         p.get("layers").and_then(Value::as_array).map(|a| a.iter().filter_map(|v| v.as_str().map(str::to_ascii_lowercase)).collect());
     let mut out = Vec::new();
     let mut found: Vec<(u32, Vec<usize>)> = Vec::new();
+    let mut bounded = Vec::new();
     for l in s.map().layers.iter().filter(|l| s.map().shows(l) && !l.is_raster() && !l.graphics) {
         if only.as_ref().is_some_and(|o| !o.contains(&l.name.to_ascii_lowercase())) {
             continue;
         }
         let rows: Vec<usize> = l.visible_rows().iter().enumerate().filter(|(_, v)| **v).map(|(i, _)| i).collect();
-        let near = network::within(&l.dataset, &rows, at, radius);
+        let (near, n) = network::within_counted(&l.dataset, &rows, at, radius);
+        if n > 0 {
+            bounded.push((l.name.clone(), n));
+        }
         if near.is_empty() {
             continue;
         }
@@ -749,6 +756,7 @@ fn neighborhood(s: &mut Session, p: &Value) -> CmdResult {
         out.push(json!({"layer": l.name, "count": near.len(), "nearestKm": near.first().map(|x| x.1 / 1000.0), "sums": sums}));
         found.push((l.id, near.into_iter().map(|x| x.0).collect()));
     }
+    super::select::warn_bounded(s, &bounded);
     if b(p, "select").unwrap_or(false) {
         s.map_mut().clear_selection();
         for (id, rows) in found {