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 map takes it straight in longitude and latitude, and their percentages and shares follow the map, 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

@@ -1238,8 +1238,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();
@@ -1256,6 +1259,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(),
@@ -1264,6 +1275,12 @@ fn topology_validate(s: &mut Session, p: &Json) -> CmdResult {
             lat: e.at.y,
         }));
     }
+    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

@@ -712,17 +712,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;
         }
@@ -738,6 +745,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 {

crates/engine/src/tests.rs

@@ -6416,6 +6416,25 @@ fn analysis_workflows_models_scripts_and_raster_functions() {
         .run("analysis.neighborhood", &json!({"lon": 2.35, "lat": 48.85, "radius": 300000, "layers": ["Cities"], "select": true}))
         .unwrap();
     assert!(n["layers"][0]["count"].as_u64().unwrap() >= 2);
+    // A center that's no place is refused, not taken as a pole (#1033).
+    assert!(s.run("analysis.neighborhood", &json!({"lon": 2.35, "lat": 1000, "radius": 300000})).is_err());
+    // Lines near a pole with more vertices than one search measures: those past its steps are
+    // counted at the least distance they could be, and it says how many (#1033); a search that
+    // measures every one says nothing.
+    let mut ds = rustgis_data::Dataset::new("Spirals", rustgis_data::GeometryType::Polyline, vec![]);
+    for _ in 0..300 {
+        let line = (0..=348).map(|i| rustgis_data::Coord::new(if i % 2 == 0 { -179.5 } else { 179.5 }, 89.95)).collect();
+        ds.push(Some(rustgis_data::Geometry::Line(vec![line])), vec![]);
+    }
+    s.map_mut().add_dataset(ds);
+    let mark = s.notice_mark();
+    let n = s.run("analysis.neighborhood", &json!({"lon": 0, "lat": 89.9, "radius": 10000, "layers": ["Spirals"]})).unwrap();
+    assert_eq!(n["layers"][0]["count"], 300);
+    let said: Vec<&str> = s.notices_since(mark).iter().map(|n| n.text.as_str()).collect();
+    assert!(said.iter().any(|t| t.starts_with("299 feature(s) near a pole (Spirals 299)") && t.contains("beyond the radius")), "{said:?}");
+    let mark = s.notice_mark();
+    s.run("analysis.neighborhood", &json!({"lon": 0, "lat": 89.9, "radius": 10000, "layers": ["Cities"]})).unwrap();
+    assert!(s.notices_since(mark).iter().all(|n| !n.text.contains("near a pole")), "{:?}", s.notices_since(mark));
     s.run("select.byAttributes", &json!({"layer": "Countries", "where": "name = 'France'"})).unwrap();
     let sim = s.run("analysis.simulate", &json!({"layer": "Countries", "count": 20, "seed": 3})).unwrap();
     assert_eq!(sim["points"], 20);
@@ -24585,3 +24604,32 @@ fn background_runs_follow_a_renamed_map_a_moved_window_and_a_loaded_layer() {
     assert!(s.map().layers[0].dataset.features.iter().any(|f| f.attrs.first().and_then(|v| v.as_f64()) == Some(-1.0)));
     let _ = std::fs::remove_file(path);
 }
+
+#[test]
+fn must_be_covered_by_says_which_polygons_it_left_out_at_a_pole() {
+    // Around the South Pole, a cap from 20°E round to 10°E and the 10° between written a
+    // billion degrees on: the check at the pole can't use that one, reports the point there and
+    // warns why; with a whole cap it says nothing of the kind (#1033).
+    use rustgis_data::{Coord, Dataset, Geometry, GeometryType};
+    let c = Coord::new;
+    let ring = |w: f64, e: f64| Geometry::Polygon(vec![vec![vec![c(w, -85.0), c(e, -85.0), c(e, -90.0), c(w, -90.0), c(w, -85.0)]]]);
+    for (caps, warned) in [(vec![ring(20.0, 370.0), ring(1e9 + 90.0, 1e9 + 100.0)], true), (vec![ring(-180.0, 180.0)], false)] {
+        let mut s = Session::new(Project::default());
+        let mut by = Dataset::new("Caps", GeometryType::Polygon, vec![]);
+        for g in caps {
+            by.push(Some(g), vec![]);
+        }
+        let mut pts = Dataset::new("Pole", GeometryType::Point, vec![]);
+        pts.push(Some(Geometry::Point(c(0.0, -90.0))), vec![]);
+        s.map_mut().add_dataset(by);
+        s.map_mut().add_dataset(pts);
+        s.run("topology.rules", &json!({"add": {"rule": "mustBeCoveredBy", "layer": "Pole", "by": "Caps"}})).unwrap();
+        let mark = s.notice_mark();
+        let v = s.run("topology.validate", &Value::Null).unwrap();
+        assert_eq!(v["errors"], u64::from(warned));
+        let said = s.notices_since(mark).iter().any(|n| {
+            n.text.contains("left out of the check at a pole") && n.text.contains("Must Be Covered By on 'Pole': 1 feature(s) of 'Caps'")
+        });
+        assert_eq!(said, warned, "{:?}", s.notices_since(mark));
+    }
+}