fix(analysis): keep polygon simplification monotonic at large tolerances

crates/analysis/src/tools.rs

@@ -892,11 +892,54 @@ fn xy_table_to_point(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
     Ok(ToolResult { output: Output::Layer(out), messages })
 }
 
+// Three fixed non-collinear anchors prevent ring collapse and keep vertex counts monotonic.
+fn simplify_ring(ring: &geo::LineString<f64>, eps: f64) -> geo::LineString<f64> {
+    use geo::{Distance, Euclidean};
+    if ring.0.len() <= 4 || ring.0.iter().any(|c| !c.x.is_finite() || !c.y.is_finite()) {
+        return ring.clone();
+    }
+    let Some(&first) = ring.0.first() else { return ring.clone() };
+    let Some((b, _)) = ring.0.iter().enumerate().map(|(i, c)| (i, Euclidean.distance(first, *c))).max_by(|a, b| a.1.total_cmp(&b.1)) else {
+        return ring.clone();
+    };
+    let Some(&second) = ring.0.get(b) else { return ring.clone() };
+    let Some((c, _)) = ring
+        .0
+        .iter()
+        .enumerate()
+        .map(|(i, c)| (i, ((second.x - first.x) * (c.y - first.y) - (second.y - first.y) * (c.x - first.x)).abs()))
+        .filter(|(_, area)| area.is_finite() && *area > 0.0)
+        .max_by(|a, b| a.1.total_cmp(&b.1))
+    else {
+        return ring.clone();
+    };
+    let mut stops = [0, b, c, ring.0.len() - 1];
+    stops.sort_unstable();
+    let mut corners = Vec::new();
+    for bounds in stops.windows(2) {
+        let [start, end] = bounds else { return ring.clone() };
+        let Some(arc) = ring.0.get(*start..=*end) else { return ring.clone() };
+        let simple = geo::LineString::new(arc.to_vec()).simplify(eps);
+        corners.extend(simple.0.into_iter().skip(usize::from(!corners.is_empty())));
+    }
+    geo::LineString::new(corners)
+}
+
 /// A geometry simplified (Douglas–Peucker) to a tolerance in metres.
 pub fn simplify_geometry(g: &Geometry, tol_m: f64) -> Geometry {
     let eps = tol_m * merc_scale(feature_lat(g));
     let sg = match g {
-        Geometry::Polygon(_) => to_multipolygon(g).and_then(|m| from_multipolygon(&m.simplify(eps))),
+        Geometry::Polygon(_) => to_multipolygon(g).and_then(|mut m| {
+            for p in &mut m.0 {
+                p.exterior_mut(|r| *r = simplify_ring(r, eps));
+                p.interiors_mut(|rings| {
+                    for r in rings {
+                        *r = simplify_ring(r, eps);
+                    }
+                });
+            }
+            from_multipolygon(&m)
+        }),
         Geometry::Line(_) => to_multiline(g).and_then(|m| from_line(&m.simplify(eps))),
         _ => None,
     };
@@ -912,14 +955,9 @@ fn simplify(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
         let Some(f) = ds.features.get(r) else { continue };
         let Some(g) = &f.geometry else { continue };
         before += g.coord_count();
-        let eps = tol * merc_scale(feature_lat(g));
-        let sg = match g {
-            Geometry::Polygon(_) => to_multipolygon(g).and_then(|m| from_multipolygon(&m.simplify(eps))),
-            Geometry::Line(_) => to_multiline(g).and_then(|m| from_line(&m.simplify(eps))),
-            other => Some(other.clone()),
-        };
-        after += sg.as_ref().map(|g| g.coord_count()).unwrap_or(0);
-        out.push(sg, f.attrs.clone());
+        let sg = simplify_geometry(g, tol);
+        after += sg.coord_count();
+        out.push(Some(sg), f.attrs.clone());
     }
     result(out, format!("Simplify: {before} → {after} vertices"))
 }
@@ -2091,6 +2129,82 @@ mod tests {
         }
     }
 
+    #[test]
+    fn simplify_rings_do_not_grow_with_tolerance() {
+        use geo::{Distance, Euclidean};
+        let ring = |x: f64, y: f64, radius: f64| {
+            (0..=80)
+                .map(|i| {
+                    let angle = (i % 80) as f64 * std::f64::consts::TAU / 80.0;
+                    let r = radius * (1.0 + 0.1 * (7.0 * angle).sin());
+                    Coord::new(x + r * angle.cos(), y + r * angle.sin())
+                })
+                .collect::<Vec<_>>()
+        };
+        let original = Geometry::Polygon(vec![vec![ring(10.5, 46.5, 0.1), ring(10.5, 46.52, 0.001)], vec![ring(11.0, 46.5, 0.05)]]);
+        let mut ds = Dataset::new("zones", GeometryType::Polygon, vec![]);
+        ds.push(Some(original.clone()), vec![]);
+        let inputs = Map(HashMap::from([("zones".into(), ds)]));
+        let mut previous = original.coord_count();
+        for tolerance in [0.0, 100.0, 500.0, 1000.0, 2000.0, 5000.0, 25000.0, 100000.0] {
+            let simple = simplify_geometry(&original, tolerance);
+            assert!(simple.coord_count() <= previous, "{tolerance} m: {} > {previous}", simple.coord_count());
+            previous = simple.coord_count();
+            let Geometry::Polygon(parts) = &simple else { panic!("polygon required") };
+            assert_eq!(parts.iter().map(Vec::len).collect::<Vec<_>>(), [2, 1]);
+            for rings in parts {
+                for r in rings {
+                    assert!(r.len() >= 4);
+                    assert_eq!(r.first(), r.last());
+                    assert!(rustgis_data::geometry::ring_area(r).abs() > 1e-14);
+                }
+            }
+            let source = to_multipolygon(&original).unwrap();
+            let target = to_multipolygon(&simple).unwrap();
+            for (a, b) in source.0.iter().zip(&target.0) {
+                for (a, b) in std::iter::once(a.exterior()).chain(a.interiors()).zip(std::iter::once(b.exterior()).chain(b.interiors())) {
+                    let deviation = a.0.iter().map(|c| Euclidean.distance(&geo::Point::from(*c), b)).fold(0.0_f64, f64::max);
+                    assert!(deviation <= tolerance * merc_scale(feature_lat(&original)) + 1e-6, "{tolerance} m: {deviation}");
+                }
+            }
+            let out =
+                layer(crate::run("simplify", &inputs, &serde_json::json!({"input":"zones", "tolerance":tolerance / 1000.0})).unwrap());
+            assert_eq!(out.features[0].geometry.as_ref(), Some(&simple));
+        }
+        assert_eq!(previous, 12, "three closed triangle rings at a very large tolerance");
+    }
+
+    #[test]
+    fn simplify_keeps_points_lines_and_degenerate_rings() {
+        use geo::Area;
+        let point = Geometry::Point(Coord::new(10.0, 46.0));
+        assert_eq!(simplify_geometry(&point, 1e6), point);
+        let line = Geometry::Line(vec![vec![Coord::new(10.0, 46.0), Coord::new(10.1, 46.0), Coord::new(10.2, 46.0)]]);
+        assert_eq!(simplify_geometry(&line, 1e6).coord_count(), 2);
+        let repeated = geo::LineString::new(vec![geo::Coord { x: 1.0, y: 1.0 }; 6]);
+        assert_eq!(simplify_ring(&repeated, 1e6), repeated);
+        let collinear = geo::LineString::from(vec![(0.0, 0.0), (1.0, 0.0), (2.0, 0.0), (3.0, 0.0), (0.0, 0.0)]);
+        assert_eq!(simplify_ring(&collinear, 1e6), collinear);
+        let both_sides = geo::LineString::from(vec![(0.0, 0.0), (3.0, 0.0), (-3.0, 0.0), (0.0, 0.0), (1.0, 0.0), (0.0, 0.0)]);
+        assert_eq!(simplify_ring(&both_sides, 1e6), both_sides);
+        let concave = geo::LineString::from(vec![(0.0, 0.0), (3.0, 0.0), (2.0, 1.0), (-2.0, 1.0), (-3.0, 0.0), (-2.0, -1.0), (0.0, 0.0)]);
+        let mut previous = concave.0.len();
+        for tolerance in [0.1, 1.1, 2.5, 1e6] {
+            use geo::{Distance, Euclidean};
+            let reduced = simplify_ring(&concave, tolerance);
+            assert!((4..=previous).contains(&reduced.0.len()));
+            previous = reduced.0.len();
+            let deviation = concave.0.iter().map(|c| Euclidean.distance(&geo::Point::from(*c), &reduced)).fold(0.0_f64, f64::max);
+            assert!(deviation <= tolerance + 1e-12);
+            assert!(geo::Polygon::new(reduced, vec![]).unsigned_area() > 0.0);
+        }
+        assert_eq!(previous, 4);
+        let triangle = geo::LineString::from(vec![(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (0.0, 0.0)]);
+        assert_eq!(simplify_ring(&triangle, 1e6), triangle);
+        let nonfinite = geo::LineString::from(vec![(0.0, 0.0), (f64::INFINITY, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]);
+        assert_eq!(simplify_ring(&nonfinite, 1e6), nonfinite);
+    }
+
     #[test]
     fn all_tools_run() {
         let m = inputs();

crates/engine/src/cmd/layout.rs

@@ -1439,10 +1439,26 @@ fn publication(s: &mut Session, p: &Value) -> CmdResult {
             // Each polygon simplified to a tenth of a millimetre on the disc (its scale in metres),
             // so the project keeps a sketch of the layers for it, not a second full copy of them.
             let tolerance_m = globe * 1e-4;
-            // A polygon too small for that (an island) would collapse to a ring of under four
-            // points: it keeps its own shape, already small.
+            // On the locator, sub-tolerance rings keep their outline rather than becoming triangles.
             let sketch = |g: &rustgis_data::Geometry| {
-                let simple = rustgis_analysis::simplify_geometry(g, tolerance_m);
+                let mut simple = rustgis_analysis::simplify_geometry(g, tolerance_m);
+                if let (rustgis_data::Geometry::Polygon(original), rustgis_data::Geometry::Polygon(reduced)) = (g, &mut simple) {
+                    for (original, reduced) in original.iter().flatten().zip(reduced.iter_mut().flatten()) {
+                        let mut extent = Extent::EMPTY;
+                        // Keep the extent continuous across the antimeridian.
+                        let mut longitude = original.iter().find(|c| c.is_finite()).map_or(0.0, |c| c.x);
+                        for &c in original.iter().filter(|c| c.is_finite()) {
+                            longitude += (c.x - longitude + 180.0).rem_euclid(360.0) - 180.0;
+                            extent.expand(Coord::new(longitude, c.y));
+                        }
+                        let projected = rustgis_geom::extent_to_merc(&extent);
+                        let size_m =
+                            projected.width().hypot(projected.height()) * extent.center().y.clamp(-85.0, 85.0).to_radians().cos().max(0.01);
+                        if size_m < tolerance_m {
+                            *reduced = original.clone();
+                        }
+                    }
+                }
                 match &simple {
                     rustgis_data::Geometry::Polygon(polys) if polys.iter().flatten().all(|ring| ring.len() >= 4) => simple,
                     _ => g.clone(),

crates/engine/src/tests.rs

@@ -16858,3 +16858,59 @@ fn raster_projection_failures_are_reported_for_zonal_statistics_and_clipping() {
     let data = s.map().layer_by_name("clipped").unwrap().raster.as_ref().unwrap().data.as_ref().unwrap();
     assert_eq!((data.width, data.height), (8, 8));
 }
+
+#[test]
+fn simplify_large_tolerances_keep_a_small_polygon_and_one_step_undo() {
+    use rustgis_data::{Dataset, Field, FieldType, Geometry, GeometryType, Value};
+    use rustgis_geom::Coord;
+    let mut s = world();
+    let ring = (0..=40)
+        .map(|i| {
+            let angle = (i % 40) as f64 * std::f64::consts::TAU / 40.0;
+            Coord::new(10.5 + 0.05 * angle.cos(), 46.5 + 0.05 * angle.sin())
+        })
+        .collect();
+    let mut zones = Dataset::new("zones", GeometryType::Polygon, vec![Field::new("code", FieldType::Text)]);
+    zones.push(Some(Geometry::Polygon(vec![vec![ring]])), vec![Value::Text("Z1".into())]);
+    s.map_mut().add_layer(rustgis_map::Layer::new(0, zones.clone()));
+    let mut previous = zones.features[0].geometry.as_ref().unwrap().coord_count();
+    for tolerance in [0.1, 0.5, 1.0, 2.0, 5.0, 25.0, 100.0] {
+        let depth = s.undo_depth();
+        s.run("gp.run", &json!({"tool":"simplify", "params":{"input":"zones", "tolerance":tolerance, "output":"simple"}})).unwrap();
+        let out = &s.map().layer_by_name("simple").unwrap().dataset;
+        assert_eq!(out.features.len(), 1);
+        assert_eq!(out.features[0].attrs, zones.features[0].attrs);
+        let count = out.features[0].geometry.as_ref().unwrap().coord_count();
+        assert!((4..=previous).contains(&count), "{tolerance} km: {count} > {previous}");
+        previous = count;
+        assert_eq!(s.map().layer_by_name("zones").unwrap().dataset, zones);
+        assert_eq!(s.undo_depth(), depth + 1);
+        assert!(s.undo());
+        assert!(s.map().layer_by_name("simple").is_none());
+        assert!(s.redo());
+        assert_eq!(s.map().layer_by_name("simple").unwrap().dataset.features[0].geometry.as_ref().unwrap().coord_count(), count);
+        assert!(s.undo());
+    }
+    assert_eq!(previous, 4);
+}
+
+#[test]
+fn publication_locator_keeps_small_islands_across_the_antimeridian() {
+    use rustgis_data::{Dataset, Geometry, GeometryType};
+    use rustgis_geom::Coord;
+    for (lon, lat) in [(0.0_f64, -5.0), (180.0, -5.0), (-180.0, 70.0)] {
+        let mut s = world();
+        let ring: Vec<_> = [(-0.0005, 0.0), (0.0005, 0.0), (0.0005, 0.001), (-0.0005, 0.001), (-0.0005, 0.0)]
+            .into_iter()
+            .map(|(x, y)| Coord::new((lon + x + 180.0).rem_euclid(360.0) - 180.0, lat + y))
+            .collect();
+        let original = Geometry::Polygon(vec![vec![ring]]);
+        let mut island = Dataset::new("island", GeometryType::Polygon, vec![]);
+        island.push(Some(original.clone()), vec![]);
+        s.map_mut().add_layer(rustgis_map::Layer::new(0, island));
+        let r = s.run("layout.publication", &json!({"extent": [-46.8, -11.6, -38.5, -2.2]})).unwrap();
+        let locator = s.project.maps.iter().find(|m| Some(m.name.as_str()) == r["locator"].as_str()).unwrap();
+        let geometry = &locator.layer_by_name("island").unwrap().dataset.features[0].geometry;
+        assert_eq!(geometry.as_ref(), Some(&original), "island at {lon}, {lat}");
+    }
+}

docs/user-guide/analysis.md

@@ -104,6 +104,13 @@ Use these ids in `gp.run`'s `tool` parameter (`gp.tools` lists the available too
 | IDW Interpolation | `idwInterpolation` | Analysis › Interpolation |
 | Ordinary Kriging | `kriging` | Analysis › Interpolation |
 
+**Simplify** uses Douglas–Peucker with a tolerance in kilometres. Each polygon ring keeps at
+least three different, non-collinear corners (plus the closing vertex), so a tolerance larger
+than a polygon produces a small ring instead of restoring its original vertices. Holes and
+multipart rings are retained. This is not a topology-preserving simplifier: check the result
+where holes or neighbouring polygons must stay aligned. The default tolerance remains 25 km;
+use a smaller value for local boundaries.
+
 Most feature tools process only the selected rows when an input has a selection. Selection tools
 instead replace the input's selection; Select Layer By Location honours the selecting layer's
 selection. Summarize Within counts the selected points against **all** input polygons, and