diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index b9ec91289..ca3eb561a 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -892,11 +892,54 @@ fn xy_table_to_point(inp: &dyn Inputs, p: &Json) -> Result { 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, eps: f64) -> geo::LineString { + 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 { 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::>() + }; + 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::>(), [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(); diff --git a/crates/engine/src/cmd/layout.rs b/crates/engine/src/cmd/layout.rs index 69b3f145a..8680650c2 100644 --- a/crates/engine/src/cmd/layout.rs +++ b/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(), diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index e3a15e5da..655c681ba 100644 --- a/crates/engine/src/tests.rs +++ b/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}"); + } +} diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 8ae08fb11..1c023f19c 100644 --- a/docs/user-guide/analysis.md +++ b/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