From 85bb76cb9acb34a1e94537b1d827cba75af2c434 Mon Sep 17 00:00:00 2001 From: mavrick Date: Mon, 26 Oct 2026 13:28:41 -0500 Subject: [PATCH 1/5] feat(analysis): add an Intersect overlay tool Intersect splits input features (points, lines or polygons) by the polygons of an overlay layer: one feature per part inside an overlay polygon, with the input's attributes followed by the overlay's (a taken name gets a suffix). An overlay edge counts as inside, as for Clip. Selections in both layers are honoured. Clip's per-geometry logic becomes part_in, shared with Intersect, and Spatial Join's extent grid becomes ExtentIndex, so Intersect prunes candidates the same way. The GEOS reference corpus gains intersect cases, a road along a zone's edge, and keys of several fields. Part of #959. Co-Authored-By: Claude Opus 5.5 --- GEOLIBRE_PARITY.md | 4 +- README.md | 2 +- ROADMAP.md | 16 +- crates/analysis/src/tools.rs | 390 ++++++++++++++---- crates/analysis/tests/reference.rs | 48 ++- crates/analysis/tests/reference/cases.json | 89 +++- .../reference/expected/clip_roads.geojson | 23 ++ .../expected/intersect_parcels.geojson | 277 +++++++++++++ .../expected/intersect_roads.geojson | 138 +++++++ crates/analysis/tests/reference/generate.py | 22 + .../tests/reference/layers/roads.geojson | 19 + .../assistant/skills/geoprocessing/SKILL.md | 3 +- crates/engine/src/cmd/workflows.rs | 15 +- crates/ui-egui/src/icons.rs | 7 + crates/ui-egui/src/panes.rs | 1 + crates/ui-egui/src/ribbon.rs | 1 + docs/comparison/index.md | 6 +- docs/user-guide/analysis.md | 20 +- 18 files changed, 954 insertions(+), 127 deletions(-) create mode 100644 crates/analysis/tests/reference/expected/intersect_parcels.geojson create mode 100644 crates/analysis/tests/reference/expected/intersect_roads.geojson diff --git a/GEOLIBRE_PARITY.md b/GEOLIBRE_PARITY.md index 125f07c96..f3010b864 100644 --- a/GEOLIBRE_PARITY.md +++ b/GEOLIBRE_PARITY.md @@ -159,8 +159,8 @@ _GeoLibre: docs/features.md › Attribute data and expressions; docs/arcgis-edit _GeoLibre: docs/features.md › Processing and analysis; docs/user-guide/processing.md, model-builder.md, segmentation.md, obia.md_ -- [ ] **Spatial join, intersect, union**. Partial: spatial join, clip and erase are built in; - intersect, union, identity and near are missing (_Roadmap: solid ground_). +- [ ] **Spatial join, intersect, union**. Partial: spatial join, clip, erase and intersect are + built in; union, identity and near are missing (_Roadmap: solid ground_). - [ ] **More built-in vector tools**: attribute join, merge, extract vertices, points along lines, fix/check geometry, random extract, buffer options. Partial (_Roadmap: solid ground_). - [ ] **Routing services**: isochrones, service areas, OD matrix, multi-stop routing, diff --git a/README.md b/README.md index 9c631fcda..a2e432fca 100644 --- a/README.md +++ b/README.md @@ -99,7 +99,7 @@ Every screenshot is RustGIS itself, rendered offscreen by `cargo run -p rustgis- | **Web tiles** | XYZ/TMS/WMS raster tiles, Mapbox Vector Tiles and MapLibre styles (sprites, patterns, expressions), cached and fetched in the background | | **Cartography** | Single symbol, unique values, graduated colours, proportional symbols; label classes with SQL queries, straight and curved line placement, installed font families, scale ranges and annotation; stacked and dashed line strokes, symbol levels, reference scales and polygons without fill; masking, display filters, aggregation bins, 2.5D extrusion; graphics and dynamic text | | **Editing** | Create, move, rotate, scale, merge, split and explode features; snapping, templates, elevation, traverse, topology rules; domains, subtypes, attribute rules, relationships, archiving; undo for everything | -| **Analysis** | Selection by attributes (SQL) and location; Buffer, Clip, Erase, Dissolve, Summarize Within, Summary Statistics, Spatial Join, Zonal Statistics, A5 grids and a geometry calculator; viewshed, line of sight, IDW, ordinary kriging, shortest path, suitability modelling; ModelBuilder on a node canvas and scripts with command, parameter and value completion; 1,000+ [GeoLibre](https://github.com/opengeos/geolibre-rust) and Whitebox raster, vector, LiDAR, hydrology and terrain tools | +| **Analysis** | Selection by attributes (SQL) and location; Buffer, Clip, Erase, Intersect, Dissolve, Summarize Within, Summary Statistics, Spatial Join, Zonal Statistics, A5 grids and a geometry calculator; viewshed, line of sight, IDW, ordinary kriging, shortest path, suitability modelling; ModelBuilder on a node canvas and scripts with command, parameter and value completion; 1,000+ [GeoLibre](https://github.com/opengeos/geolibre-rust) and Whitebox raster, vector, LiDAR, hydrology and terrain tools | | **Imagery** | Ortho and reality mapping, pixel classification (k-means, minimum distance, a trainable neural network), mensuration, hosted imagery tiles | | **Layouts and reports** | Page layouts with map frames, legends, scale bars, grids and map series; publication sheets and figures in one command, format panes, guides, snapping and element locks; reports with groups, statistics, charts and subreports; presentations (HTML slideshow, PDF); notebooks (Jupyter in and out) | | **Sharing** | PNG export at any DPI; PDF; web map pages, OGC API Features collections, MapLibre styles, glTF scenes, locators; project, map and mobile map packages; offline areas with sync | diff --git a/ROADMAP.md b/ROADMAP.md index 75f973741..2ac014de1 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -60,7 +60,8 @@ Correct, robust and recoverable, with the foundations other themes stand on. administrator can fix for a managed install (no AI Assistant, plugins or web services). - Extend the OS keychain secret store to service connections (database passwords, S3 keys, AI provider keys and the ChatGPT sign-in's tokens are already kept there). -- Overlay tools beside Clip and Erase: intersect, union, identity and near _(open)_. +- Overlay tools beside Clip, Erase and Intersect: union, identity and near + ([#959](https://github.com/opengeos/RustGIS/issues/959)). - Spatial Join by distance, closest, centre in and largest overlap, with merge rules (sum, mean, first…) for one-to-one joins ([#851](https://github.com/opengeos/RustGIS/issues/851)). - Vector tools, each writing a new layer: join field (a one-off attribute join), merge (append @@ -121,8 +122,7 @@ And alongside it: a line layer). - Hot-spot statistics (Getis-Ord Gi\*), and emerging hot spots on the time slider. - Charts linked to the selection: histogram, scatter, bar and box plots of a layer's fields. -- Tabulate intersection: one layer's counts and totals shared out to another by area overlap, - once intersect exists. +- Tabulate intersection: one layer's counts and totals shared out to another by area overlap. - Reverse geocoding, and the US Census geocoder beside Nominatim and Photon; a DEM downloader so viewshed works without a loaded elevation raster _(open)_. @@ -220,11 +220,11 @@ Open questions rather than plans; discuss before building. templates, elevation, topology; domains, subtypes, attribute rules, contingent values, relationships, archiving. - Analysis: selection by attributes and location; polygon, lasso and radius on the map; overlay - (clip, erase), proximity (buffer) and summary tools, spatial join, zonal statistics, A5 grids, - kriging and Calculate Geometry with metric or imperial units; guided workflows (viewshed, line - of sight, IDW, shortest path, suitability); ModelBuilder on a node canvas, scripts with command, - parameter and value completion, and history copied as code; the GeoLibre and Whitebox suite - (1,000+ tools) natively and in the browser. + (clip, erase, intersect), proximity (buffer) and summary tools, spatial join, zonal statistics, + A5 grids, kriging and Calculate Geometry with metric or imperial units; guided workflows + (viewshed, line of sight, IDW, shortest path, suitability); ModelBuilder on a node canvas, + scripts with command, parameter and value completion, and history copied as code; the GeoLibre + and Whitebox suite (1,000+ tools) natively and in the browser. - Layouts (map frames, surrounds, map series, publication sheets and figures in one command, format panes, guides, snapping, element locking and PDF), reports, presentations and notebooks; the web app embeddable in a page (`?embed`). diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index 956361a52..f934da209 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -58,6 +58,18 @@ pub static TOOLS: &[ToolSpec] = &[ ], run: erase, }, + ToolSpec { + id: "intersect", + name: "Intersect", + toolbox: "Analysis › Overlay", + summary: "Splits input features where they overlap overlay polygons, keeping both layers' attributes.", + params: &[ + P { name: "input", label: "Input Features", kind: K::Layer, required: true, default: "", help: "" }, + P { name: "overlay", label: "Overlay Features", kind: K::PolygonLayer, required: true, default: "", help: "" }, + P { name: "output", label: "Output Name", kind: K::Text, required: false, default: "", help: "" }, + ], + run: intersect, + }, ToolSpec { id: "dissolve", name: "Dissolve", @@ -679,37 +691,43 @@ fn union_of(inp: &dyn Inputs, p: &Json, k: &str) -> Result, St Ok(union_all(polys.iter().flat_map(|m| m.0.iter()))) } +/// The part of a feature inside `mask` (polygons in Web Mercator, oriented as `union_all` leaves +/// them): a polygon's overlap with it, a line's stretches in it, the points in it or on its edge. +/// `None` when nothing is left, or only a shared edge or corner. +fn part_in(g: &Geometry, mask: &MultiPolygon) -> Option { + match g { + Geometry::Polygon(_) => to_multipolygon(g).and_then(|m| from_multipolygon(&m.intersection(mask))), + Geometry::Line(_) => to_multiline(g).and_then(|m| from_line(&mask.clip(&m, false))), + Geometry::Point(_) | Geometry::MultiPoint(_) => { + let pts: &[Coord] = match g { + Geometry::Point(c) => std::slice::from_ref(c), + Geometry::MultiPoint(v) => v, + _ => &[], + }; + let inside: Vec = pts + .iter() + .copied() + .filter(|c| { + let m = rustgis_geom::lonlat_to_merc(*c); + mask.intersects(&geo::Point::new(m.x, m.y)) + }) + .collect(); + match inside.len() { + 0 => None, + 1 => inside.first().map(|c| Geometry::Point(*c)), + _ => Some(Geometry::MultiPoint(inside)), + } + } + } +} + fn clip(inp: &dyn Inputs, p: &Json) -> Result { let ds = input(inp, p, "input")?; let mask = union_of(inp, p, "clip")?; let mut out = copy_schema(ds, out_name(p, format!("{}_Clip", ds.name)), ds.geometry_type); for r in rows(inp, p, "input", ds) { let Some(f) = ds.features.get(r) else { continue }; - let Some(g) = &f.geometry else { continue }; - let clipped = match g { - Geometry::Polygon(_) => to_multipolygon(g).and_then(|m| from_multipolygon(&m.intersection(&mask))), - Geometry::Line(_) => to_multiline(g).and_then(|m| from_line(&mask.clip(&m, false))), - Geometry::Point(_) | Geometry::MultiPoint(_) => { - let pts: Vec = match g { - Geometry::Point(c) => vec![*c], - Geometry::MultiPoint(v) => v.clone(), - _ => vec![], - }; - let inside: Vec = pts - .into_iter() - .filter(|c| { - let m = rustgis_geom::lonlat_to_merc(*c); - mask.intersects(&geo::Point::new(m.x, m.y)) - }) - .collect(); - match inside.len() { - 0 => None, - 1 => inside.first().map(|c| Geometry::Point(*c)), - _ => Some(Geometry::MultiPoint(inside)), - } - } - }; - if let Some(c) = clipped { + if let Some(c) = f.geometry.as_ref().and_then(|g| part_in(g, &mask)) { out.push(Some(c), f.attrs.clone()); } } @@ -731,6 +749,81 @@ fn erase(inp: &dyn Inputs, p: &Json) -> Result { result(out, format!("Erase: {n} features remain")) } +/// An overlay polygon: its row in the overlay layer, its extent and its polygons in Web Mercator, +/// oriented for the overlay. +struct Zone { + row: usize, + ext: Extent, + poly: MultiPolygon, +} + +/// The overlay layer of an overlay tool (its `overlay` parameter), which must be polygons, and its +/// (selected) features with geometry as [`Zone`]s, in order. +fn zones<'a>(inp: &'a dyn Inputs, p: &Json) -> Result<(&'a Dataset, Vec), String> { + let ov = input(inp, p, "overlay")?; + if ov.geometry_type != GeometryType::Polygon { + return Err(format!("'{}' is not a polygon layer", ov.name)); + } + let zones = rows(inp, p, "overlay", ov) + .into_iter() + .filter_map(|row| { + let g = ov.features.get(row)?.geometry.as_ref()?; + Some(Zone { row, ext: g.extent(), poly: to_multipolygon(g)?.orient(Direction::Default) }) + }) + .collect(); + Ok((ov, zones)) +} + +/// An overlay tool's output layer: the input's fields, then the overlay's, an overlay field whose +/// name is taken getting a suffix (`name_2`). It keeps the input's source CRS for export. +fn overlay_schema(ds: &Dataset, ov: &Dataset, name: String, g: GeometryType) -> Dataset { + let mut out = copy_schema(ds, name, g); + out.crs = ds.crs; + for f in &ov.fields { + let name = unique_field(&out, &f.name); + out.fields.push(Field { name, ..f.clone() }); + } + out +} + +/// A row of an overlay tool's output: the input feature's values, then the overlay feature's, each +/// padded to its layer's fields (all empty for a side with no feature). +fn overlay_row(ds: &Dataset, a: Option<&[Value]>, ov: &Dataset, b: Option<&[Value]>) -> Vec { + let mut row = Vec::with_capacity(ds.fields.len() + ov.fields.len()); + for (values, n) in [(a, ds.fields.len()), (b, ov.fields.len())] { + row.extend((0..n).map(|i| values.and_then(|v| v.get(i)).cloned().unwrap_or(Value::Null))); + } + row +} + +/// Intersect: a feature for each part of an input feature inside an overlay polygon, with the +/// attributes of both, in input order and then overlay order. Parts that only touch are left out. +fn intersect(inp: &dyn Inputs, p: &Json) -> Result { + let ds = input(inp, p, "input")?; + if ds.geometry_type == GeometryType::None { + return Err(format!("'{}' has no geometry", ds.name)); + } + let (ov, zones) = zones(inp, p)?; + let index = ExtentIndex::new(zones.iter().map(|z| z.ext)); + let mut out = overlay_schema(ds, ov, out_name(p, format!("{}_Intersect", ds.name)), ds.geometry_type); + for r in rows(inp, p, "input", ds) { + let Some(f) = ds.features.get(r) else { continue }; + let Some(g) = &f.geometry else { continue }; + let ext = g.extent(); + for z in index.candidates(&ext).into_iter().filter_map(|i| zones.get(i)) { + if !ext.intersects(&z.ext) { + continue; + } + if let Some(part) = part_in(g, &z.poly) { + let attrs = overlay_row(ds, Some(&f.attrs), ov, ov.features.get(z.row).map(|o| o.attrs.as_slice())); + out.push_next(Some(part), attrs); + } + } + } + let n = out.features.len(); + result(out, format!("Intersect: {n} features")) +} + fn dissolve(inp: &dyn Inputs, p: &Json) -> Result { let ds = input(inp, p, "input")?; if ds.geometry_type != GeometryType::Polygon { @@ -1472,6 +1565,94 @@ fn unique_field(out: &Dataset, base: &str) -> String { name } +/// Features' extents on a grid, to find those whose extents may meet another extent without +/// testing every one. Up to 16 features, or when the grid would hold more than a million +/// references (large features on a fine grid), there is no grid and every feature is a candidate. +struct ExtentIndex { + /// The union of the extents. + total: Extent, + cell_w: f64, + cell_h: f64, + /// The features each cell meets, in order; `None` without a grid. + grid: Option>>, + /// How many features there are. + len: usize, +} + +impl ExtentIndex { + /// The most cell references the grid may hold. + const MAX_REFERENCES: usize = 1_000_000; + + /// An index of these extents, numbered in order from 0. + fn new(extents: impl Iterator) -> Self { + let extents: Vec = extents.collect(); + let total = extents.iter().fold(Extent::EMPTY, |t, e| t.union(e)); + let n_cells = (extents.len() as f64).sqrt().clamp(4.0, 256.0); + let mut index = ExtentIndex { + total, + cell_w: (total.width() / n_cells).max(1e-9), + cell_h: (total.height() / n_cells).max(1e-9), + grid: None, + len: extents.len(), + }; + if extents.len() <= 16 || total.is_empty() { + return index; + } + let mut references = 0usize; + for e in &extents { + let ((x0, y0), (x1, y1)) = (index.key(e.min_x, e.min_y), index.key(e.max_x, e.max_y)); + references = references.saturating_add(((x1 - x0 + 1) as usize).saturating_mul((y1 - y0 + 1) as usize)); + if references > Self::MAX_REFERENCES { + return index; + } + } + let mut grid: std::collections::HashMap<(i64, i64), Vec> = std::collections::HashMap::new(); + for (i, e) in extents.iter().enumerate() { + let ((x0, y0), (x1, y1)) = (index.key(e.min_x, e.min_y), index.key(e.max_x, e.max_y)); + for x in x0..=x1 { + for y in y0..=y1 { + grid.entry((x, y)).or_default().push(i); + } + } + } + index.grid = Some(grid); + index + } + + /// The cell holding a position, clamped to the grid (a non-finite one goes in the first). + fn key(&self, x: f64, y: f64) -> (i64, i64) { + let kx = if x.is_finite() { ((x - self.total.min_x) / self.cell_w).floor() } else { 0.0 }; + let ky = if y.is_finite() { ((y - self.total.min_y) / self.cell_h).floor() } else { 0.0 }; + (kx.clamp(0.0, 512.0) as i64, ky.clamp(0.0, 512.0) as i64) + } + + /// The larger side of a cell in degrees, when there is a grid. + fn cell_size(&self) -> Option { + self.grid.as_ref().map(|_| self.cell_w.max(self.cell_h)) + } + + /// The features whose extents may meet `ext`, in order and once each: every feature without a + /// grid. Each still needs its own test. + fn candidates(&self, ext: &Extent) -> Vec { + let Some(grid) = &self.grid else { return (0..self.len).collect() }; + if !ext.intersects(&self.total) { + return Vec::new(); + } + let ((x0, y0), (x1, y1)) = (self.key(ext.min_x, ext.min_y), self.key(ext.max_x, ext.max_y)); + let mut found = Vec::new(); + for x in x0..=x1 { + for y in y0..=y1 { + if let Some(indices) = grid.get(&(x, y)) { + found.extend_from_slice(indices); + } + } + } + found.sort_unstable(); + found.dedup(); + found + } +} + /// Joins normalized lon/lat geometries, retaining the target's source CRS for export. fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result { let target = input(inp, p, "target")?; @@ -1561,71 +1742,9 @@ fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result { let target_fields_len = target.fields.len(); let join_fields_len = join.fields.len(); - // Spatial grid index over join geometries for efficient candidate pruning - let mut total_ext = Extent::EMPTY; - for j in &join_geoms { - total_ext = total_ext.union(&j.ext); - } - let mut use_grid = join_geoms.len() > 16 && !total_ext.is_empty(); - let n_cells = (join_geoms.len() as f64).sqrt().clamp(4.0, 256.0); - let cell_w = (total_ext.width() / n_cells).max(1e-9); - let cell_h = (total_ext.height() / n_cells).max(1e-9); - let key = |x: f64, y: f64| -> (i64, i64) { - let kx = if x.is_finite() { ((x - total_ext.min_x) / cell_w).floor() } else { 0.0 }; - let ky = if y.is_finite() { ((y - total_ext.min_y) / cell_h).floor() } else { 0.0 }; - (kx.clamp(0.0, 512.0) as i64, ky.clamp(0.0, 512.0) as i64) - }; - - const MAX_GRID_REFERENCES: usize = 1_000_000; - if use_grid { - let mut projected_refs = 0usize; - for j in &join_geoms { - let (gx0, gy0) = key(j.ext.min_x, j.ext.min_y); - let (gx1, gy1) = key(j.ext.max_x, j.ext.max_y); - let cells = ((gx1 - gx0 + 1) as usize).saturating_mul((gy1 - gy0 + 1) as usize); - projected_refs = projected_refs.saturating_add(cells); - if projected_refs > MAX_GRID_REFERENCES { - use_grid = false; - break; - } - } - } - - let mut grid: std::collections::HashMap<(i64, i64), Vec> = std::collections::HashMap::new(); - if use_grid { - for (idx, j) in join_geoms.iter().enumerate() { - let (gx0, gy0) = key(j.ext.min_x, j.ext.min_y); - let (gx1, gy1) = key(j.ext.max_x, j.ext.max_y); - for gx in gx0..=gx1 { - for gy in gy0..=gy1 { - grid.entry((gx, gy)).or_default().push(idx); - } - } - } - } - + let index = ExtentIndex::new(join_geoms.iter().map(|j| j.ext)); // The join features whose extents may meet `ext`, in join order. - let candidates = |ext: &Extent| -> Vec { - if !use_grid { - return (0..join_geoms.len()).collect(); - } - if !ext.intersects(&total_ext) { - return Vec::new(); - } - let (tx0, ty0) = key(ext.min_x, ext.min_y); - let (tx1, ty1) = key(ext.max_x, ext.max_y); - let mut found = Vec::new(); - for gx in tx0..=tx1 { - for gy in ty0..=ty1 { - if let Some(indices) = grid.get(&(gx, gy)) { - found.extend_from_slice(indices); - } - } - } - found.sort_unstable(); - found.dedup(); - found - }; + let candidates = |ext: &Extent| index.candidates(ext); let mut total_matches = 0usize; let mut written_features = 0usize; @@ -1712,8 +1831,8 @@ fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result { // (A grid of coincident features has cells of no size: the first search // reaches at least `MIN_REACH_M`.) _ => { - let mut reach = if use_grid { - (cell_w.max(cell_h).to_radians() * rustgis_geom::EARTH_RADIUS_M).max(MIN_REACH_M) + let mut reach = if let Some(cell) = index.cell_size() { + (cell.to_radians() * rustgis_geom::EARTH_RADIUS_M).max(MIN_REACH_M) } else { f64::INFINITY }; @@ -3013,6 +3132,9 @@ mod tests { serde_json::json!({"target": "polys", "join": "pts", "match": "Within a distance", "distance": 1e308, "units": "Feet", "merge": "[{\"field\": 7}]"}), serde_json::json!({"target": "pts", "join": "polys", "match": "Largest overlap", "merge": [{"field": "k", "rule": "sum"}, null, "x"]}), serde_json::json!({"target": "polys", "join": "polys", "match": "Have their center in", "merge": [{"field": "k", "rule": "join", "delimiter": null}], "distanceField": ""}), + serde_json::json!({"input": "polys", "overlay": "polys"}), + serde_json::json!({"input": "pts", "overlay": "pts", "output": ""}), + serde_json::json!({"input": 3, "overlay": ["polys"]}), ]; for t in TOOLS { for p in &junk { @@ -3023,6 +3145,92 @@ mod tests { } } + /// Parcels and zones for Intersect: `p1` is a 2° square that zone `z1` covers the right half of; + /// `z2` overlaps `z1` but only touches `p1` along its right edge, and `p2` along its left. + /// Both layers have an `id` field. + fn parcels_and_zones() -> Map { + let text = |t: &str| Value::Text(t.into()); + let fields = |second: &str| vec![Field::new("id", FieldType::Text), Field::new(second, FieldType::Text)]; + let mut parcels = Dataset::new("parcels", GeometryType::Polygon, fields("use")); + parcels.push(Some(sq(0.0, 0.0, 2.0)), vec![text("p1"), text("home")]); + parcels.push(Some(sq(4.0, 0.0, 1.0)), vec![text("p2"), text("shop")]); + let mut zones = Dataset::new("zones", GeometryType::Polygon, fields("zone")); + zones.push(Some(sq(1.0, 0.0, 2.0)), vec![text("z1"), text("R1")]); + zones.push(Some(sq(2.0, 0.0, 2.0)), vec![text("z2"), text("C2")]); + Map(HashMap::from([("parcels".to_string(), parcels), ("zones".to_string(), zones)])) + } + + #[test] + fn intersect_keeps_the_overlap_with_both_layers_attributes() { + let m = parcels_and_zones(); + let out = layer(crate::run("intersect", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + assert_eq!(out.name, "parcels_Intersect"); + assert_eq!(out.geometry_type, GeometryType::Polygon); + let names: Vec<&str> = out.fields.iter().map(|f| f.name.as_str()).collect(); + assert_eq!(names, ["id", "use", "id_2", "zone"], "the zones' id is taken"); + // Only p1 and z1 overlap; touching along an edge makes no feature. + assert_eq!(out.features.len(), 1); + let f = &out.features[0]; + let text = |v: &[&str]| v.iter().map(|t| Value::Text((*t).into())).collect::>(); + assert_eq!(f.attrs, text(&["p1", "home", "z1", "R1"])); + let e = f.geometry.as_ref().unwrap().extent(); + for (got, want) in [(e.min_x, 1.0), (e.max_x, 2.0), (e.min_y, 0.0), (e.max_y, 2.0)] { + assert!((got - want).abs() < 1e-9, "{e:?}"); + } + } + + #[test] + fn intersect_splits_points_by_each_polygon_they_fall_in() { + let m = inputs(); + let out = layer(crate::run("intersect", &m, &serde_json::json!({"input": "pts", "overlay": "polys"})).unwrap()); + assert_eq!(out.geometry_type, GeometryType::Point); + // (0.5, 0.5) is in the first square only, (1.5, 1.5) in both overlapping ones, (50, 50) in none. + let rows: Vec<(f64, &str)> = out + .features + .iter() + .map(|f| (f.attrs[0].as_f64().unwrap(), if let Value::Text(t) = &f.attrs[1] { t.as_str() } else { "" })) + .collect(); + assert_eq!(rows, [(1.0, "a"), (2.0, "a"), (2.0, "a")]); + } + + #[test] + fn intersect_needs_an_overlay_of_polygons_and_an_input_with_geometry() { + let mut m = inputs(); + let e = crate::run("intersect", &m, &serde_json::json!({"input": "polys", "overlay": "pts"})).unwrap_err(); + assert!(e.contains("not a polygon layer"), "{e}"); + m.0.insert("table".into(), Dataset::new("table", GeometryType::None, vec![])); + let e = crate::run("intersect", &m, &serde_json::json!({"input": "table", "overlay": "polys"})).unwrap_err(); + assert!(e.contains("no geometry"), "{e}"); + } + + /// Layers with a selection in each, as the Tools pane passes them. + struct Selected(Map, HashMap>); + impl Inputs for Selected { + fn dataset(&self, layer: &str) -> Option<&Dataset> { + self.0.dataset(layer) + } + fn selection(&self, layer: &str) -> &[usize] { + self.1.get(layer).map_or(&[], Vec::as_slice) + } + } + + #[test] + fn intersect_honours_both_layers_selections() { + let mut parcels_and_zones = parcels_and_zones(); + // A zone covering p2, which only the selection would leave out. + let z3 = Some(sq(4.0, 0.0, 1.0)); + parcels_and_zones.0.get_mut("zones").unwrap().push(z3, vec![Value::Text("z3".into()), Value::Null]); + let m = Selected(parcels_and_zones, HashMap::from([("parcels".to_string(), vec![0]), ("zones".to_string(), vec![1, 2])])); + let out = layer(crate::run("intersect", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + // p1 meets only z1 (not selected); p2 meets z3 but isn't selected. + assert!(out.features.is_empty(), "{:?}", out.features); + // Selecting p2 too gives its part in z3, and only that. + let m = Selected(m.0, HashMap::from([("parcels".to_string(), vec![0, 1]), ("zones".to_string(), vec![1, 2])])); + let out = layer(crate::run("intersect", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + let ids: Vec<&Value> = out.features.iter().flat_map(|f| [&f.attrs[0], &f.attrs[2]]).collect(); + assert_eq!(ids, [&Value::Text("p2".into()), &Value::Text("z3".into())]); + } + #[test] fn select_by_location_invert_excludes_features_without_geometry() { let mut polys = Dataset::new("polys", GeometryType::Polygon, vec![]); diff --git a/crates/analysis/tests/reference.rs b/crates/analysis/tests/reference.rs index 5b28b96f8..9158eb002 100644 --- a/crates/analysis/tests/reference.rs +++ b/crates/analysis/tests/reference.rs @@ -65,15 +65,31 @@ fn json_text(v: &Json) -> Option { } } -/// The features of a layer by key: the `key` field, or one entry `"all"` when there is no key. -fn by_key(ds: &Dataset, key: Option<&str>) -> BTreeMap> { +/// A feature's key: its key fields' values as text, joined by `|` as `generate.py` writes a key of +/// several fields ("P1|Z1"); `"all"` when there are no key fields. +fn key_text(ds: &Dataset, row: usize, fields: Option<&[String]>) -> String { + match fields { + Some(fields) => fields.iter().map(|f| text(ds, row, f).unwrap_or_default()).collect::>().join("|"), + None => "all".into(), + } +} + +/// The fields a case's output is keyed by: its `key`, one field or a list; `id` without one; none +/// when it is null. +fn key_fields(expect: &Json) -> Option> { + match expect.get("key") { + Some(Json::Null) => None, + Some(Json::Array(fields)) => Some(fields.iter().filter_map(|f| f.as_str().map(String::from)).collect()), + Some(Json::String(field)) => Some(vec![field.clone()]), + _ => Some(vec!["id".into()]), + } +} + +/// The features of a layer by key (see [`key_text`]). +fn by_key(ds: &Dataset, key: Option<&[String]>) -> BTreeMap> { let mut out: BTreeMap> = BTreeMap::new(); for (row, f) in ds.features.iter().enumerate() { - let k = match key { - Some(key) => text(ds, row, key).unwrap_or_default(), - None => "all".into(), - }; - out.entry(k).or_default().extend(f.geometry.clone()); + out.entry(key_text(ds, row, key)).or_default().extend(f.geometry.clone()); } out } @@ -146,9 +162,11 @@ fn check(case: &Json, layers: &Layers) -> Vec { if let Some(path) = expect["features"].as_str() { let Some(ds) = layer_of(result.output) else { return vec!["no layer came out".into()] }; - let key = if expect.get("key").is_some_and(Json::is_null) { None } else { Some(expect["key"].as_str().unwrap_or("id")) }; - let want = by_key(&read_layer(&dir().join(path), "expected"), key.map(|_| "id")); - let got = by_key(&ds, key); + let key = key_fields(expect); + // The reference writes every key, of however many fields, as its `id`. + let id = ["id".to_string()]; + let want = by_key(&read_layer(&dir().join(path), "expected"), key.as_ref().map(|_| id.as_slice())); + let got = by_key(&ds, key.as_deref()); for (k, w) in &want { match got.get(k) { None => problems.push(format!("{k}: missing (the reference has it)")), @@ -161,10 +179,7 @@ fn check(case: &Json, layers: &Layers) -> Vec { if let Some(counts) = expect["counts"].as_object() { for (row, _) in ds.features.iter().enumerate() { // The same key as `by_key` gives the feature. - let k = match key { - Some(key) => text(&ds, row, key).unwrap_or_default(), - None => "all".into(), - }; + let k = key_text(&ds, row, key.as_deref()); let (n, want) = (text(&ds, row, "count"), counts.get(&k).and_then(json_text)); if n != want { problems.push(format!("{k}: count {n:?}, the reference {want:?}")); @@ -179,7 +194,10 @@ fn check(case: &Json, layers: &Layers) -> Vec { None => return vec![format!("a selection in '{layer}', which isn't a reference layer")], }, other => match layer_of(other) { - Some(ds) => (0..ds.features.len()).filter_map(|r| text(&ds, r, "id")).collect(), + Some(ds) => { + let key = key_fields(expect); + (0..ds.features.len()).map(|r| key_text(&ds, r, key.as_deref())).collect() + } None => return vec!["neither a selection nor a layer came out".into()], }, }; diff --git a/crates/analysis/tests/reference/cases.json b/crates/analysis/tests/reference/cases.json index 868fce829..0adf2a891 100644 --- a/crates/analysis/tests/reference/cases.json +++ b/crates/analysis/tests/reference/cases.json @@ -54,6 +54,57 @@ "features": "expected/erase_parcels.geojson" } }, + { + "name": "intersect polygons", + "tool": "intersect", + "params": { + "input": "parcels", + "overlay": "zones" + }, + "expect": { + "features": "expected/intersect_parcels.geojson", + "key": [ + "id", + "id_2" + ] + } + }, + { + "name": "intersect lines", + "tool": "intersect", + "params": { + "input": "roads", + "overlay": "zones" + }, + "expect": { + "features": "expected/intersect_roads.geojson", + "key": [ + "id", + "id_2" + ] + } + }, + { + "name": "intersect points", + "tool": "intersect", + "params": { + "input": "points", + "overlay": "zones" + }, + "expect": { + "ids": [ + "b|Z1", + "d|Z4", + "e|Z1", + "e|Z4", + "j|Z5" + ], + "key": [ + "id", + "id_2" + ] + } + }, { "name": "dissolve by field", "tool": "dissolve", @@ -151,7 +202,8 @@ "r1", "r2", "r4", - "r6" + "r6", + "r7" ] } }, @@ -187,7 +239,8 @@ "r1", "r2", "r3", - "r6" + "r6", + "r7" ] } }, @@ -428,6 +481,18 @@ [ "r6", "P8" + ], + [ + "r7", + "P1" + ], + [ + "r7", + "P2" + ], + [ + "r7", + "P6" ] ] } @@ -554,6 +619,11 @@ "r2", 28983.146568103246 ], + [ + "a", + "r7", + 29413.348290383856 + ], [ "b", "r2", @@ -564,11 +634,21 @@ "r1", 14084.263048419794 ], + [ + "b", + "r7", + 27864.913213951313 + ], [ "c", "r2", 6947.529600616946 ], + [ + "c", + "r7", + 13819.003210362416 + ], [ "c", "r1", @@ -584,6 +664,11 @@ "r3", 7906.554524361339 ], + [ + "e", + "r7", + 28851.844441780486 + ], [ "f", null, diff --git a/crates/analysis/tests/reference/expected/clip_roads.geojson b/crates/analysis/tests/reference/expected/clip_roads.geojson index e3ba16c2c..0c4b84a7d 100644 --- a/crates/analysis/tests/reference/expected/clip_roads.geojson +++ b/crates/analysis/tests/reference/expected/clip_roads.geojson @@ -76,6 +76,29 @@ ] ] } + }, + { + "type": "Feature", + "properties": { + "id": "r7" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.875, + 46.37499999999999 + ], + [ + 11.125, + 46.37499999999999 + ] + ] + } } ] } diff --git a/crates/analysis/tests/reference/expected/intersect_parcels.geojson b/crates/analysis/tests/reference/expected/intersect_parcels.geojson new file mode 100644 index 000000000..fbc23127e --- /dev/null +++ b/crates/analysis/tests/reference/expected/intersect_parcels.geojson @@ -0,0 +1,277 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "P1|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.250000000000002, + 45.999999999999986 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 10.500000000000002, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 11.0, + 46.37499999999999 + ], + [ + 11.0, + 45.999999999999986 + ], + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.500000000000002, + 46.37499999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2|Z6" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.0, + 46.50000000000001 + ], + [ + 11.0, + 46.12500000000001 + ], + [ + 10.875, + 46.12500000000001 + ], + [ + 10.875, + 46.50000000000001 + ], + [ + 11.0, + 46.50000000000001 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 46.37499999999999 + ], + [ + 11.500000000000002, + 46.37499999999999 + ], + [ + 11.500000000000002, + 46.25 + ], + [ + 11.500000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3|Z6" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 46.625 + ], + [ + 11.375, + 46.625 + ], + [ + 11.375, + 46.12500000000001 + ], + [ + 11.250000000000002, + 46.12500000000001 + ], + [ + 11.250000000000002, + 46.625 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P4|Z2" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.0, + 47.0 + ], + [ + 10.0, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.0 + ], + [ + 10.0, + 47.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P6|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.75, + 46.25 + ], + [ + 10.250000000000002, + 46.25 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.75, + 46.37499999999999 + ], + [ + 10.75, + 46.25 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P8|Z5" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 47.5 + ], + [ + 11.500000000000002, + 47.2505900350133 + ], + [ + 11.500000000000002, + 47.12499999999999 + ], + [ + 11.125, + 47.12499999999999 + ], + [ + 11.125, + 47.2505900350133 + ], + [ + 11.250000000000002, + 47.5 + ] + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/expected/intersect_roads.geojson b/crates/analysis/tests/reference/expected/intersect_roads.geojson new file mode 100644 index 000000000..cece9a3ad --- /dev/null +++ b/crates/analysis/tests/reference/expected/intersect_roads.geojson @@ -0,0 +1,138 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "r1|Z1" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.250000000000002, + 46.25 + ], + [ + 11.125, + 46.25 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r1|Z6" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.875, + 46.25 + ], + [ + 11.125, + 46.25 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r2|Z1" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.500000000000002, + 46.0625 + ], + [ + 10.500000000000002, + 46.37499999999999 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r3|Z4" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 11.5625, + 46.31250000000001 + ], + [ + 11.6875, + 46.4375 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r6|Z5" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 11.206763423774959, + 47.12499999999999 + ], + [ + 11.500000000000002, + 47.34499451815567 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r7|Z1" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 11.125, + 46.37499999999999 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r7|Z6" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.875, + 46.37499999999999 + ], + [ + 11.125, + 46.37499999999999 + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/generate.py b/crates/analysis/tests/reference/generate.py index 20b574d17..0593b5174 100644 --- a/crates/analysis/tests/reference/generate.py +++ b/crates/analysis/tests/reference/generate.py @@ -75,6 +75,8 @@ def ring(*pts): ("r4", LineString([(12, 46.5), (12.5, 46.5)])), ("r5", LineString([(15, 49), (15.5, 49.5)])), ("r6", LineString([(10.875, 46.875), (11.875, 47.625)])), + # Along Z1's top edge, and on into Z6. + ("r7", LineString([(10.25, 46.375), (11.125, 46.375)])), ] STATIONS = [ ("s1", Point(10.75, 46.0)), @@ -230,6 +232,26 @@ def case(name, tool, params, expect): case("erase", "erase", {"input": "parcels", "erase": "zones"}, {"features": overlay("erase_parcels", parcels, lambda g: polygonal(g.difference(zones_m)))}) + # Intersect: each input feature's part in each zone, keyed by both ids ("P1|Z1"), in input + # order, then zone order. Parts that only touch (a shared edge or corner) are left out. The + # zones' `id` is taken by the parcels', so it comes out as `id_2`. + def pieces(name, rows, part): + out = [] + for i, g in rows: + for zi, z in ZONES: + r = part(m(g).intersection(m(z))) + if r is not None and not r.is_empty: + out.append((f"{i}|{zi}", ll(r))) + return write_expected(name, out) + + both = ["id", "id_2"] + case("intersect polygons", "intersect", {"input": "parcels", "overlay": "zones"}, + {"features": pieces("intersect_parcels", parcels, polygonal), "key": both}) + case("intersect lines", "intersect", {"input": "roads", "overlay": "zones"}, + {"features": pieces("intersect_roads", ROADS, linear), "key": both}) + case("intersect points", "intersect", {"input": "points", "overlay": "zones"}, + {"ids": [f"{i}|{zi}" for i, g in POINTS for zi, z in ZONES if m(z).intersects(m(g))], "key": both}) + # Dissolve: one feature per code, the union of its parcels. codes = sorted({c for _, c, _ in PARCELS}) rows = [(c, ll(unary_union([m(g) for _, k, g in PARCELS if k == c]))) for c in codes] diff --git a/crates/analysis/tests/reference/layers/roads.geojson b/crates/analysis/tests/reference/layers/roads.geojson index c7d2d2056..b08b67d2a 100644 --- a/crates/analysis/tests/reference/layers/roads.geojson +++ b/crates/analysis/tests/reference/layers/roads.geojson @@ -114,6 +114,25 @@ ] ] } + }, + { + "type": "Feature", + "properties": { + "id": "r7" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.25, + 46.375 + ], + [ + 11.125, + 46.375 + ] + ] + } } ] } diff --git a/crates/assistant/skills/geoprocessing/SKILL.md b/crates/assistant/skills/geoprocessing/SKILL.md index e44879306..df374da38 100644 --- a/crates/assistant/skills/geoprocessing/SKILL.md +++ b/crates/assistant/skills/geoprocessing/SKILL.md @@ -11,7 +11,8 @@ Based on the RustGIS user guide, Analysis › Geoprocessing tools 1. **Find the tool:** `search_tools` with words ("buffer", "overlay", "count points"), then `describe_tool` for its parameters. Never guess parameter names. The built-in tools are: - proximity: `buffer` (`input`, `distance`, `units`, `dissolve`) - - overlay: `clip` (`input`, `clip`), `erase` (`input`, `erase`) + - overlay: `clip` (`input`, `clip`), `erase` (`input`, `erase`), `intersect` (`input`, + `overlay`: the parts in each overlay polygon, with both layers' attributes) - generalization: `dissolve` (`input`, `field`), `simplify` (`input`, `tolerance`) - geometry: `convexHull`, `featureToPoint`, `xyTableToPoint` (`input`, `x`, `y`, `crs`), `calculateGeometry` diff --git a/crates/engine/src/cmd/workflows.rs b/crates/engine/src/cmd/workflows.rs index 5532b05ae..5f6800d40 100644 --- a/crates/engine/src/cmd/workflows.rs +++ b/crates/engine/src/cmd/workflows.rs @@ -185,9 +185,18 @@ fn gallery(_: &mut Session, p: &Value) -> CmdResult { pub fn gallery_list(kind: &str) -> CmdResult { let ids: &[&str] = match kind { "ready" => &["buffer", "dissolve", "countPoints", "clip", "summaryStatistics", "selectByLocation"], - "feature" => { - &["buffer", "clip", "erase", "dissolve", "convexHull", "featureToPoint", "simplify", "countPoints", "calculateGeometry"] - } + "feature" => &[ + "buffer", + "clip", + "erase", + "intersect", + "dissolve", + "convexHull", + "featureToPoint", + "simplify", + "countPoints", + "calculateGeometry", + ], "raster" => &[], k => return Err(format!("unknown gallery '{k}'")), }; diff --git a/crates/ui-egui/src/icons.rs b/crates/ui-egui/src/icons.rs index 3bffb4393..1f1bc3a9e 100644 --- a/crates/ui-egui/src/icons.rs +++ b/crates/ui-egui/src/icons.rs @@ -502,6 +502,12 @@ pub fn paint(painter: &Painter, r: Rect, name: &str, ink: Color32, accent: Color p.polygon_shape(&[(2.0, 4.0), (16.0, 3.0), (20.0, 14.0), (9.0, 20.0)], 110); p.circle(16.0, 16.0, 6.0, Some(Color32::WHITE), Some((1.3, red))); } + "intersect" => { + // Two overlapping polygons, only their shared part filled. + p.rect(9.0, 9.0, 6.0, 6.0, Some(accent), None, 0.0); + p.rect(3.0, 3.0, 12.0, 12.0, None, Some((1.3, ink)), 0.0); + p.rect(9.0, 9.0, 12.0, 12.0, None, Some((1.3, ink)), 0.0); + } "dissolve" => { p.polygon_shape(&[(2.0, 6.0), (12.0, 3.0), (22.0, 6.0), (20.0, 19.0), (4.0, 19.0)], 110); p.line(&[(12.0, 3.5), (12.0, 18.5)], 1.0, ink.gamma_multiply(0.35)); @@ -1202,6 +1208,7 @@ pub const ALL: &[&str] = &[ "buffer", "clip", "erase", + "intersect", "dissolve", "hull", "centroid", diff --git a/crates/ui-egui/src/panes.rs b/crates/ui-egui/src/panes.rs index 8dc0eee03..9ae7cf119 100644 --- a/crates/ui-egui/src/panes.rs +++ b/crates/ui-egui/src/panes.rs @@ -3389,6 +3389,7 @@ pub fn tool_icon(id: &str) -> &'static str { "buffer" => "buffer", "clip" => "clip", "erase" => "erase", + "intersect" => "intersect", "dissolve" => "dissolve", "convexHull" => "hull", "featureToPoint" => "centroid", diff --git a/crates/ui-egui/src/ribbon.rs b/crates/ui-egui/src/ribbon.rs index 60757baf6..01eabce2e 100644 --- a/crates/ui-egui/src/ribbon.rs +++ b/crates/ui-egui/src/ribbon.rs @@ -1499,6 +1499,7 @@ const GALLERY: &[(&str, &str, &str)] = &[ ("buffer", "buffer", "Buffer"), ("clip", "clip", "Clip"), ("erase", "erase", "Erase"), + ("intersect", "intersect", "Intersect"), ("dissolve", "dissolve", "Dissolve"), ("countPoints", "count", "Summarize\nWithin"), ("convexHull", "hull", "Bounding\nGeometry"), diff --git a/docs/comparison/index.md b/docs/comparison/index.md index 543647bf4..63b4bdafa 100644 --- a/docs/comparison/index.md +++ b/docs/comparison/index.md @@ -78,7 +78,7 @@ data and AI providers, send requests over the network. See |---|---|---|---|---|---| | **Geoprocessing tools** | Built-in tools plus the optional **1,000+** GeoLibre/Whitebox suite (the available list comes from the runner) | **1,000+** (GeoLibre and Whitebox) | **1,000+** (native, GDAL, GRASS, SAGA) | **Most mature and complete**, plus paid extensions | A useful subset, credit-metered | | **Where analysis runs** | Your machine; in the browser on WebAssembly | Your browser (WASM); optional Python sidecar on desktop | Your machine | Your machine, or a server | Vendor cloud | -| **Vector analysis** | Geodesic buffer, clip, erase, dissolve, spatial join, summarize within, statistics, selection by location, A5 grids and aggregation; optional GeoLibre vector tools | Buffer, overlay, dissolve, joins, topology checks | Comprehensive | Comprehensive | Common tools | +| **Vector analysis** | Geodesic buffer, clip, erase, intersect, dissolve, spatial join, summarize within, statistics, selection by location, A5 grids and aggregation; optional GeoLibre vector tools | Buffer, overlay, dissolve, joins, topology checks | Comprehensive | Comprehensive | Common tools | | **Raster analysis** | Raster functions, band math, indices, change detection, classification, viewshed, IDW and ordinary kriging; time-series summaries and trends; optional GeoLibre terrain and hydrology tools | Hillshade, slope, zonal stats, raster calculator, reclassify | Comprehensive (GDAL/GRASS/SAGA) | Comprehensive; Spatial Analyst | Limited | | **Query language** | SQL `WHERE` clauses for selection and definition queries; DuckDB Spatial queries over files, web addresses and databases in an SQL pane (desktop) | DuckDB Spatial, PGlite/PostGIS, Sedona, all in the browser | Virtual layers, PostGIS | SQL against geodatabases | — | | **Model / batch chaining** | ModelBuilder (a node canvas; models from history), tasks, scripts | Visual Model Builder, batch runner | Graphical Model Designer | ModelBuilder | — | @@ -142,8 +142,8 @@ Stated plainly, so the tables above are worth trusting: - **Coordinate systems.** Vector data is stored in longitude/latitude; there is no full PROJ. - **Editing depth.** Vertex move, insert and delete are supported, but not versioned or multi-user database editing. Edits are kept in the project and exported, not written back to databases. -- **Specialized analysis.** Built-in spatial joins and shortest paths on line networks are - available, but full network datasets, overlay union/intersect tools and advanced spatial +- **Specialized analysis.** Built-in spatial joins, intersect and shortest paths on line networks + are available, but full network datasets, overlay union and identity tools and advanced spatial statistics depend on the external suite or another GIS. - **Python.** The [`rustgis` package](../python.md) runs commands headlessly, but there is no Python console inside the app. Python process plugins are supported on desktop and CLI, without a sandbox. diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 9e24ca0c9..460df9b3b 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -78,7 +78,7 @@ written some other way (a token typed into a SQL query, say) can still be there: before sharing it. A point buffered where the map was clicked is copied as the command that buffered it, with the point. -There are 20 built-in tools; GeoLibre and plugin tools join them when available. +There are 21 built-in tools; GeoLibre and plugin tools join them when available. Use these ids in `gp.run`'s `tool` parameter (`gp.tools` lists the available tools and their parameters): | Tool | Id | Toolbox | @@ -86,6 +86,7 @@ Use these ids in `gp.run`'s `tool` parameter (`gp.tools` lists the available too | Buffer | `buffer` | Analysis › Proximity | | Clip | `clip` | Analysis › Extract | | Erase | `erase` | Analysis › Overlay | +| Intersect | `intersect` | Analysis › Overlay | | Spatial Join | `spatialJoin` | Analysis › Overlay | | Dissolve | `dissolve` | Data Management › Generalization | | Minimum Bounding Geometry | `convexHull` | Data Management › Features | @@ -152,6 +153,23 @@ lines and polygons contribute their bounding-box centres. At least two distinct numeric values are needed. It adds **two separate Float32 GeoTIFF layers**: the prediction, named by `output`, and `_variance`. +**Intersect** (`intersect`) splits `input` features by the polygons of `overlay`: each part of an +input feature inside an overlay polygon becomes a feature with the attributes of both, the input's +first. A parcel two flood zones cross gives two features, one for each zone, and a point in two +overlapping polygons gives two. Input features outside every overlay polygon are left out. An +overlay polygon's edge counts as inside it, as for Clip: a point on the edge and a stretch of line +along it are kept, but polygons that only share an edge or a corner, and a line that only touches +one, give nothing. The input can be points, lines or polygons and keeps its geometry type; the +overlay must be polygons. An overlay field whose name the input already has gets a suffix (`id_2`). +Overlapping features within one layer are not split against each other. Like Clip and Erase, the +overlay works on the map's Web Mercator coordinates, so the parts' edges are straight on the map. +Selections in both layers are honoured. To measure each part, run Calculate Geometry Attributes on +the result: + +```text +gp.run={"tool":"intersect","params":{"input":"parcels","overlay":"floodplain","output":"parcels_in_floodplain"}} +``` + **Spatial Join** (`spatialJoin`, also `gp.spatialJoin`) copies attributes from `join` features to `target` features. Layers may have different source coordinate systems: their geometry is stored in longitude/latitude before the join. `match` is one of: From 3a6edb0ceeb88c1523bb7b84e5ab1120d65b1282 Mon Sep 17 00:00:00 2001 From: mavrick Date: Mon, 26 Oct 2026 14:31:32 -0500 Subject: [PATCH 2/5] feat(analysis): add Identity and Union overlay tools Identity keeps every part of the input: each feature's parts in the overlay polygons, with both layers' attributes, then its part outside them all. Union, for two polygon layers, also adds each overlay polygon's part outside every input polygon. Both follow Intersect's rules for fields, edges and selections. Intersect, Identity and Union are now one overlay function, and Clip's helper becomes part(g, mask, inside), which also gives the part outside a mask. GEOS reference cases cover identity of polygons, lines and points, and union. Part of #959. Co-Authored-By: Claude Opus 5.5 --- GEOLIBRE_PARITY.md | 4 +- README.md | 2 +- ROADMAP.md | 13 +- crates/analysis/src/tools.rs | 207 ++++- crates/analysis/tests/reference/cases.json | 72 ++ .../expected/identity_parcels.geojson | 597 ++++++++++++ .../reference/expected/identity_roads.geojson | 245 +++++ .../reference/expected/union_parcels.geojson | 857 ++++++++++++++++++ crates/analysis/tests/reference/generate.py | 29 +- .../assistant/skills/geoprocessing/SKILL.md | 4 +- crates/engine/src/cmd/workflows.rs | 2 + docs/comparison/index.md | 4 +- docs/user-guide/analysis.md | 12 +- 13 files changed, 2010 insertions(+), 38 deletions(-) create mode 100644 crates/analysis/tests/reference/expected/identity_parcels.geojson create mode 100644 crates/analysis/tests/reference/expected/identity_roads.geojson create mode 100644 crates/analysis/tests/reference/expected/union_parcels.geojson diff --git a/GEOLIBRE_PARITY.md b/GEOLIBRE_PARITY.md index f3010b864..43ab92678 100644 --- a/GEOLIBRE_PARITY.md +++ b/GEOLIBRE_PARITY.md @@ -159,8 +159,8 @@ _GeoLibre: docs/features.md › Attribute data and expressions; docs/arcgis-edit _GeoLibre: docs/features.md › Processing and analysis; docs/user-guide/processing.md, model-builder.md, segmentation.md, obia.md_ -- [ ] **Spatial join, intersect, union**. Partial: spatial join, clip, erase and intersect are - built in; union, identity and near are missing (_Roadmap: solid ground_). +- [ ] **Spatial join, intersect, union**. Partial: spatial join, clip, erase, intersect, identity + and union are built in; near is missing (_Roadmap: solid ground_). - [ ] **More built-in vector tools**: attribute join, merge, extract vertices, points along lines, fix/check geometry, random extract, buffer options. Partial (_Roadmap: solid ground_). - [ ] **Routing services**: isochrones, service areas, OD matrix, multi-stop routing, diff --git a/README.md b/README.md index a2e432fca..1e4b02da3 100644 --- a/README.md +++ b/README.md @@ -99,7 +99,7 @@ Every screenshot is RustGIS itself, rendered offscreen by `cargo run -p rustgis- | **Web tiles** | XYZ/TMS/WMS raster tiles, Mapbox Vector Tiles and MapLibre styles (sprites, patterns, expressions), cached and fetched in the background | | **Cartography** | Single symbol, unique values, graduated colours, proportional symbols; label classes with SQL queries, straight and curved line placement, installed font families, scale ranges and annotation; stacked and dashed line strokes, symbol levels, reference scales and polygons without fill; masking, display filters, aggregation bins, 2.5D extrusion; graphics and dynamic text | | **Editing** | Create, move, rotate, scale, merge, split and explode features; snapping, templates, elevation, traverse, topology rules; domains, subtypes, attribute rules, relationships, archiving; undo for everything | -| **Analysis** | Selection by attributes (SQL) and location; Buffer, Clip, Erase, Intersect, Dissolve, Summarize Within, Summary Statistics, Spatial Join, Zonal Statistics, A5 grids and a geometry calculator; viewshed, line of sight, IDW, ordinary kriging, shortest path, suitability modelling; ModelBuilder on a node canvas and scripts with command, parameter and value completion; 1,000+ [GeoLibre](https://github.com/opengeos/geolibre-rust) and Whitebox raster, vector, LiDAR, hydrology and terrain tools | +| **Analysis** | Selection by attributes (SQL) and location; Buffer, Clip, Erase, Intersect, Identity, Union, Dissolve, Summarize Within, Summary Statistics, Spatial Join, Zonal Statistics, A5 grids and a geometry calculator; viewshed, line of sight, IDW, ordinary kriging, shortest path, suitability modelling; ModelBuilder on a node canvas and scripts with command, parameter and value completion; 1,000+ [GeoLibre](https://github.com/opengeos/geolibre-rust) and Whitebox raster, vector, LiDAR, hydrology and terrain tools | | **Imagery** | Ortho and reality mapping, pixel classification (k-means, minimum distance, a trainable neural network), mensuration, hosted imagery tiles | | **Layouts and reports** | Page layouts with map frames, legends, scale bars, grids and map series; publication sheets and figures in one command, format panes, guides, snapping and element locks; reports with groups, statistics, charts and subreports; presentations (HTML slideshow, PDF); notebooks (Jupyter in and out) | | **Sharing** | PNG export at any DPI; PDF; web map pages, OGC API Features collections, MapLibre styles, glTF scenes, locators; project, map and mobile map packages; offline areas with sync | diff --git a/ROADMAP.md b/ROADMAP.md index 2ac014de1..90875e9c2 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -60,8 +60,7 @@ Correct, robust and recoverable, with the foundations other themes stand on. administrator can fix for a managed install (no AI Assistant, plugins or web services). - Extend the OS keychain secret store to service connections (database passwords, S3 keys, AI provider keys and the ChatGPT sign-in's tokens are already kept there). -- Overlay tools beside Clip, Erase and Intersect: union, identity and near - ([#959](https://github.com/opengeos/RustGIS/issues/959)). +- A Near tool beside the overlay tools ([#959](https://github.com/opengeos/RustGIS/issues/959)). - Spatial Join by distance, closest, centre in and largest overlap, with merge rules (sum, mean, first…) for one-to-one joins ([#851](https://github.com/opengeos/RustGIS/issues/851)). - Vector tools, each writing a new layer: join field (a one-off attribute join), merge (append @@ -220,11 +219,11 @@ Open questions rather than plans; discuss before building. templates, elevation, topology; domains, subtypes, attribute rules, contingent values, relationships, archiving. - Analysis: selection by attributes and location; polygon, lasso and radius on the map; overlay - (clip, erase, intersect), proximity (buffer) and summary tools, spatial join, zonal statistics, - A5 grids, kriging and Calculate Geometry with metric or imperial units; guided workflows - (viewshed, line of sight, IDW, shortest path, suitability); ModelBuilder on a node canvas, - scripts with command, parameter and value completion, and history copied as code; the GeoLibre - and Whitebox suite (1,000+ tools) natively and in the browser. + (clip, erase, intersect, identity, union), proximity (buffer) and summary tools, spatial join, + zonal statistics, A5 grids, kriging and Calculate Geometry with metric or imperial units; guided + workflows (viewshed, line of sight, IDW, shortest path, suitability); ModelBuilder on a node + canvas, scripts with command, parameter and value completion, and history copied as code; the + GeoLibre and Whitebox suite (1,000+ tools) natively and in the browser. - Layouts (map frames, surrounds, map series, publication sheets and figures in one command, format panes, guides, snapping, element locking and PDF), reports, presentations and notebooks; the web app embeddable in a page (`?embed`). diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index f934da209..05d646588 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -70,6 +70,30 @@ pub static TOOLS: &[ToolSpec] = &[ ], run: intersect, }, + ToolSpec { + id: "identity", + name: "Identity", + toolbox: "Analysis › Overlay", + summary: "Splits input features by overlay polygons, keeping every part: those inside get the overlay's attributes.", + params: &[ + P { name: "input", label: "Input Features", kind: K::Layer, required: true, default: "", help: "" }, + P { name: "overlay", label: "Identity Features", kind: K::PolygonLayer, required: true, default: "", help: "" }, + P { name: "output", label: "Output Name", kind: K::Text, required: false, default: "", help: "" }, + ], + run: identity, + }, + ToolSpec { + id: "union", + name: "Union", + toolbox: "Analysis › Overlay", + summary: "Combines two polygon layers into every piece either covers, with the attributes of both where they overlap.", + params: &[ + P { name: "input", label: "Input Features", kind: K::PolygonLayer, required: true, default: "", help: "" }, + P { name: "overlay", label: "Union Features", kind: K::PolygonLayer, required: true, default: "", help: "" }, + P { name: "output", label: "Output Name", kind: K::Text, required: false, default: "", help: "" }, + ], + run: union, + }, ToolSpec { id: "dissolve", name: "Dissolve", @@ -692,30 +716,33 @@ fn union_of(inp: &dyn Inputs, p: &Json, k: &str) -> Result, St } /// The part of a feature inside `mask` (polygons in Web Mercator, oriented as `union_all` leaves -/// them): a polygon's overlap with it, a line's stretches in it, the points in it or on its edge. -/// `None` when nothing is left, or only a shared edge or corner. -fn part_in(g: &Geometry, mask: &MultiPolygon) -> Option { +/// them), or outside it when `inside` is false. Inside, a polygon's overlap with it, a line's +/// stretches in it, the points in it or on its edge; outside, the rest. `None` when nothing is +/// left, or only a shared edge or corner. +fn part(g: &Geometry, mask: &MultiPolygon, inside: bool) -> Option { match g { - Geometry::Polygon(_) => to_multipolygon(g).and_then(|m| from_multipolygon(&m.intersection(mask))), - Geometry::Line(_) => to_multiline(g).and_then(|m| from_line(&mask.clip(&m, false))), + Geometry::Polygon(_) => { + to_multipolygon(g).and_then(|m| from_multipolygon(&if inside { m.intersection(mask) } else { m.difference(mask) })) + } + Geometry::Line(_) => to_multiline(g).and_then(|m| from_line(&mask.clip(&m, !inside))), Geometry::Point(_) | Geometry::MultiPoint(_) => { let pts: &[Coord] = match g { Geometry::Point(c) => std::slice::from_ref(c), Geometry::MultiPoint(v) => v, _ => &[], }; - let inside: Vec = pts + let kept: Vec = pts .iter() .copied() .filter(|c| { let m = rustgis_geom::lonlat_to_merc(*c); - mask.intersects(&geo::Point::new(m.x, m.y)) + mask.intersects(&geo::Point::new(m.x, m.y)) == inside }) .collect(); - match inside.len() { + match kept.len() { 0 => None, - 1 => inside.first().map(|c| Geometry::Point(*c)), - _ => Some(Geometry::MultiPoint(inside)), + 1 => kept.first().map(|c| Geometry::Point(*c)), + _ => Some(Geometry::MultiPoint(kept)), } } } @@ -727,7 +754,7 @@ fn clip(inp: &dyn Inputs, p: &Json) -> Result { let mut out = copy_schema(ds, out_name(p, format!("{}_Clip", ds.name)), ds.geometry_type); for r in rows(inp, p, "input", ds) { let Some(f) = ds.features.get(r) else { continue }; - if let Some(c) = f.geometry.as_ref().and_then(|g| part_in(g, &mask)) { + if let Some(c) = f.geometry.as_ref().and_then(|g| part(g, &mask, true)) { out.push(Some(c), f.attrs.clone()); } } @@ -796,32 +823,107 @@ fn overlay_row(ds: &Dataset, a: Option<&[Value]>, ov: &Dataset, b: Option<&[Valu row } -/// Intersect: a feature for each part of an input feature inside an overlay polygon, with the -/// attributes of both, in input order and then overlay order. Parts that only touch are left out. -fn intersect(inp: &dyn Inputs, p: &Json) -> Result { +/// What an overlay tool keeps. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +enum Overlay { + /// The parts of input features inside overlay polygons. + Intersect, + /// Those, and the parts of input features outside every overlay polygon. + Identity, + /// Those, and the parts of overlay polygons outside every input polygon. + Union, +} + +impl Overlay { + /// The tool's name, for its messages and its output's default name. + fn name(self) -> &'static str { + match self { + Overlay::Intersect => "Intersect", + Overlay::Identity => "Identity", + Overlay::Union => "Union", + } + } +} + +/// Intersect, Identity and Union. Each input feature in turn gives its part in each overlay +/// polygon, in overlay order, with the attributes of both; then (Identity, Union) its part outside +/// them all, with empty overlay attributes. Union then adds each overlay polygon's part outside +/// every input polygon, with empty input attributes. Parts that only touch are left out, and +/// features of one layer that overlap each other are not split against each other. +fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result { let ds = input(inp, p, "input")?; if ds.geometry_type == GeometryType::None { return Err(format!("'{}' has no geometry", ds.name)); } + if op == Overlay::Union && ds.geometry_type != GeometryType::Polygon { + return Err(format!("'{}' is not a polygon layer", ds.name)); + } let (ov, zones) = zones(inp, p)?; let index = ExtentIndex::new(zones.iter().map(|z| z.ext)); - let mut out = overlay_schema(ds, ov, out_name(p, format!("{}_Intersect", ds.name)), ds.geometry_type); + let mut out = overlay_schema(ds, ov, out_name(p, format!("{}_{}", ds.name, op.name())), ds.geometry_type); + let zone_attrs = |z: &Zone| ov.features.get(z.row).map(|o| o.attrs.as_slice()); + // For Union, the input polygons that overlap each zone. + let mut covering: Vec>> = vec![Vec::new(); if op == Overlay::Union { zones.len() } else { 0 }]; for r in rows(inp, p, "input", ds) { let Some(f) = ds.features.get(r) else { continue }; let Some(g) = &f.geometry else { continue }; let ext = g.extent(); - for z in index.candidates(&ext).into_iter().filter_map(|i| zones.get(i)) { - if !ext.intersects(&z.ext) { - continue; + let mut hit: Vec = Vec::new(); + for zi in index.candidates(&ext) { + let Some(z) = zones.get(zi).filter(|z| ext.intersects(&z.ext)) else { continue }; + if let Some(inside) = part(g, &z.poly, true) { + out.push_next(Some(inside), overlay_row(ds, Some(&f.attrs), ov, zone_attrs(z))); + hit.push(zi); } - if let Some(part) = part_in(g, &z.poly) { - let attrs = overlay_row(ds, Some(&f.attrs), ov, ov.features.get(z.row).map(|o| o.attrs.as_slice())); - out.push_next(Some(part), attrs); + } + if op == Overlay::Intersect { + continue; + } + let rest = if hit.is_empty() { + Some(g.clone()) + } else { + part(g, &union_all(hit.iter().filter_map(|zi| zones.get(*zi)).flat_map(|z| z.poly.0.iter())), false) + }; + if let Some(rest) = rest { + out.push_next(Some(rest), overlay_row(ds, Some(&f.attrs), ov, None)); + } + if op == Overlay::Union + && let Some(poly) = to_multipolygon(g) + { + for zi in hit { + if let Some(c) = covering.get_mut(zi) { + c.push(poly.clone()); + } } } } + for (z, inputs) in zones.iter().zip(&covering) { + let rest = if inputs.is_empty() { + from_multipolygon(&z.poly) + } else { + from_multipolygon(&z.poly.difference(&union_all(inputs.iter().flat_map(|m| m.0.iter())))) + }; + if let Some(rest) = rest { + out.push_next(Some(rest), overlay_row(ds, None, ov, zone_attrs(z))); + } + } let n = out.features.len(); - result(out, format!("Intersect: {n} features")) + result(out, format!("{}: {n} features", op.name())) +} + +/// Intersect: see [`overlay`]. +fn intersect(inp: &dyn Inputs, p: &Json) -> Result { + overlay(inp, p, Overlay::Intersect) +} + +/// Identity: see [`overlay`]. +fn identity(inp: &dyn Inputs, p: &Json) -> Result { + overlay(inp, p, Overlay::Identity) +} + +/// Union: see [`overlay`]. +fn union(inp: &dyn Inputs, p: &Json) -> Result { + overlay(inp, p, Overlay::Union) } fn dissolve(inp: &dyn Inputs, p: &Json) -> Result { @@ -3203,6 +3305,49 @@ mod tests { assert!(e.contains("no geometry"), "{e}"); } + /// Each feature's `id` and the overlay's (`id_2`), "" for an empty one: "p1|z1", "p1|", "|z2". + fn overlay_keys(out: &Dataset) -> Vec { + let text = |v: Option<&Value>| if let Some(Value::Text(t)) = v { t.clone() } else { String::new() }; + out.features.iter().map(|f| format!("{}|{}", text(f.attrs.first()), text(f.attrs.get(2)))).collect() + } + + #[test] + fn identity_keeps_every_part_of_the_input() { + let m = parcels_and_zones(); + let out = layer(crate::run("identity", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + assert_eq!(out.name, "parcels_Identity"); + // p1's half in z1, then its other half; p2 only touches z2, so it is kept whole. + assert_eq!(overlay_keys(&out), ["p1|z1", "p1|", "p2|"]); + let rest = out.features[1].geometry.as_ref().unwrap().extent(); + assert!((rest.min_x - 0.0).abs() < 1e-9 && (rest.max_x - 1.0).abs() < 1e-9, "{rest:?}"); + assert_eq!(out.features[1].attrs[2..], [Value::Null, Value::Null]); + assert_eq!(out.features[2].geometry, m.0["parcels"].features[1].geometry); + } + + #[test] + fn identity_keeps_points_outside_every_polygon() { + let m = inputs(); + let out = layer(crate::run("identity", &m, &serde_json::json!({"input": "pts", "overlay": "polys"})).unwrap()); + let rows: Vec<(f64, Value)> = out.features.iter().map(|f| (f.attrs[0].as_f64().unwrap(), f.attrs[1].clone())).collect(); + let a = Value::Text("a".into()); + assert_eq!(rows, [(1.0, a.clone()), (2.0, a.clone()), (2.0, a), (3.0, Value::Null)]); + } + + #[test] + fn union_adds_the_overlay_left_uncovered() { + let m = parcels_and_zones(); + let out = layer(crate::run("union", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + assert_eq!(out.name, "parcels_Union"); + // z1 less p1's half of it; z2, which no parcel overlaps, whole (overlapping z1 still). + assert_eq!(overlay_keys(&out), ["p1|z1", "p1|", "p2|", "|z1", "|z2"]); + let z1_rest = out.features[3].geometry.as_ref().unwrap().extent(); + assert!((z1_rest.min_x - 2.0).abs() < 1e-9 && (z1_rest.max_x - 3.0).abs() < 1e-9, "{z1_rest:?}"); + assert_eq!(out.features[3].attrs[..2], [Value::Null, Value::Null]); + assert_eq!(out.features[4].attrs[3], Value::Text("C2".into())); + let e = crate::run("union", &inputs(), &serde_json::json!({"input": "pts", "overlay": "polys"})).unwrap_err(); + assert!(e.contains("not a polygon layer"), "{e}"); + } + /// Layers with a selection in each, as the Tools pane passes them. struct Selected(Map, HashMap>); impl Inputs for Selected { @@ -3231,6 +3376,24 @@ mod tests { assert_eq!(ids, [&Value::Text("p2".into()), &Value::Text("z3".into())]); } + #[test] + fn identity_and_union_honour_both_layers_selections() { + // Only z2 is selected: it touches p1 without overlapping it, so z1 doesn't cut p1. + let m = Selected(parcels_and_zones(), HashMap::from([("zones".to_string(), vec![1])])); + let out = layer(crate::run("identity", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + assert_eq!(overlay_keys(&out), ["p1|", "p2|"]); + assert_eq!(out.features[0].geometry, m.0.0["parcels"].features[0].geometry); + // Only p2 is selected: it touches z2 without overlapping it, so p1 doesn't cover z1. + let m = Selected(parcels_and_zones(), HashMap::from([("parcels".to_string(), vec![1])])); + let out = layer(crate::run("union", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + assert_eq!(overlay_keys(&out), ["p2|", "|z1", "|z2"]); + // z1 whole (back from Web Mercator, so not to the last bit). + let (got, want) = (out.features[1].geometry.as_ref().unwrap().extent(), sq(1.0, 0.0, 2.0).extent()); + for (g, w) in [(got.min_x, want.min_x), (got.max_x, want.max_x), (got.min_y, want.min_y), (got.max_y, want.max_y)] { + assert!((g - w).abs() < 1e-9, "{got:?}"); + } + } + #[test] fn select_by_location_invert_excludes_features_without_geometry() { let mut polys = Dataset::new("polys", GeometryType::Polygon, vec![]); diff --git a/crates/analysis/tests/reference/cases.json b/crates/analysis/tests/reference/cases.json index 0adf2a891..7a196cb00 100644 --- a/crates/analysis/tests/reference/cases.json +++ b/crates/analysis/tests/reference/cases.json @@ -105,6 +105,78 @@ ] } }, + { + "name": "identity polygons", + "tool": "identity", + "params": { + "input": "parcels", + "overlay": "zones" + }, + "expect": { + "features": "expected/identity_parcels.geojson", + "key": [ + "id", + "id_2" + ] + } + }, + { + "name": "identity lines", + "tool": "identity", + "params": { + "input": "roads", + "overlay": "zones" + }, + "expect": { + "features": "expected/identity_roads.geojson", + "key": [ + "id", + "id_2" + ] + } + }, + { + "name": "identity points", + "tool": "identity", + "params": { + "input": "points", + "overlay": "zones" + }, + "expect": { + "ids": [ + "a|", + "b|Z1", + "c|", + "d|Z4", + "e|Z1", + "e|Z4", + "f|", + "g|", + "h|", + "i|", + "j|Z5" + ], + "key": [ + "id", + "id_2" + ] + } + }, + { + "name": "union", + "tool": "union", + "params": { + "input": "parcels", + "overlay": "zones" + }, + "expect": { + "features": "expected/union_parcels.geojson", + "key": [ + "id", + "id_2" + ] + } + }, { "name": "dissolve by field", "tool": "dissolve", diff --git a/crates/analysis/tests/reference/expected/identity_parcels.geojson b/crates/analysis/tests/reference/expected/identity_parcels.geojson new file mode 100644 index 000000000..8be37e0f6 --- /dev/null +++ b/crates/analysis/tests/reference/expected/identity_parcels.geojson @@ -0,0 +1,597 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "P1|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.250000000000002, + 45.999999999999986 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 10.500000000000002, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P1|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.0, + 45.999999999999986 + ], + [ + 10.0, + 46.50000000000001 + ], + [ + 10.500000000000002, + 46.50000000000001 + ], + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.250000000000002, + 45.999999999999986 + ], + [ + 10.0, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 11.0, + 46.37499999999999 + ], + [ + 11.0, + 45.999999999999986 + ], + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.500000000000002, + 46.37499999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2|Z6" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.0, + 46.50000000000001 + ], + [ + 11.0, + 46.12500000000001 + ], + [ + 10.875, + 46.12500000000001 + ], + [ + 10.875, + 46.50000000000001 + ], + [ + 11.0, + 46.50000000000001 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 46.50000000000001 + ], + [ + 10.875, + 46.50000000000001 + ], + [ + 10.875, + 46.37499999999999 + ], + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 10.500000000000002, + 46.50000000000001 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 46.37499999999999 + ], + [ + 11.500000000000002, + 46.37499999999999 + ], + [ + 11.500000000000002, + 46.25 + ], + [ + 11.500000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3|Z6" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 46.625 + ], + [ + 11.375, + 46.625 + ], + [ + 11.375, + 46.12500000000001 + ], + [ + 11.250000000000002, + 46.12500000000001 + ], + [ + 11.250000000000002, + 46.625 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 45.999999999999986 + ], + [ + 11.500000000000002, + 45.999999999999986 + ], + [ + 11.500000000000002, + 46.25 + ], + [ + 11.749999999999998, + 46.25 + ], + [ + 11.749999999999998, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.37499999999999 + ], + [ + 11.375, + 46.37499999999999 + ], + [ + 11.375, + 46.625 + ], + [ + 11.250000000000002, + 46.625 + ], + [ + 11.250000000000002, + 46.74999999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P4|Z2" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.0, + 47.0 + ], + [ + 10.0, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.0 + ], + [ + 10.0, + 47.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P4|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.75, + 47.250000000000014 + ], + [ + 10.75, + 47.0 + ], + [ + 10.500000000000002, + 47.0 + ], + [ + 10.500000000000002, + 47.250000000000014 + ], + [ + 10.75, + 47.250000000000014 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P5|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 12.499999999999998, + 45.999999999999986 + ], + [ + 13.250000000000002, + 45.999999999999986 + ], + [ + 13.250000000000002, + 46.25 + ], + [ + 12.749999999999998, + 46.25 + ], + [ + 12.749999999999998, + 46.74999999999999 + ], + [ + 12.499999999999998, + 46.74999999999999 + ], + [ + 12.499999999999998, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P6|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.75, + 46.25 + ], + [ + 10.250000000000002, + 46.25 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.75, + 46.37499999999999 + ], + [ + 10.75, + 46.25 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P6|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.250000000000002, + 46.74999999999999 + ], + [ + 10.75, + 46.74999999999999 + ], + [ + 10.75, + 46.37499999999999 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.250000000000002, + 46.74999999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P7|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 12.250000000000002, + 46.74999999999999 + ], + [ + 12.250000000000002, + 47.0 + ], + [ + 12.000000000000002, + 47.0 + ], + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.250000000000002, + 46.74999999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P8|Z5" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 47.5 + ], + [ + 11.500000000000002, + 47.2505900350133 + ], + [ + 11.500000000000002, + 47.12499999999999 + ], + [ + 11.125, + 47.12499999999999 + ], + [ + 11.125, + 47.2505900350133 + ], + [ + 11.250000000000002, + 47.5 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P8|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.749999999999998, + 47.0 + ], + [ + 11.0, + 47.0 + ], + [ + 11.125, + 47.2505900350133 + ], + [ + 11.125, + 47.12499999999999 + ], + [ + 11.500000000000002, + 47.12499999999999 + ], + [ + 11.500000000000002, + 47.2505900350133 + ], + [ + 11.749999999999998, + 47.0 + ] + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/expected/identity_roads.geojson b/crates/analysis/tests/reference/expected/identity_roads.geojson new file mode 100644 index 000000000..276df012e --- /dev/null +++ b/crates/analysis/tests/reference/expected/identity_roads.geojson @@ -0,0 +1,245 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "r1|Z1" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.250000000000002, + 46.25 + ], + [ + 11.125, + 46.25 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r1|Z6" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.875, + 46.25 + ], + [ + 11.125, + 46.25 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r1|" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 9.875, + 46.25 + ], + [ + 10.250000000000002, + 46.25 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r2|Z1" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.500000000000002, + 46.0625 + ], + [ + 10.500000000000002, + 46.37499999999999 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r2|" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 10.500000000000002, + 46.4375 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r3|Z4" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 11.5625, + 46.31250000000001 + ], + [ + 11.6875, + 46.4375 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r4|" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 12.000000000000002, + 46.50000000000001 + ], + [ + 12.499999999999998, + 46.50000000000001 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r5|" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 14.999999999999998, + 49.00000000000001 + ], + [ + 15.5, + 49.5 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r6|Z5" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 11.206763423774959, + 47.12499999999999 + ], + [ + 11.500000000000002, + 47.34499451815567 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r6|" + }, + "geometry": { + "type": "MultiLineString", + "coordinates": [ + [ + [ + 10.875, + 46.87500000000001 + ], + [ + 11.206763423774959, + 47.12499999999999 + ] + ], + [ + [ + 11.500000000000002, + 47.34499451815567 + ], + [ + 11.875, + 47.625000000000014 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r7|Z1" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 11.125, + 46.37499999999999 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r7|Z6" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.875, + 46.37499999999999 + ], + [ + 11.125, + 46.37499999999999 + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/expected/union_parcels.geojson b/crates/analysis/tests/reference/expected/union_parcels.geojson new file mode 100644 index 000000000..e9586ebf4 --- /dev/null +++ b/crates/analysis/tests/reference/expected/union_parcels.geojson @@ -0,0 +1,857 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "P1|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.250000000000002, + 45.999999999999986 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 10.500000000000002, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P1|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.0, + 45.999999999999986 + ], + [ + 10.0, + 46.50000000000001 + ], + [ + 10.500000000000002, + 46.50000000000001 + ], + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.250000000000002, + 45.999999999999986 + ], + [ + 10.0, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 11.0, + 46.37499999999999 + ], + [ + 11.0, + 45.999999999999986 + ], + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.500000000000002, + 46.37499999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2|Z6" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.0, + 46.50000000000001 + ], + [ + 11.0, + 46.12500000000001 + ], + [ + 10.875, + 46.12500000000001 + ], + [ + 10.875, + 46.50000000000001 + ], + [ + 11.0, + 46.50000000000001 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 46.50000000000001 + ], + [ + 10.875, + 46.50000000000001 + ], + [ + 10.875, + 46.37499999999999 + ], + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 10.500000000000002, + 46.50000000000001 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 46.37499999999999 + ], + [ + 11.500000000000002, + 46.37499999999999 + ], + [ + 11.500000000000002, + 46.25 + ], + [ + 11.500000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3|Z6" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 46.625 + ], + [ + 11.375, + 46.625 + ], + [ + 11.375, + 46.12500000000001 + ], + [ + 11.250000000000002, + 46.12500000000001 + ], + [ + 11.250000000000002, + 46.625 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 45.999999999999986 + ], + [ + 11.500000000000002, + 45.999999999999986 + ], + [ + 11.500000000000002, + 46.25 + ], + [ + 11.749999999999998, + 46.25 + ], + [ + 11.749999999999998, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.37499999999999 + ], + [ + 11.375, + 46.37499999999999 + ], + [ + 11.375, + 46.625 + ], + [ + 11.250000000000002, + 46.625 + ], + [ + 11.250000000000002, + 46.74999999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P4|Z2" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.0, + 47.0 + ], + [ + 10.0, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.0 + ], + [ + 10.0, + 47.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P4|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.75, + 47.250000000000014 + ], + [ + 10.75, + 47.0 + ], + [ + 10.500000000000002, + 47.0 + ], + [ + 10.500000000000002, + 47.250000000000014 + ], + [ + 10.75, + 47.250000000000014 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P5|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 12.499999999999998, + 45.999999999999986 + ], + [ + 13.250000000000002, + 45.999999999999986 + ], + [ + 13.250000000000002, + 46.25 + ], + [ + 12.749999999999998, + 46.25 + ], + [ + 12.749999999999998, + 46.74999999999999 + ], + [ + 12.499999999999998, + 46.74999999999999 + ], + [ + 12.499999999999998, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P6|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.75, + 46.25 + ], + [ + 10.250000000000002, + 46.25 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.75, + 46.37499999999999 + ], + [ + 10.75, + 46.25 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P6|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.250000000000002, + 46.74999999999999 + ], + [ + 10.75, + 46.74999999999999 + ], + [ + 10.75, + 46.37499999999999 + ], + [ + 10.250000000000002, + 46.37499999999999 + ], + [ + 10.250000000000002, + 46.74999999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P7|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 12.250000000000002, + 46.74999999999999 + ], + [ + 12.250000000000002, + 47.0 + ], + [ + 12.000000000000002, + 47.0 + ], + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.250000000000002, + 46.74999999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P8|Z5" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 47.5 + ], + [ + 11.500000000000002, + 47.2505900350133 + ], + [ + 11.500000000000002, + 47.12499999999999 + ], + [ + 11.125, + 47.12499999999999 + ], + [ + 11.125, + 47.2505900350133 + ], + [ + 11.250000000000002, + 47.5 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P8|" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.749999999999998, + 47.0 + ], + [ + 11.0, + 47.0 + ], + [ + 11.125, + 47.2505900350133 + ], + [ + 11.125, + 47.12499999999999 + ], + [ + 11.500000000000002, + 47.12499999999999 + ], + [ + 11.500000000000002, + 47.2505900350133 + ], + [ + 11.749999999999998, + 47.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "|Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.500000000000002, + 45.75000000000001 + ], + [ + 10.250000000000002, + 45.75000000000001 + ], + [ + 10.250000000000002, + 45.999999999999986 + ], + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 11.0, + 45.999999999999986 + ], + [ + 11.0, + 46.37499999999999 + ], + [ + 11.250000000000002, + 46.37499999999999 + ], + [ + 11.250000000000002, + 45.999999999999986 + ], + [ + 11.500000000000002, + 45.999999999999986 + ], + [ + 11.500000000000002, + 45.75000000000001 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "|Z2" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 46.87500000000001 + ], + [ + 9.75, + 46.87500000000001 + ], + [ + 9.75, + 47.5 + ], + [ + 10.500000000000002, + 47.5 + ], + [ + 10.500000000000002, + 47.250000000000014 + ], + [ + 10.500000000000002, + 47.0 + ], + [ + 10.500000000000002, + 46.87500000000001 + ] + ], + [ + [ + 10.250000000000002, + 47.0 + ], + [ + 10.250000000000002, + 47.250000000000014 + ], + [ + 10.0, + 47.250000000000014 + ], + [ + 10.0, + 47.0 + ], + [ + 10.250000000000002, + 47.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "|Z3" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 14.0, + 45.999999999999986 + ], + [ + 14.0, + 46.50000000000001 + ], + [ + 13.499999999999998, + 46.50000000000001 + ], + [ + 13.499999999999998, + 45.999999999999986 + ], + [ + 14.0, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "|Z4" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.749999999999998, + 46.25 + ], + [ + 11.749999999999998, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.25 + ], + [ + 11.749999999999998, + 46.25 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "|Z5" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.125, + 47.75000000000001 + ], + [ + 11.500000000000002, + 47.75000000000001 + ], + [ + 11.500000000000002, + 47.2505900350133 + ], + [ + 11.250000000000002, + 47.5 + ], + [ + 11.125, + 47.2505900350133 + ], + [ + 11.125, + 47.75000000000001 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "|Z6" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.875, + 46.625 + ], + [ + 11.250000000000002, + 46.625 + ], + [ + 11.250000000000002, + 46.12500000000001 + ], + [ + 11.0, + 46.12500000000001 + ], + [ + 11.0, + 46.50000000000001 + ], + [ + 10.875, + 46.50000000000001 + ], + [ + 10.875, + 46.625 + ] + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/generate.py b/crates/analysis/tests/reference/generate.py index 0593b5174..935a365f4 100644 --- a/crates/analysis/tests/reference/generate.py +++ b/crates/analysis/tests/reference/generate.py @@ -235,13 +235,27 @@ def case(name, tool, params, expect): # Intersect: each input feature's part in each zone, keyed by both ids ("P1|Z1"), in input # order, then zone order. Parts that only touch (a shared edge or corner) are left out. The # zones' `id` is taken by the parcels', so it comes out as `id_2`. - def pieces(name, rows, part): - out = [] + # Identity adds each input feature's part outside the zones it has a part in, keyed "P1|" (no + # zone id); Union adds each zone's part outside the parcels that have a part in it, "|Z1". + def pieces(name, rows, part, rest=False, zone_rest=False): + out, covering = [], {zi: [] for zi, _ in ZONES} for i, g in rows: + hit = [] for zi, z in ZONES: r = part(m(g).intersection(m(z))) if r is not None and not r.is_empty: out.append((f"{i}|{zi}", ll(r))) + hit.append(m(z)) + covering[zi].append(m(g)) + if rest: + r = part(m(g).difference(unary_union(hit))) if hit else m(g) + if r is not None and not r.is_empty: + out.append((f"{i}|", ll(r))) + if zone_rest: + for zi, z in ZONES: + r = polygonal(m(z).difference(unary_union(covering[zi]))) if covering[zi] else m(z) + if r is not None and not r.is_empty: + out.append((f"|{zi}", ll(r))) return write_expected(name, out) both = ["id", "id_2"] @@ -251,6 +265,17 @@ def pieces(name, rows, part): {"features": pieces("intersect_roads", ROADS, linear), "key": both}) case("intersect points", "intersect", {"input": "points", "overlay": "zones"}, {"ids": [f"{i}|{zi}" for i, g in POINTS for zi, z in ZONES if m(z).intersects(m(g))], "key": both}) + case("identity polygons", "identity", {"input": "parcels", "overlay": "zones"}, + {"features": pieces("identity_parcels", parcels, polygonal, rest=True), "key": both}) + case("identity lines", "identity", {"input": "roads", "overlay": "zones"}, + {"features": pieces("identity_roads", ROADS, linear, rest=True), "key": both}) + ids = [] + for i, g in POINTS: + hits = [f"{i}|{zi}" for zi, z in ZONES if m(z).intersects(m(g))] + ids += hits or [f"{i}|"] + case("identity points", "identity", {"input": "points", "overlay": "zones"}, {"ids": ids, "key": both}) + case("union", "union", {"input": "parcels", "overlay": "zones"}, + {"features": pieces("union_parcels", parcels, polygonal, rest=True, zone_rest=True), "key": both}) # Dissolve: one feature per code, the union of its parcels. codes = sorted({c for _, c, _ in PARCELS}) diff --git a/crates/assistant/skills/geoprocessing/SKILL.md b/crates/assistant/skills/geoprocessing/SKILL.md index df374da38..604e466e3 100644 --- a/crates/assistant/skills/geoprocessing/SKILL.md +++ b/crates/assistant/skills/geoprocessing/SKILL.md @@ -12,7 +12,9 @@ Based on the RustGIS user guide, Analysis › Geoprocessing tools `describe_tool` for its parameters. Never guess parameter names. The built-in tools are: - proximity: `buffer` (`input`, `distance`, `units`, `dissolve`) - overlay: `clip` (`input`, `clip`), `erase` (`input`, `erase`), `intersect` (`input`, - `overlay`: the parts in each overlay polygon, with both layers' attributes) + `overlay`: the parts in each overlay polygon, with both layers' attributes), `identity` + (the same, plus the parts outside every overlay polygon), `union` (both polygon layers: + every piece either covers) - generalization: `dissolve` (`input`, `field`), `simplify` (`input`, `tolerance`) - geometry: `convexHull`, `featureToPoint`, `xyTableToPoint` (`input`, `x`, `y`, `crs`), `calculateGeometry` diff --git a/crates/engine/src/cmd/workflows.rs b/crates/engine/src/cmd/workflows.rs index 5f6800d40..967a6ec75 100644 --- a/crates/engine/src/cmd/workflows.rs +++ b/crates/engine/src/cmd/workflows.rs @@ -190,6 +190,8 @@ pub fn gallery_list(kind: &str) -> CmdResult { "clip", "erase", "intersect", + "identity", + "union", "dissolve", "convexHull", "featureToPoint", diff --git a/docs/comparison/index.md b/docs/comparison/index.md index 63b4bdafa..e651ce71c 100644 --- a/docs/comparison/index.md +++ b/docs/comparison/index.md @@ -142,8 +142,8 @@ Stated plainly, so the tables above are worth trusting: - **Coordinate systems.** Vector data is stored in longitude/latitude; there is no full PROJ. - **Editing depth.** Vertex move, insert and delete are supported, but not versioned or multi-user database editing. Edits are kept in the project and exported, not written back to databases. -- **Specialized analysis.** Built-in spatial joins, intersect and shortest paths on line networks - are available, but full network datasets, overlay union and identity tools and advanced spatial +- **Specialized analysis.** Built-in spatial joins, overlays (intersect, identity, union) and + shortest paths on line networks are available, but full network datasets and advanced spatial statistics depend on the external suite or another GIS. - **Python.** The [`rustgis` package](../python.md) runs commands headlessly, but there is no Python console inside the app. Python process plugins are supported on desktop and CLI, without a sandbox. diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 460df9b3b..00332e543 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -78,7 +78,7 @@ written some other way (a token typed into a SQL query, say) can still be there: before sharing it. A point buffered where the map was clicked is copied as the command that buffered it, with the point. -There are 21 built-in tools; GeoLibre and plugin tools join them when available. +There are 23 built-in tools; GeoLibre and plugin tools join them when available. Use these ids in `gp.run`'s `tool` parameter (`gp.tools` lists the available tools and their parameters): | Tool | Id | Toolbox | @@ -87,6 +87,8 @@ Use these ids in `gp.run`'s `tool` parameter (`gp.tools` lists the available too | Clip | `clip` | Analysis › Extract | | Erase | `erase` | Analysis › Overlay | | Intersect | `intersect` | Analysis › Overlay | +| Identity | `identity` | Analysis › Overlay | +| Union | `union` | Analysis › Overlay | | Spatial Join | `spatialJoin` | Analysis › Overlay | | Dissolve | `dissolve` | Data Management › Generalization | | Minimum Bounding Geometry | `convexHull` | Data Management › Features | @@ -170,6 +172,14 @@ the result: gp.run={"tool":"intersect","params":{"input":"parcels","overlay":"floodplain","output":"parcels_in_floodplain"}} ``` +**Identity** (`identity`) takes the same parameters. Each input feature gives its Intersect +features, then its part outside every overlay polygon, with empty overlay attributes. Every parcel +is then in the result, in or out of the floodplain. **Union** (`union`) takes two polygon +layers and adds, after Identity's features, each overlay polygon's part outside every input polygon, +with empty input attributes: together the features cover everything either layer covers. Both +follow Intersect's rules for field names, touching parts and selections. Where overlay polygons +overlap each other, an input feature's part in both comes out once for each. + **Spatial Join** (`spatialJoin`, also `gp.spatialJoin`) copies attributes from `join` features to `target` features. Layers may have different source coordinate systems: their geometry is stored in longitude/latitude before the join. `match` is one of: From 6708f9e5c771b108e21b45952310696d1bce066b Mon Sep 17 00:00:00 2001 From: mavrick Date: Mon, 26 Oct 2026 22:18:23 -0500 Subject: [PATCH 3/5] fix(analysis): Identity and Union review follow-ups - An overlay feature of several overlapping parts is unioned before the overlay, as Clip's mask is: Intersect, Identity and Union lost the overlap's area (x 1-2 of two squares 0-2 and 1-3 came out as two pieces). A single polygon is only oriented, so a point exactly on its edge stays on it. - Identity and Union keep input features without geometry, with empty overlay attributes, as Near will; Intersect still leaves them out. - Identity and Union get icons in the Tools pane, drawn in code. - The user guide counts 24 built-in tools. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/tools.rs | 55 +++++++++++++++++++++++++++++++++--- crates/ui-egui/src/icons.rs | 17 +++++++++++ crates/ui-egui/src/panes.rs | 2 ++ docs/user-guide/analysis.md | 3 +- 4 files changed, 72 insertions(+), 5 deletions(-) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index bc4879fa1..c22be0bfe 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -777,7 +777,7 @@ fn erase(inp: &dyn Inputs, p: &Json) -> Result { } /// An overlay polygon: its row in the overlay layer, its extent and its polygons in Web Mercator, -/// oriented for the overlay. +/// oriented for the overlay, its parts unioned when it has several (as Clip's mask is). struct Zone { row: usize, ext: Extent, @@ -795,7 +795,11 @@ fn overlay_zones<'a>(inp: &'a dyn Inputs, p: &Json) -> Result<(&'a Dataset, Vec< .into_iter() .filter_map(|row| { let g = ov.features.get(row)?.geometry.as_ref()?; - Some(Zone { row, ext: g.extent(), poly: to_multipolygon(g)?.orient(Direction::Default) }) + // Parts of one feature may overlap, which the overlay can't take as they are; one + // polygon is only oriented, so its edges stay exactly where they are. + let poly = to_multipolygon(g)?; + let poly = if poly.0.len() > 1 { union_all(poly.0.iter()) } else { poly.orient(Direction::Default) }; + Some(Zone { row, ext: g.extent(), poly }) }) .collect(); Ok((ov, zones)) @@ -849,7 +853,8 @@ impl Overlay { /// polygon, in overlay order, with the attributes of both; then (Identity, Union) its part outside /// them all, with empty overlay attributes. Union then adds each overlay polygon's part outside /// every input polygon, with empty input attributes. Parts that only touch are left out, and -/// features of one layer that overlap each other are not split against each other. +/// features of one layer that overlap each other are not split against each other. Intersect +/// leaves out input features without geometry; Identity and Union keep them, as they are. fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result { let ds = input(inp, p, "input")?; if ds.geometry_type == GeometryType::None { @@ -866,7 +871,13 @@ fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result>> = vec![Vec::new(); if op == Overlay::Union { zones.len() } else { 0 }]; for r in rows(inp, p, "input", ds) { let Some(f) = ds.features.get(r) else { continue }; - let Some(g) = &f.geometry else { continue }; + let Some(g) = &f.geometry else { + // Identity and Union keep every input feature: one without geometry, as it is. + if op != Overlay::Intersect { + out.push_next(None, overlay_row(ds, Some(&f.attrs), ov, None)); + } + continue; + }; let ext = g.extent(); let mut hit: Vec = Vec::new(); for zi in index.candidates(&ext) { @@ -3285,6 +3296,27 @@ mod tests { } } + #[test] + fn intersect_unions_an_overlay_feature_whose_parts_overlap() { + // One zone of two overlapping squares (x 0–2 and 1–3), as a file can hold it. + let ring = + |x: f64| vec![Coord::new(x, 0.0), Coord::new(x + 2.0, 0.0), Coord::new(x + 2.0, 2.0), Coord::new(x, 2.0), Coord::new(x, 0.0)]; + let mut zones = Dataset::new("zones", GeometryType::Polygon, vec![]); + zones.push(Some(Geometry::Polygon(vec![vec![ring(0.0)], vec![ring(1.0)]])), vec![]); + let mut parcels = Dataset::new("parcels", GeometryType::Polygon, vec![]); + parcels.push(Some(sq(-1.0, 0.0, 5.0)), vec![]); + let m = Map(HashMap::from([("zones".to_string(), zones), ("parcels".to_string(), parcels)])); + let out = layer(crate::run("intersect", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + assert_eq!(out.features.len(), 1); + // The part is the zone's whole outline, x 0–3, once: the overlap is neither lost nor doubled. + use geo::Area; + let part = to_multipolygon(out.features[0].geometry.as_ref().unwrap()).unwrap(); + let outline = vec![Coord::new(0.0, 0.0), Coord::new(3.0, 0.0), Coord::new(3.0, 2.0), Coord::new(0.0, 2.0), Coord::new(0.0, 0.0)]; + let want = to_multipolygon(&Geometry::Polygon(vec![vec![outline]])).unwrap(); + assert_eq!(part.0.len(), 1, "{part:?}"); + assert!((part.unsigned_area() - want.unsigned_area()).abs() < 1e-6 * want.unsigned_area(), "{part:?}"); + } + #[test] fn intersect_splits_points_by_each_polygon_they_fall_in() { let m = inputs(); @@ -3341,6 +3373,21 @@ mod tests { assert_eq!(out.features[2].geometry, m.0["parcels"].features[1].geometry); } + #[test] + fn identity_and_union_keep_input_features_without_geometry() { + let mut m = parcels_and_zones(); + m.0.get_mut("parcels").unwrap().push(None, vec![Value::Text("p3".into()), Value::Text("lot".into())]); + for (tool, keys) in [("identity", vec!["p1|z1", "p1|", "p2|", "p3|"]), ("union", vec!["p1|z1", "p1|", "p2|", "p3|", "|z1", "|z2"])] + { + let out = layer(crate::run(tool, &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + assert_eq!(overlay_keys(&out), keys, "{tool}"); + let p3 = out.features.iter().find(|f| f.attrs.first() == Some(&Value::Text("p3".into()))).unwrap(); + assert_eq!((p3.geometry.is_none(), &p3.attrs[1]), (true, &Value::Text("lot".into())), "{tool}"); + } + let out = layer(crate::run("intersect", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + assert_eq!(overlay_keys(&out), ["p1|z1"]); + } + #[test] fn identity_keeps_points_outside_every_polygon() { let m = inputs(); diff --git a/crates/ui-egui/src/icons.rs b/crates/ui-egui/src/icons.rs index 1f1bc3a9e..5a1c8d1aa 100644 --- a/crates/ui-egui/src/icons.rs +++ b/crates/ui-egui/src/icons.rs @@ -508,6 +508,21 @@ pub fn paint(painter: &Painter, r: Rect, name: &str, ink: Color32, accent: Color p.rect(3.0, 3.0, 12.0, 12.0, None, Some((1.3, ink)), 0.0); p.rect(9.0, 9.0, 12.0, 12.0, None, Some((1.3, ink)), 0.0); } + "identity" => { + // Intersect's squares, the input's square all kept: its part outside lighter. + p.rect(3.0, 3.0, 12.0, 12.0, Some(p.tint(60)), None, 0.0); + p.rect(9.0, 9.0, 6.0, 6.0, Some(accent), None, 0.0); + p.rect(3.0, 3.0, 12.0, 12.0, None, Some((1.3, ink)), 0.0); + p.rect(9.0, 9.0, 12.0, 12.0, None, Some((1.3, ink)), 0.0); + } + "union" => { + // Both squares kept, their shared part marked. + p.rect(3.0, 3.0, 12.0, 12.0, Some(p.tint(60)), None, 0.0); + p.rect(9.0, 9.0, 12.0, 12.0, Some(p.tint(60)), None, 0.0); + p.rect(9.0, 9.0, 6.0, 6.0, Some(accent), None, 0.0); + p.rect(3.0, 3.0, 12.0, 12.0, None, Some((1.3, ink)), 0.0); + p.rect(9.0, 9.0, 12.0, 12.0, None, Some((1.3, ink)), 0.0); + } "dissolve" => { p.polygon_shape(&[(2.0, 6.0), (12.0, 3.0), (22.0, 6.0), (20.0, 19.0), (4.0, 19.0)], 110); p.line(&[(12.0, 3.5), (12.0, 18.5)], 1.0, ink.gamma_multiply(0.35)); @@ -1209,6 +1224,8 @@ pub const ALL: &[&str] = &[ "clip", "erase", "intersect", + "identity", + "union", "dissolve", "hull", "centroid", diff --git a/crates/ui-egui/src/panes.rs b/crates/ui-egui/src/panes.rs index fae7d8c0c..83ffd3686 100644 --- a/crates/ui-egui/src/panes.rs +++ b/crates/ui-egui/src/panes.rs @@ -3394,6 +3394,8 @@ pub fn tool_icon(id: &str) -> &'static str { "clip" => "clip", "erase" => "erase", "intersect" => "intersect", + "identity" => "identity", + "union" => "union", "dissolve" => "dissolve", "convexHull" => "hull", "featureToPoint" => "centroid", diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 2e60a3486..eb537eef2 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -175,7 +175,8 @@ gp.run={"tool":"intersect","params":{"input":"parcels","overlay":"floodplain","o **Identity** (`identity`) takes the same parameters. Each input feature gives its Intersect features, then its part outside every overlay polygon, with empty overlay attributes. Every parcel -is then in the result, in or out of the floodplain. **Union** (`union`) takes two polygon +is then in the result, in or out of the floodplain; one without geometry is kept as it is (Intersect +leaves it out). **Union** (`union`) takes two polygon layers and adds, after Identity's features, each overlay polygon's part outside every input polygon, with empty input attributes: together the features cover everything either layer covers. Both follow Intersect's rules for field names, touching parts and selections. Where overlay polygons From 7622e7787746a22e06c65681cda1f3c598a2ebdf Mon Sep 17 00:00:00 2001 From: mavrick Date: Tue, 27 Oct 2026 00:29:06 -0500 Subject: [PATCH 4/5] fix(analysis): union an input feature's overlapping parts too An input feature of several overlapping polygons went into the overlay as it was, and geo's even-odd fill dropped the overlap: Identity's part outside a zone over one edge of two squares (x 0-2 and 1-3) lost x 1-2. part() now unions a feature of several polygons first, as the overlay zones are; a single polygon is passed as it is, in either winding. Union's zones keep the rows of the input polygons that cover them, not a copy of each polygon per zone. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/tools.rs | 58 ++++++++++++++++++++++++++---------- 1 file changed, 43 insertions(+), 15 deletions(-) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index c22be0bfe..19f887ce7 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -719,10 +719,15 @@ fn union_of(inp: &dyn Inputs, p: &Json, k: &str) -> Result, St /// them), or outside it when `inside` is false. Inside, a polygon's overlap with it, a line's /// stretches in it, the points in it or on its edge; outside, the rest. `None` when nothing is /// left, or only a shared edge or corner. +/// +/// The feature's polygons may wind either way (the reference corpus has a clockwise parcel); a +/// feature of several polygons is unioned first, as the overlay can't take parts that overlap. fn part(g: &Geometry, mask: &MultiPolygon, inside: bool) -> Option { match g { Geometry::Polygon(_) => { - to_multipolygon(g).and_then(|m| from_multipolygon(&if inside { m.intersection(mask) } else { m.difference(mask) })) + let m = to_multipolygon(g)?; + let m = if m.0.len() > 1 { union_all(m.0.iter()) } else { m }; + from_multipolygon(&if inside { m.intersection(mask) } else { m.difference(mask) }) } Geometry::Line(_) => to_multiline(g).and_then(|m| from_line(&mask.clip(&m, !inside))), Geometry::Point(_) | Geometry::MultiPoint(_) => { @@ -867,8 +872,8 @@ fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result>> = vec![Vec::new(); if op == Overlay::Union { zones.len() } else { 0 }]; + // For Union, the rows of the input polygons that overlap each zone. + let mut covering: Vec> = vec![Vec::new(); zones.len()]; for r in rows(inp, p, "input", ds) { let Some(f) = ds.features.get(r) else { continue }; let Some(g) = &f.geometry else { @@ -898,24 +903,26 @@ fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result> = + inputs.iter().filter_map(|r| ds.features.get(*r)?.geometry.as_ref()).filter_map(to_multipolygon).collect(); + from_multipolygon(&z.poly.difference(&union_all(polys.iter().flat_map(|m| m.0.iter())))) + }; + if let Some(rest) = rest { + out.push_next(Some(rest), overlay_row(ds, None, ov, zone_attrs(z))); + } } } let n = out.features.len(); @@ -3388,6 +3395,27 @@ mod tests { assert_eq!(overlay_keys(&out), ["p1|z1"]); } + #[test] + fn identity_unions_an_input_feature_whose_parts_overlap() { + // One parcel of two overlapping squares (x 0–2 and 1–3), and a zone over its right edge. + let ring = + |x: f64| vec![Coord::new(x, 0.0), Coord::new(x + 2.0, 0.0), Coord::new(x + 2.0, 2.0), Coord::new(x, 2.0), Coord::new(x, 0.0)]; + let mut parcels = Dataset::new("parcels", GeometryType::Polygon, vec![]); + parcels.push(Some(Geometry::Polygon(vec![vec![ring(0.0)], vec![ring(1.0)]])), vec![]); + let mut zones = Dataset::new("zones", GeometryType::Polygon, vec![]); + zones.push(Some(sq(2.5, 0.0, 2.0)), vec![]); + let m = Map(HashMap::from([("parcels".to_string(), parcels), ("zones".to_string(), zones)])); + let out = layer(crate::run("identity", &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + assert_eq!(out.features.len(), 2); + // The part outside is x 0–2.5 whole, the overlap (x 1–2) included. + use geo::Area; + let rest = to_multipolygon(out.features[1].geometry.as_ref().unwrap()).unwrap(); + let outline = vec![Coord::new(0.0, 0.0), Coord::new(2.5, 0.0), Coord::new(2.5, 2.0), Coord::new(0.0, 2.0), Coord::new(0.0, 0.0)]; + let want = to_multipolygon(&Geometry::Polygon(vec![vec![outline]])).unwrap(); + assert_eq!(rest.0.len(), 1, "{rest:?}"); + assert!((rest.unsigned_area() - want.unsigned_area()).abs() < 1e-6 * want.unsigned_area(), "{rest:?}"); + } + #[test] fn identity_keeps_points_outside_every_polygon() { let m = inputs(); From 28164b83ce2d2027babd45a4e90bbbb794194567 Mon Sep 17 00:00:00 2001 From: mavrick Date: Tue, 27 Oct 2026 02:34:08 -0500 Subject: [PATCH 5/5] fix(analysis): union an input feature once; cut both sides with one mask - A feature of several polygons is unioned once, by dissolved(), for all the zones it meets, not again in part() for each one. - Identity's part outside one zone is cut with that zone itself, the mask its inside part was cut with, not a union of it, which moves its edges by a rounding error. A road along a zone's edge comes out once, inside (new test). - part()'s doc says which polygons must be oriented (the mask's) and that the feature's must not overlap; the guide says Union's remainders of overlapping overlay polygons each keep their overlap. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/tools.rs | 64 +++++++++++++++++++++++++++++------- docs/user-guide/analysis.md | 4 ++- 2 files changed, 55 insertions(+), 13 deletions(-) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index 19f887ce7..df8def229 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -715,19 +715,28 @@ fn union_of(inp: &dyn Inputs, p: &Json, k: &str) -> Result, St Ok(union_all(polys.iter().flat_map(|m| m.0.iter()))) } -/// The part of a feature inside `mask` (polygons in Web Mercator, oriented as `union_all` leaves -/// them), or outside it when `inside` is false. Inside, a polygon's overlap with it, a line's -/// stretches in it, the points in it or on its edge; outside, the rest. `None` when nothing is -/// left, or only a shared edge or corner. +/// A feature as the overlay can take it: a feature of several polygons, which may overlap, unioned +/// into one; any other as it is, so a single polygon's edges stay exactly where they are. +fn dissolved(g: &Geometry) -> std::borrow::Cow<'_, Geometry> { + match g { + Geometry::Polygon(polys) if polys.len() > 1 => to_multipolygon(g) + .and_then(|m| from_multipolygon(&union_all(m.0.iter()))) + .map_or(std::borrow::Cow::Borrowed(g), std::borrow::Cow::Owned), + _ => std::borrow::Cow::Borrowed(g), + } +} + +/// The part of a feature inside `mask`, or outside it when `inside` is false. Inside, a polygon's +/// overlap with it, a line's stretches in it, the points in it or on its edge; outside, the rest. +/// `None` when nothing is left, or only a shared edge or corner. /// -/// The feature's polygons may wind either way (the reference corpus has a clockwise parcel); a -/// feature of several polygons is unioned first, as the overlay can't take parts that overlap. +/// The mask's polygons are in Web Mercator, oriented as `union_all` leaves them. The feature's may +/// wind either way (the reference corpus has a clockwise parcel), but must not overlap each other: +/// pass a feature of several polygons through [`dissolved`] first, once for all its masks. fn part(g: &Geometry, mask: &MultiPolygon, inside: bool) -> Option { match g { Geometry::Polygon(_) => { - let m = to_multipolygon(g)?; - let m = if m.0.len() > 1 { union_all(m.0.iter()) } else { m }; - from_multipolygon(&if inside { m.intersection(mask) } else { m.difference(mask) }) + to_multipolygon(g).and_then(|m| from_multipolygon(&if inside { m.intersection(mask) } else { m.difference(mask) })) } Geometry::Line(_) => to_multiline(g).and_then(|m| from_line(&mask.clip(&m, !inside))), Geometry::Point(_) | Geometry::MultiPoint(_) => { @@ -759,7 +768,7 @@ fn clip(inp: &dyn Inputs, p: &Json) -> Result { let mut out = copy_schema(ds, out_name(p, format!("{}_Clip", ds.name)), ds.geometry_type); for r in rows(inp, p, "input", ds) { let Some(f) = ds.features.get(r) else { continue }; - if let Some(c) = f.geometry.as_ref().and_then(|g| part(g, &mask, true)) { + if let Some(c) = f.geometry.as_ref().and_then(|g| part(&dissolved(g), &mask, true)) { out.push(Some(c), f.attrs.clone()); } } @@ -884,10 +893,12 @@ fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result = Vec::new(); for zi in index.candidates(&ext) { let Some(z) = zones.get(zi).filter(|z| ext.intersects(&z.ext)) else { continue }; - if let Some(inside) = part(g, &z.poly, true) { + if let Some(inside) = part(&whole, &z.poly, true) { out.push_next(Some(inside), overlay_row(ds, Some(&f.attrs), ov, zone_attrs(z))); hit.push(zi); } @@ -898,7 +909,12 @@ fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result zones.get(*only).and_then(|z| part(&whole, &z.poly, false)), + _ => part(&whole, &union_all(hit.iter().filter_map(|zi| zones.get(*zi)).flat_map(|z| z.poly.0.iter())), false), + } }; if let Some(rest) = rest { out.push_next(Some(rest), overlay_row(ds, Some(&f.attrs), ov, None)); @@ -3416,6 +3432,30 @@ mod tests { assert!((rest.unsigned_area() - want.unsigned_area()).abs() < 1e-6 * want.unsigned_area(), "{rest:?}"); } + #[test] + fn identity_keeps_a_line_along_a_zone_edge_once() { + // A road along the top edge of a zone, and on past its corner. + let mut roads = Dataset::new("roads", GeometryType::Polyline, vec![Field::new("id", FieldType::Text)]); + roads.push(Some(Geometry::Line(vec![vec![Coord::new(0.0, 2.0), Coord::new(3.0, 2.0)]])), vec![Value::Text("r".into())]); + let mut zones = Dataset::new("zones", GeometryType::Polygon, vec![Field::new("zone", FieldType::Text)]); + zones.push(Some(sq(0.0, 0.0, 2.0)), vec![Value::Text("z".into())]); + let m = Map(HashMap::from([("roads".to_string(), roads), ("zones".to_string(), zones)])); + let out = layer(crate::run("identity", &m, &serde_json::json!({"input": "roads", "overlay": "zones"})).unwrap()); + // The stretch on the edge is inside (x 0–2), once; the rest (x 2–3) outside. + let spans: Vec<(f64, f64, &Value)> = out + .features + .iter() + .map(|f| { + let e = f.geometry.as_ref().unwrap().extent(); + (e.min_x, e.max_x, &f.attrs[1]) + }) + .collect(); + assert_eq!(spans.len(), 2, "{spans:?}"); + let near = |a: f64, b: f64| (a - b).abs() < 1e-9; + assert!(near(spans[0].0, 0.0) && near(spans[0].1, 2.0) && *spans[0].2 == Value::Text("z".into()), "{spans:?}"); + assert!(near(spans[1].0, 2.0) && near(spans[1].1, 3.0) && *spans[1].2 == Value::Null, "{spans:?}"); + } + #[test] fn identity_keeps_points_outside_every_polygon() { let m = inputs(); diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index eb537eef2..7c22c3486 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -180,7 +180,9 @@ leaves it out). **Union** (`union`) takes two polygon layers and adds, after Identity's features, each overlay polygon's part outside every input polygon, with empty input attributes: together the features cover everything either layer covers. Both follow Intersect's rules for field names, touching parts and selections. Where overlay polygons -overlap each other, an input feature's part in both comes out once for each. +overlap each other, an input feature's part in both comes out once for each, and Union adds each +overlay polygon's own remainder, so the overlay's overlap outside the input comes out once for each +too. **Spatial Join** (`spatialJoin`, also `gp.spatialJoin`) copies attributes from `join` features to `target` features. Layers may have different source coordinate systems: their geometry is