From 85bb76cb9acb34a1e94537b1d827cba75af2c434 Mon Sep 17 00:00:00 2001 From: mavrick Date: Mon, 26 Oct 2026 13:28:41 -0500 Subject: [PATCH 1/8] 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/8] 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 ccf224e2696290ba931f94b66608d3876309eedb Mon Sep 17 00:00:00 2001 From: mavrick Date: Mon, 26 Oct 2026 16:59:16 -0500 Subject: [PATCH 3/8] feat(analysis): add a Near tool MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Near copies the input and adds, for each feature, the nearest feature of another layer: its OID, the ground distance to it in the chosen units, and the point on it nearest the input feature (a point they share when they intersect), within an optional search distance. With one layer as both, a feature's nearest is another feature. Spatial Join's Closest search moves, unchanged, into Nearby, which both tools use, and join::nearest_points returns the pair distance_m measures. Longitudes of the point are brought back within ±180°. Closes #959. Co-Authored-By: Claude Opus 5.5 --- GEOLIBRE_PARITY.md | 4 +- README.md | 2 +- ROADMAP.md | 11 +- crates/analysis/src/join.rs | 19 +- crates/analysis/src/tools.rs | 481 +++++++++++++++--- .../assistant/skills/geoprocessing/SKILL.md | 3 +- crates/engine/src/cmd/workflows.rs | 1 + docs/user-guide/analysis.md | 19 +- 8 files changed, 456 insertions(+), 84 deletions(-) diff --git a/GEOLIBRE_PARITY.md b/GEOLIBRE_PARITY.md index 43ab92678..33fc921f7 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, intersect, identity - and union are built in; near is missing (_Roadmap: solid ground_). +- [x] **Spatial join, intersect, union** ([#959](https://github.com/opengeos/RustGIS/issues/959)): + spatial join, clip, erase, intersect, identity, union and near are built in. - [ ] **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 1e4b02da3..edc3d89bc 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, 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 | +| **Analysis** | Selection by attributes (SQL) and location; Buffer, Near, 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 90875e9c2..5078c4aa9 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -60,7 +60,6 @@ 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). -- 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 @@ -219,11 +218,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, 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. + (clip, erase, intersect, identity, union), proximity (buffer, near) 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/join.rs b/crates/analysis/src/join.rs index 14bc2eac6..19f88e08f 100644 --- a/crates/analysis/src/join.rs +++ b/crates/analysis/src/join.rs @@ -1,6 +1,7 @@ -//! Spatial Join's measures and merge rules: the ground distance between two features, a test that -//! rules out features too far away to matter, the overlap a feature has with a polygon, and the -//! rules that roll the values of several matches up into one. +//! Spatial Join's and Near's measures, and Spatial Join's merge rules: the ground distance between +//! two features and their nearest points, a test that rules out features too far away to matter, +//! the overlap a feature has with a polygon, and the rules that roll the values of several matches +//! up into one. //! //! Distances are on the WGS 84 ellipsoid, as buffers' are. The nearest points of two features are //! found in an azimuthal equidistant plane, centred first between the features and then between @@ -35,12 +36,18 @@ pub fn first_coord(g: &Geometry) -> Option { /// in one centred between the pair found, where the plane is true near the answer however large /// the features are. pub fn distance_m(a: &Geometry, b: &Geometry) -> Option { + let (p, q) = nearest_points(a, b)?; + let d = wgs84_distance_m(p, q); + d.is_finite().then_some(d) +} + +/// The nearest points (lon/lat) of two features that don't intersect, as [`distance_m`] finds +/// them: a point of `a`, then one of `b`. `None` for a geometry without a finite coordinate. +pub fn nearest_points(a: &Geometry, b: &Geometry) -> Option<(Coord, Coord)> { // The features' first vertices are a pair of their points: the search starts from them. let (fa, fb) = (first_coord(a)?, first_coord(b)?); let (p, q) = nearest_pair(a, b, &Aeqd::new(geodesic_interpolate(fa, fb, 0.5)), Some((fa, fb)))?; - let (p, q) = nearest_pair(a, b, &Aeqd::new(geodesic_interpolate(p, q, 0.5)), Some((p, q))).unwrap_or((p, q)); - let d = wgs84_distance_m(p, q); - d.is_finite().then_some(d) + Some(nearest_pair(a, b, &Aeqd::new(geodesic_interpolate(p, q, 0.5)), Some((p, q))).unwrap_or((p, q))) } /// The nearest points of two geometries (lon/lat) as found in `plane`: a point of `a`, one of `b`. diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index 05d646588..0534d0bc9 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -106,6 +106,27 @@ pub static TOOLS: &[ToolSpec] = &[ ], run: dissolve, }, + ToolSpec { + id: "near", + name: "Near", + toolbox: "Analysis › Proximity", + summary: "Adds the distance to the nearest feature of another layer, its id and the nearest point on it.", + params: &[ + P { name: "input", label: "Input Features", kind: K::Layer, required: true, default: "", help: "" }, + P { name: "near", label: "Near Features", kind: K::Layer, required: true, default: "", help: "" }, + P { + name: "distance", + label: "Search Distance", + kind: K::Number, + required: false, + default: "", + help: "Features with nothing this near get empty values. Leave empty to search without a limit.", + }, + P { name: "units", label: "Distance Units", kind: K::Choice(UNITS), required: false, default: "\"Kilometers\"", help: "" }, + P { name: "output", label: "Output Name", kind: K::Text, required: false, default: "", help: "" }, + ], + run: near, + }, ToolSpec { id: "convexHull", name: "Minimum Bounding Geometry", @@ -1755,6 +1776,250 @@ impl ExtentIndex { } } +/// A feature to join or measure to: its row in its layer, its extent, and its geometry in Web +/// Mercator (`geo`) and in lon/lat (`geom`). +struct Located<'a> { + row: usize, + ext: Extent, + geo: geo::Geometry, + geom: &'a Geometry, +} + +/// Features to measure ground distances to (Spatial Join's join features, Near's near features), +/// looked up through an index of their extents. Features are named by their place in `feats`. +struct Nearby<'a, 'b> { + feats: &'b [Located<'a>], + index: &'b ExtentIndex, +} + +impl Nearby<'_, '_> { + /// The features within `reach` metres of a target (`t`, with its extent, and its geometry in Web + /// Mercator), in order, with their distances. With `nearest`, only the nearest is wanted: a + /// feature that can't be as near as the nearest so far isn't measured, so some farther ones + /// may be missing. The feature in row `exclude` is left out (the target itself, when both are + /// one layer). `measured` keeps each feature's distance, measured once per target. + fn within( + &self, + t: &Geometry, + t_ext: &Extent, + t_geo: &geo::Geometry, + reach: f64, + nearest: bool, + exclude: Option, + measured: &mut std::collections::HashMap>, + ) -> Vec<(usize, f64)> { + let looked = match crate::join::search_box(t_ext, reach) { + Some(b) => self.index.candidates(&b), + None => (0..self.feats.len()).collect(), + }; + let mut near = Vec::new(); + let mut bound = reach; + for idx in looked { + let Some(j) = self.feats.get(idx).filter(|j| Some(j.row) != exclude) else { continue }; + if !crate::join::may_be_within(t_ext, &j.ext, bound) { + continue; + } + let d = *measured.entry(idx).or_insert_with(|| ground_distance(t, t_geo, j.geom, &j.geo)); + if let Some(d) = d.filter(|d| *d <= reach * (1.0 + 1e-9) + 1e-6) { + near.push((idx, d)); + if nearest { + bound = bound.min(d); + } + } + } + near + } + + /// Every feature within `reach` metres of a target, nearest first; equal distances keep the + /// features' order. + fn all_within(&self, t: &Geometry, t_ext: &Extent, t_geo: &geo::Geometry, reach: f64) -> Vec<(usize, f64)> { + let mut near = self.within(t, t_ext, t_geo, reach, false, None, &mut std::collections::HashMap::new()); + near.sort_by(|a, b| a.1.total_cmp(&b.1)); + near + } + + /// The feature nearest a target, within `limit` metres when there is one, and its distance; + /// of equally near features, the first. The feature in row `exclude` is left out. + /// + /// Searches widen from about a grid cell, doubling until one finds a feature within its + /// distance (no nearer one can lie outside it), up to the limit. Without an index, one search + /// of everything. (A grid of coincident features has cells of no size: the first search reaches + /// at least `MIN_REACH_M`.) + fn closest( + &self, + t: &Geometry, + t_ext: &Extent, + t_geo: &geo::Geometry, + limit: Option, + exclude: Option, + ) -> Option<(usize, f64)> { + let mut measured = std::collections::HashMap::new(); + let mut reach = match self.index.cell_size() { + Some(cell) => (cell.to_radians() * rustgis_geom::EARTH_RADIUS_M).max(MIN_REACH_M), + None => f64::INFINITY, + }; + let mut near = loop { + let r = limit.map_or(reach, |l| reach.min(l)); + // A search reaching a pole or across the antimeridian looks at every feature anyway: it + // measures them out to the limit. + if crate::join::search_box(t_ext, r).is_none() { + break self.within(t, t_ext, t_geo, limit.unwrap_or(f64::INFINITY), true, exclude, &mut measured); + } + let found = self.within(t, t_ext, t_geo, r, true, exclude, &mut measured); + if !found.is_empty() || limit.is_some_and(|l| r >= l) { + break found; + } + reach *= 2.0; + }; + near.sort_by(|a, b| a.1.total_cmp(&b.1)); + near.into_iter().next() + } +} + +/// A point (lon/lat) two intersecting features share: one of the points of either on the other; +/// else a point of their overlap, when one is polygons; else a point where their edges meet. The +/// geometries come with their Web Mercator forms, where the test and the overlay work. +fn shared_point(a: &Geometry, a_geo: &geo::Geometry, b: &Geometry, b_geo: &geo::Geometry) -> Option { + use geo::LinesIter; + let back = |c: geo::Coord| rustgis_geom::merc_to_lonlat(Coord::new(c.x, c.y)); + // Points answer with their own coordinates, not ones back from Web Mercator. + let point_on = |pts: &Geometry, other: &geo::Geometry| -> Option> { + let pts: &[Coord] = match pts { + Geometry::Point(c) => std::slice::from_ref(c), + Geometry::MultiPoint(v) => v, + _ => return None, + }; + Some(pts.iter().copied().find(|c| { + let m = rustgis_geom::lonlat_to_merc(*c); + other.intersects(&geo::Point::new(m.x, m.y)) + })) + }; + if let Some(found) = point_on(a, b_geo).or_else(|| point_on(b, a_geo)) { + return found; + } + let overlap = match (a_geo, b_geo) { + (_, geo::Geometry::MultiPolygon(m)) => part(a, &m.orient(Direction::Default), true), + (geo::Geometry::MultiPolygon(m), _) => part(b, &m.orient(Direction::Default), true), + _ => None, + }; + if let Some(c) = overlap.as_ref().and_then(crate::join::first_coord) { + return Some(c); + } + // Lines, or polygons that only touch: edges by their west end, so each of `a`'s is tried only + // against those of `b` whose boxes may meet it. + let edges = |g: &geo::Geometry| -> Vec> { + let mut v: Vec> = match g { + geo::Geometry::MultiLineString(m) => m.lines_iter().collect(), + geo::Geometry::MultiPolygon(m) => m.lines_iter().collect(), + _ => Vec::new(), + }; + v.sort_by(|x, y| x.start.x.min(x.end.x).total_cmp(&y.start.x.min(y.end.x))); + v + }; + let (ea, eb) = (edges(a_geo), edges(b_geo)); + for la in &ea { + let (x0, x1) = (la.start.x.min(la.end.x), la.start.x.max(la.end.x)); + let (y0, y1) = (la.start.y.min(la.end.y), la.start.y.max(la.end.y)); + for lb in &eb { + if lb.start.x.min(lb.end.x) > x1 { + break; + } + if lb.start.x.max(lb.end.x) < x0 || lb.start.y.max(lb.end.y) < y0 || lb.start.y.min(lb.end.y) > y1 { + continue; + } + match geo::line_intersection::line_intersection(*la, *lb) { + Some(geo::LineIntersection::SinglePoint { intersection, .. }) => return Some(back(intersection)), + Some(geo::LineIntersection::Collinear { intersection }) => return Some(back(intersection.start)), + None => {} + } + } + } + None +} + +/// Near: a copy of the input features with the nearest near feature's OID (`near_fid`), the ground +/// distance to it (`near_dist`, in `units`) and the point on it nearest the input feature +/// (`near_x`, `near_y`, lon/lat; a point they share when they intersect). Within a search distance +/// when one is given; features with nothing that near, or without geometry, get empty values. +/// With one layer as both, a feature's nearest is another feature. +fn near(inp: &dyn Inputs, p: &Json) -> Result { + let ds = input(inp, p, "input")?; + let near = input(inp, p, "near")?; + if ds.geometry_type == GeometryType::None || near.geometry_type == GeometryType::None { + return Err(format!("'{}' has no geometry", if ds.geometry_type == GeometryType::None { &ds.name } else { &near.name })); + } + let given = p.get("distance").filter(|v| !v.is_null() && v.as_str().is_none_or(|t| !t.trim().is_empty())); + let limit = match given { + None => None, + Some(_) => { + let d = num(p, "distance").ok_or("Search Distance must be a number")?; + if d < 0.0 { + return Err("Search Distance can't be negative".into()); + } + Some(to_metres(d, s(p, "units"))?) + } + }; + let unit_m = to_metres(1.0, s(p, "units"))?; + // The same layer as both: leave each feature out of its own search. + let same = std::ptr::eq(ds, near); + + // In layer order, whatever order the selection was made in: of equally near features, the + // first in the layer wins. + let mut near_rows = rows(inp, p, "near", near); + near_rows.sort_unstable(); + let feats: Vec> = near_rows + .into_iter() + .filter_map(|row| { + let g = near.features.get(row)?.geometry.as_ref()?; + Some(Located { row, ext: g.extent(), geo: to_geo(g), geom: g }) + }) + .collect(); + let index = ExtentIndex::new(feats.iter().map(|f| f.ext)); + let nearby = Nearby { feats: &feats, index: &index }; + + let mut out = copy_schema(ds, out_name(p, format!("{}_Near", ds.name)), ds.geometry_type); + out.crs = ds.crs; + for (name, ty) in + [("near_fid", FieldType::Integer), ("near_dist", FieldType::Double), ("near_x", FieldType::Double), ("near_y", FieldType::Double)] + { + let name = unique_field(&out, name); + out.fields.push(Field::new(name, ty)); + } + let (mut found, mut total) = (0usize, 0usize); + for r in rows(inp, p, "input", ds) { + let Some(f) = ds.features.get(r) else { continue }; + total += 1; + let hit = f.geometry.as_ref().and_then(|g| { + let (t_ext, t_geo) = (g.extent(), to_geo(g)); + let (idx, d) = nearby.closest(g, &t_ext, &t_geo, limit, same.then_some(r))?; + let j = feats.get(idx)?; + let at = if d > 0.0 { + // The plane's longitudes run on past ±180°: back within them. + crate::join::nearest_points(g, j.geom).map(|(_, q)| Coord::new((q.x + 180.0).rem_euclid(360.0) - 180.0, q.y)) + } else { + shared_point(g, &t_geo, j.geom, &j.geo) + }; + Some((near.features.get(j.row)?.id, d, at)) + }); + let mut attrs: Vec = (0..ds.fields.len()).map(|i| f.attrs.get(i).cloned().unwrap_or(Value::Null)).collect(); + match hit { + Some((fid, d, at)) => { + found += 1; + let xy = |v: Option| v.map_or(Value::Null, Value::Float); + attrs.extend([ + i64::try_from(fid).map_or(Value::Null, Value::Int), + Value::Float(d / unit_m), + xy(at.map(|c| c.x)), + xy(at.map(|c| c.y)), + ]); + } + None => attrs.extend([Value::Null, Value::Null, Value::Null, Value::Null]), + } + out.push_next(f.geometry.clone(), attrs); + } + result(out, format!("Near: {found} of {total} features have a near feature")) +} + /// 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")?; @@ -1805,19 +2070,13 @@ fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result { let target_rows = rows(inp, p, "target", target); let join_rows = rows(inp, p, "join", join); - struct JoinGeom<'a> { - row: usize, - ext: Extent, - geo: geo::Geometry, - geom: &'a Geometry, - } - let mut join_geoms: Vec> = Vec::with_capacity(join_rows.len()); + let mut join_geoms: Vec> = Vec::with_capacity(join_rows.len()); for &jr in &join_rows { if let Some(f) = join.features.get(jr) && let Some(g) = f.geometry.as_ref() { let geo = if mode == JoinMatch::LargestOverlap { oriented(g) } else { to_geo(g) }; - join_geoms.push(JoinGeom { row: jr, ext: g.extent(), geo, geom: g }); + join_geoms.push(Located { row: jr, ext: g.extent(), geo, geom: g }); } } @@ -1847,6 +2106,7 @@ fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result { 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| index.candidates(ext); + let nearby = Nearby { feats: &join_geoms, index: &index }; let mut total_matches = 0usize; let mut written_features = 0usize; @@ -1898,67 +2158,12 @@ fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result { matches.extend(best.map(|(row, _)| (row, None))); } JoinMatch::Closest | JoinMatch::WithinDistance => { - // The join features within `reach` metres, in join order, measured once each. - let mut measured: std::collections::HashMap> = std::collections::HashMap::new(); - let mut within = |reach: f64| -> Vec<(usize, f64)> { - let looked = match crate::join::search_box(&t_ext, reach) { - Some(b) => candidates(&b), - None => (0..join_geoms.len()).collect(), - }; - let mut near = Vec::new(); - // Closest wants only the nearest: a feature that can't be as near as the - // nearest so far isn't measured. - let mut bound = reach; - for idx in looked { - let Some(j) = at(idx) else { continue }; - if !crate::join::may_be_within(&t_ext, &j.ext, bound) { - continue; - } - let d = *measured.entry(idx).or_insert_with(|| ground_distance(tg, &t_geo, j.geom, &j.geo)); - if let Some(d) = d.filter(|d| *d <= reach * (1.0 + 1e-9) + 1e-6) { - near.push((j.row, d)); - if mode == JoinMatch::Closest { - bound = bound.min(d); - } - } - } - near - }; - let mut near = match (mode, limit) { - (JoinMatch::WithinDistance, Some(l)) => within(l), + let near = match (mode, limit) { + (JoinMatch::WithinDistance, Some(l)) => nearby.all_within(tg, &t_ext, &t_geo, l), (JoinMatch::WithinDistance, None) => Vec::new(), - // Closest: searches widening from about a grid cell, doubling until one - // finds a feature within its distance (no nearer one can lie outside it), - // up to the search distance. Without an index, one search of everything. - // (A grid of coincident features has cells of no size: the first search - // reaches at least `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 - }; - loop { - let r = limit.map_or(reach, |l| reach.min(l)); - // A search reaching a pole or across the antimeridian looks at every - // feature anyway: it measures them out to the search distance. - if crate::join::search_box(&t_ext, r).is_none() { - break within(limit.unwrap_or(f64::INFINITY)); - } - let found = within(r); - if !found.is_empty() || limit.is_some_and(|l| r >= l) { - break found; - } - reach *= 2.0; - } - } + _ => nearby.closest(tg, &t_ext, &t_geo, limit, None).into_iter().collect(), }; - // Nearest first; equal distances keep the join order. - near.sort_by(|a, b| a.1.total_cmp(&b.1)); - if mode == JoinMatch::Closest { - near.truncate(1); - } - matches.extend(near.into_iter().map(|(row, d)| (row, Some(d)))); + matches.extend(near.into_iter().filter_map(|(idx, d)| Some((at(idx)?.row, Some(d))))); } } } @@ -3237,6 +3442,8 @@ mod tests { serde_json::json!({"input": "polys", "overlay": "polys"}), serde_json::json!({"input": "pts", "overlay": "pts", "output": ""}), serde_json::json!({"input": 3, "overlay": ["polys"]}), + serde_json::json!({"input": "pts", "near": "pts", "distance": -0.0, "units": "Feet"}), + serde_json::json!({"input": "polys", "near": "pts", "distance": 1e308}), ]; for t in TOOLS { for p in &junk { @@ -3348,6 +3555,69 @@ mod tests { assert!(e.contains("not a polygon layer"), "{e}"); } + /// Near's four fields of row `r`: the near feature's id, the distance and the point. + fn near_fields(out: &Dataset, r: usize) -> (Option, Option, Option, Option) { + let n = out.fields.len(); + let v = |k: usize| out.features[r].attrs[n - 4 + k].clone(); + (v(0).as_f64().map(|f| f as i64), v(1).as_f64(), v(2).as_f64(), v(3).as_f64()) + } + + #[test] + fn near_finds_the_nearest_feature_its_distance_and_point() { + let m = inputs(); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "pts", "near": "polys", "units": "Meters"})).unwrap()); + assert_eq!(out.name, "pts_Near"); + let names: Vec<&str> = out.fields.iter().map(|f| f.name.as_str()).collect(); + assert_eq!(names, ["v", "near_fid", "near_dist", "near_x", "near_y"]); + assert_eq!(out.features.len(), 3); + // Inside the first square: no distance, and the point itself is the nearest point. + assert_eq!(near_fields(&out, 0), (Some(1), Some(0.0), Some(0.5), Some(0.5))); + // (1.5, 1.5) is in both overlapping squares: the first wins the tie. + assert_eq!(near_fields(&out, 1).0, Some(1)); + // (50, 50) is nearest the third square's far corner, (11, 11), on the ground. + let (fid, d, x, y) = near_fields(&out, 2); + assert_eq!(fid, Some(3)); + let want = rustgis_geom::wgs84_distance_m(Coord::new(50.0, 50.0), Coord::new(11.0, 11.0)); + assert!((d.unwrap() - want).abs() < 1e-4 * want, "{d:?} vs {want}"); + assert!((x.unwrap() - 11.0).abs() < 1e-7 && (y.unwrap() - 11.0).abs() < 1e-7, "{x:?} {y:?}"); + } + + #[test] + fn near_leaves_features_beyond_the_search_distance_empty() { + let m = inputs(); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "pts", "near": "polys", "distance": 100})).unwrap()); + assert_eq!(near_fields(&out, 0).0, Some(1)); + assert_eq!(near_fields(&out, 2), (None, None, None, None), "(50, 50) is thousands of km away"); + for bad in [serde_json::json!(-1), serde_json::json!("far")] { + assert!(crate::run("near", &m, &serde_json::json!({"input": "pts", "near": "polys", "distance": bad})).is_err()); + } + } + + #[test] + fn near_within_one_layer_skips_the_feature_itself() { + let m = inputs(); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "pts", "near": "pts"})).unwrap()); + // Each point's nearest is another point, not itself at no distance. + assert_eq!(near_fields(&out, 0).0, Some(2)); + assert_eq!(near_fields(&out, 1).0, Some(1)); + assert!(near_fields(&out, 0).1.is_some_and(|d| d > 0.0)); + } + + #[test] + fn near_gives_a_shared_point_for_crossing_lines() { + let line = |a: (f64, f64), b: (f64, f64)| Geometry::Line(vec![vec![Coord::new(a.0, a.1), Coord::new(b.0, b.1)]]); + let mut a = Dataset::new("a", GeometryType::Polyline, vec![]); + a.push(Some(line((0.0, 0.0), (2.0, 2.0))), vec![]); + let mut b = Dataset::new("b", GeometryType::Polyline, vec![]); + b.push(Some(line((0.0, 2.0), (2.0, 0.0))), vec![]); + let m = Map(HashMap::from([("a".to_string(), a), ("b".to_string(), b)])); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "a", "near": "b"})).unwrap()); + let (fid, d, x, y) = near_fields(&out, 0); + assert_eq!((fid, d), (Some(1), Some(0.0))); + // Where they cross: (1, 1), give or take the map's straight lines. + assert!((x.unwrap() - 1.0).abs() < 1e-3 && (y.unwrap() - 1.0).abs() < 1e-3, "{x:?} {y:?}"); + } + /// Layers with a selection in each, as the Tools pane passes them. struct Selected(Map, HashMap>); impl Inputs for Selected { @@ -3394,6 +3664,83 @@ mod tests { } } + /// A point layer named `name` with a feature at each lon/lat in `at`, ids 1, 2, …. + fn points(name: &str, at: &[(f64, f64)]) -> Dataset { + let mut ds = Dataset::new(name, GeometryType::Point, vec![]); + for (x, y) in at { + ds.push_next(Some(Geometry::Point(Coord::new(*x, *y))), vec![]); + } + ds + } + + #[test] + fn near_searches_the_grid_and_breaks_ties_in_layer_order() { + // Enough points for a grid: a row 1° apart on the equator, and one far along it. + let mut at: Vec<(f64, f64)> = (0..20).map(|i| (f64::from(i), 0.0)).collect(); + at.push((100.0, 0.0)); + let m = Map(HashMap::from([("row".to_string(), points("row", &at))])); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "row", "near": "row"})).unwrap()); + // The first point's nearest is the second, not itself. + assert_eq!(near_fields(&out, 0).0, Some(2)); + // The eleventh (10°) is as near the tenth as the twelfth: the tenth wins. + let (fid, d, _, _) = near_fields(&out, 10); + assert_eq!(fid, Some(10)); + let degree = rustgis_geom::wgs84_distance_m(Coord::new(0.0, 0.0), Coord::new(1.0, 0.0)) / 1000.0; + assert!((d.unwrap() - degree).abs() < 1e-6 * degree, "{d:?} km vs {degree} km"); + // The far point's search widens until it reaches the last of the row, 81° away. + assert_eq!(near_fields(&out, 20).0, Some(20)); + } + + #[test] + fn near_measures_across_the_antimeridian() { + let m = Map(HashMap::from([ + ("a".to_string(), points("a", &[(179.9, 0.0)])), + ("b".to_string(), points("b", &[(178.0, 0.0), (-179.9, 0.0)])), + ])); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "a", "near": "b", "units": "Meters"})).unwrap()); + let (fid, d, x, y) = near_fields(&out, 0); + assert_eq!(fid, Some(2), "0.2° away across ±180°, not 1.9° on this side"); + let want = rustgis_geom::wgs84_distance_m(Coord::new(179.9, 0.0), Coord::new(-179.9, 0.0)); + assert!((d.unwrap() - want).abs() < 1e-6 * want, "{d:?} vs {want}"); + assert!((x.unwrap() + 179.9).abs() < 1e-7 && y.unwrap().abs() < 1e-7, "{x:?} {y:?}"); + } + + #[test] + fn near_finds_the_nearest_point_along_a_line_within_the_search_distance() { + let line = |a: (f64, f64), b: (f64, f64)| Geometry::Line(vec![vec![Coord::new(a.0, a.1), Coord::new(b.0, b.1)]]); + let mut roads = Dataset::new("roads", GeometryType::Polyline, vec![]); + roads.push(Some(line((-1.0, 0.0), (1.0, 0.0))), vec![]); + roads.push(Some(line((0.0, 5.0), (1.0, 5.0))), vec![]); + let m = Map(HashMap::from([("pts".to_string(), points("pts", &[(0.0, 1.0)])), ("roads".to_string(), roads)])); + let run = |p: serde_json::Value| layer(crate::run("near", &m, &p).unwrap()); + // A degree north of the first road's middle: the point below it, about 110.6 km away. + let (fid, d, x, y) = near_fields(&run(serde_json::json!({"input": "pts", "near": "roads"})), 0); + assert_eq!(fid, Some(1)); + let want = rustgis_geom::wgs84_distance_m(Coord::new(0.0, 1.0), Coord::new(0.0, 0.0)) / 1000.0; + assert!((d.unwrap() - want).abs() < 1e-4 * want, "{d:?} vs {want}"); + // Found on the plane's straight edge, within 0.01 % of the distance of the geodesic's point. + assert!(x.unwrap().abs() < 1e-4 && y.unwrap().abs() < 1e-4, "{x:?} {y:?}"); + let far = run(serde_json::json!({"input": "pts", "near": "roads", "distance": 100})); + assert_eq!(near_fields(&far, 0), (None, None, None, None)); + // Only the second road selected: it is the nearest. + let m = Selected(m, HashMap::from([("roads".to_string(), vec![1])])); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "pts", "near": "roads"})).unwrap()); + assert_eq!(near_fields(&out, 0).0, Some(2)); + } + + #[test] + fn near_gives_a_point_in_the_overlap_of_polygons() { + let m = inputs(); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "polys", "near": "polys"})).unwrap()); + // The first two squares overlap (x 1–2): each is the other's nearest, at a point they share. + for r in [0, 1] { + let (fid, d, x, y) = near_fields(&out, r); + assert_eq!((fid, d), (Some(2 - r as i64), Some(0.0))); + let (x, y) = (x.unwrap(), y.unwrap()); + assert!((1.0 - 1e-9..=2.0 + 1e-9).contains(&x) && (-1e-9..=2.0 + 1e-9).contains(&y), "{x} {y}"); + } + } + #[test] fn select_by_location_invert_excludes_features_without_geometry() { let mut polys = Dataset::new("polys", GeometryType::Polygon, vec![]); diff --git a/crates/assistant/skills/geoprocessing/SKILL.md b/crates/assistant/skills/geoprocessing/SKILL.md index 604e466e3..f9c1d5426 100644 --- a/crates/assistant/skills/geoprocessing/SKILL.md +++ b/crates/assistant/skills/geoprocessing/SKILL.md @@ -10,7 +10,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`) + - proximity: `buffer` (`input`, `distance`, `units`, `dissolve`), `near` (`input`, `near`, + `distance`, `units`: adds the nearest feature's id, its distance and the nearest point) - overlay: `clip` (`input`, `clip`), `erase` (`input`, `erase`), `intersect` (`input`, `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: diff --git a/crates/engine/src/cmd/workflows.rs b/crates/engine/src/cmd/workflows.rs index 967a6ec75..0b5200f52 100644 --- a/crates/engine/src/cmd/workflows.rs +++ b/crates/engine/src/cmd/workflows.rs @@ -187,6 +187,7 @@ pub fn gallery_list(kind: &str) -> CmdResult { "ready" => &["buffer", "dissolve", "countPoints", "clip", "summaryStatistics", "selectByLocation"], "feature" => &[ "buffer", + "near", "clip", "erase", "intersect", diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 00332e543..c7b8d7194 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -78,12 +78,13 @@ 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 23 built-in tools; GeoLibre and plugin tools join them when available. +There are 24 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 | | --- | --- | --- | | Buffer | `buffer` | Analysis › Proximity | +| Near | `near` | Analysis › Proximity | | Clip | `clip` | Analysis › Extract | | Erase | `erase` | Analysis › Overlay | | Intersect | `intersect` | Analysis › Overlay | @@ -180,6 +181,22 @@ with empty input attributes: together the features cover everything either layer 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. +**Near** (`near`) copies `input` and adds, for each feature, the nearest feature of `near`: its +OID as the attribute table shows it (`near_fid`), the distance to it (`near_dist`, in `units`, +Kilometers by default) and the point on it nearest the input feature (`near_x` and `near_y`, in +longitude and latitude). Distances are measured as Spatial Join's **Closest** measures them: on the +ground, between the features' nearest points, and zero for features that intersect, whose point is +one they share. The point is found as the distance is, in a plane centred between the features, +so on a long edge it can sit off the geodesic by up to about 0.01 % of the distance (5 m at 50 km). +A `distance` limits how far it looks; features with nothing that near, or without +geometry, get empty values. Of equally near features, the first in the near layer wins. With the +same layer as both, each feature's nearest is another feature, not itself. A field name the input +already has gets a suffix (`near_dist_2`). Selections in both layers are honoured: + +```text +gp.run={"tool":"near","params":{"input":"addresses","near":"hydrants","units":"Feet"}} +``` + **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 86dae9b28cb80f808e3b694828720bcadeb7b663 Mon Sep 17 00:00:00 2001 From: mavrick Date: Tue, 27 Oct 2026 07:33:36 -0500 Subject: [PATCH 4/8] fix(analysis): Near review follow-ups, and #982's held ones - Near's near layer is near_ds, not a local named like the function. - Near gets an icon in the Tools pane, drawn in code; the guide counts 25 built-in tools. - The overlay's masks share mask_of(); Near's shared point for a polygon unions a multi-part feature and mask as the overlay does. - From #982's review: an input feature no overlay polygon meets is documented as copied unchanged, and a multi-part feature is unioned only when an overlay polygon may meet it; the Identity icon's comment reads better. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/tools.rs | 49 ++++++++++++++++++++++-------------- crates/ui-egui/src/icons.rs | 11 +++++++- crates/ui-egui/src/panes.rs | 1 + docs/user-guide/analysis.md | 6 ++--- 4 files changed, 44 insertions(+), 23 deletions(-) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index 9d08ebb58..9c903fee8 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -819,6 +819,13 @@ struct Zone { poly: MultiPolygon, } +/// A feature's polygons as a mask for the overlay: parts of one feature may overlap, which the +/// overlay can't take as they are, so several are unioned; one polygon is only oriented, so its +/// edges stay exactly where they are. +fn mask_of(poly: MultiPolygon) -> MultiPolygon { + if poly.0.len() > 1 { union_all(poly.0.iter()) } else { poly.orient(Direction::Default) } +} + /// 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 overlay_zones<'a>(inp: &'a dyn Inputs, p: &Json) -> Result<(&'a Dataset, Vec), String> { @@ -830,11 +837,7 @@ 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()?; - // 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 }) + Some(Zone { row, ext: g.extent(), poly: mask_of(to_multipolygon(g)?) }) }) .collect(); Ok((ov, zones)) @@ -889,7 +892,9 @@ impl Overlay { /// 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. Intersect -/// leaves out input features without geometry; Identity and Union keep them, as they are. +/// leaves out input features without geometry; Identity and Union keep them, as they are, and an +/// input feature no overlay polygon meets too, coordinates untouched (a feature that is cut comes +/// out rebuilt, its overlapping parts unioned). fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result { let ds = input(inp, p, "input")?; if ds.geometry_type == GeometryType::None { @@ -914,11 +919,17 @@ fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result = index + .candidates(&ext) + .into_iter() + .filter_map(|zi| Some((zi, zones.get(zi)?))) + .filter(|(_, z)| ext.intersects(&z.ext)) + .collect(); + // A feature of overlapping parts is unioned once, for all the zones it may meet, and only + // when it may meet one. + let whole = if near.is_empty() { std::borrow::Cow::Borrowed(g) } else { dissolved(g) }; 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 }; + for (zi, z) in near { 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); @@ -1936,8 +1947,8 @@ fn shared_point(a: &Geometry, a_geo: &geo::Geometry, b: &Geometry, b_geo: & return found; } let overlap = match (a_geo, b_geo) { - (_, geo::Geometry::MultiPolygon(m)) => part(a, &m.orient(Direction::Default), true), - (geo::Geometry::MultiPolygon(m), _) => part(b, &m.orient(Direction::Default), true), + (_, geo::Geometry::MultiPolygon(m)) => part(&dissolved(a), &mask_of(m.clone()), true), + (geo::Geometry::MultiPolygon(m), _) => part(&dissolved(b), &mask_of(m.clone()), true), _ => None, }; if let Some(c) = overlap.as_ref().and_then(crate::join::first_coord) { @@ -1982,9 +1993,9 @@ fn shared_point(a: &Geometry, a_geo: &geo::Geometry, b: &Geometry, b_geo: & /// With one layer as both, a feature's nearest is another feature. fn near(inp: &dyn Inputs, p: &Json) -> Result { let ds = input(inp, p, "input")?; - let near = input(inp, p, "near")?; - if ds.geometry_type == GeometryType::None || near.geometry_type == GeometryType::None { - return Err(format!("'{}' has no geometry", if ds.geometry_type == GeometryType::None { &ds.name } else { &near.name })); + let near_ds = input(inp, p, "near")?; + if ds.geometry_type == GeometryType::None || near_ds.geometry_type == GeometryType::None { + return Err(format!("'{}' has no geometry", if ds.geometry_type == GeometryType::None { &ds.name } else { &near_ds.name })); } let given = p.get("distance").filter(|v| !v.is_null() && v.as_str().is_none_or(|t| !t.trim().is_empty())); let limit = match given { @@ -1999,16 +2010,16 @@ fn near(inp: &dyn Inputs, p: &Json) -> Result { }; let unit_m = to_metres(1.0, s(p, "units"))?; // The same layer as both: leave each feature out of its own search. - let same = std::ptr::eq(ds, near); + let same = std::ptr::eq(ds, near_ds); // In layer order, whatever order the selection was made in: of equally near features, the // first in the layer wins. - let mut near_rows = rows(inp, p, "near", near); + let mut near_rows = rows(inp, p, "near", near_ds); near_rows.sort_unstable(); let feats: Vec> = near_rows .into_iter() .filter_map(|row| { - let g = near.features.get(row)?.geometry.as_ref()?; + let g = near_ds.features.get(row)?.geometry.as_ref()?; Some(Located { row, ext: g.extent(), geo: to_geo(g), geom: g }) }) .collect(); @@ -2037,7 +2048,7 @@ fn near(inp: &dyn Inputs, p: &Json) -> Result { } else { shared_point(g, &t_geo, j.geom, &j.geo) }; - Some((near.features.get(j.row)?.id, d, at)) + Some((near_ds.features.get(j.row)?.id, d, at)) }); let mut attrs: Vec = (0..ds.fields.len()).map(|i| f.attrs.get(i).cloned().unwrap_or(Value::Null)).collect(); match hit { diff --git a/crates/ui-egui/src/icons.rs b/crates/ui-egui/src/icons.rs index 5a1c8d1aa..83b42355c 100644 --- a/crates/ui-egui/src/icons.rs +++ b/crates/ui-egui/src/icons.rs @@ -509,7 +509,8 @@ pub fn paint(painter: &Painter, r: Rect, name: &str, ink: Color32, accent: Color 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. + // Intersect's two squares, with all of the input's square kept: the overlap marked, the + // rest of the input 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); @@ -523,6 +524,13 @@ 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); } + "near" => { + // A point, and the shortest line from it to a polygon. + p.polygon_shape(&[(12.0, 3.0), (21.0, 4.0), (20.0, 12.0), (13.0, 11.0)], 60); + p.line(&[(6.0, 18.0), (13.0, 11.0)], 1.2, ink); + p.circle(13.0, 11.0, 1.6, Some(ink), None); + p.circle(6.0, 18.0, 2.8, Some(red), Some((1.0, Color32::WHITE))); + } "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)); @@ -1226,6 +1234,7 @@ pub const ALL: &[&str] = &[ "intersect", "identity", "union", + "near", "dissolve", "hull", "centroid", diff --git a/crates/ui-egui/src/panes.rs b/crates/ui-egui/src/panes.rs index 690a48121..343713e14 100644 --- a/crates/ui-egui/src/panes.rs +++ b/crates/ui-egui/src/panes.rs @@ -3415,6 +3415,7 @@ pub fn tool_icon(id: &str) -> &'static str { "intersect" => "intersect", "identity" => "identity", "union" => "union", + "near" => "near", "dissolve" => "dissolve", "convexHull" => "hull", "featureToPoint" => "centroid", diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 11b4424db..7b2545c35 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 24 built-in tools; GeoLibre and plugin tools join them when available. +There are 25 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 | @@ -176,8 +176,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; one without geometry is kept as it is (Intersect -leaves it out). **Union** (`union`) takes two polygon +is then in the result, in or out of the floodplain: one no overlay polygon meets, or without +geometry, as it is (Intersect leaves those 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 5b783c1d59afd345f809b38f3b30ac9028547156 Mon Sep 17 00:00:00 2001 From: giswqs <5016453+giswqs@users.noreply.github.com> Date: Tue, 27 Oct 2026 13:31:34 +0000 Subject: [PATCH 5/8] fix(analysis): Near gives every point longitude in -180 to 180 A shared point of intersecting features came back as the input's own longitude, so a point at 200 gave near_x 200, while the nearest-point path gave -160 for the same meridian. Both paths now normalize to [-180, 180), and the user guide says so. Co-Authored-By: Claude Haiku 5.5 Applied by the claude-autofix workflow. --- crates/analysis/src/tools.rs | 18 ++++++++++++++++-- docs/user-guide/analysis.md | 2 +- 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index 9c903fee8..acf6e5852 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -2043,11 +2043,13 @@ fn near(inp: &dyn Inputs, p: &Json) -> Result { let (idx, d) = nearby.closest(g, &t_ext, &t_geo, limit, same.then_some(r))?; let j = feats.get(idx)?; let at = if d > 0.0 { - // The plane's longitudes run on past ±180°: back within them. - crate::join::nearest_points(g, j.geom).map(|(_, q)| Coord::new((q.x + 180.0).rem_euclid(360.0) - 180.0, q.y)) + crate::join::nearest_points(g, j.geom).map(|(_, q)| q) } else { shared_point(g, &t_geo, j.geom, &j.geo) }; + // The plane's longitudes run on past ±180°, and so can a shared point's: back within + // [−180°, 180°). + let at = at.map(|c| Coord::new((c.x + 180.0).rem_euclid(360.0) - 180.0, c.y)); Some((near_ds.features.get(j.row)?.id, d, at)) }); let mut attrs: Vec = (0..ds.fields.len()).map(|i| f.attrs.get(i).cloned().unwrap_or(Value::Null)).collect(); @@ -3848,6 +3850,18 @@ mod tests { assert!((x.unwrap() + 179.9).abs() < 1e-7 && y.unwrap().abs() < 1e-7, "{x:?} {y:?}"); } + #[test] + fn near_gives_a_shared_point_within_180_degrees_of_the_meridian() { + let mut polys = Dataset::new("polys", GeometryType::Polygon, vec![]); + polys.push(Some(sq(190.0, -1.0, 20.0)), vec![]); + let m = Map(HashMap::from([("a".to_string(), points("a", &[(200.0, 0.0)])), ("b".to_string(), polys)])); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "a", "near": "b"})).unwrap()); + let (fid, d, x, y) = near_fields(&out, 0); + assert_eq!((fid, d), (Some(1), Some(0.0))); + // 200° is the meridian of −160°: the shared point comes back in [−180°, 180°) as the nearest points do. + assert!((x.unwrap() + 160.0).abs() < 1e-7 && y.unwrap().abs() < 1e-7, "{x:?} {y:?}"); + } + #[test] fn near_finds_the_nearest_point_along_a_line_within_the_search_distance() { let line = |a: (f64, f64), b: (f64, f64)| Geometry::Line(vec![vec![Coord::new(a.0, a.1), Coord::new(b.0, b.1)]]); diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 7b2545c35..907f7ce70 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -188,7 +188,7 @@ too. **Near** (`near`) copies `input` and adds, for each feature, the nearest feature of `near`: its OID as the attribute table shows it (`near_fid`), the distance to it (`near_dist`, in `units`, Kilometers by default) and the point on it nearest the input feature (`near_x` and `near_y`, in -longitude and latitude). Distances are measured as Spatial Join's **Closest** measures them: on the +longitude from −180° up to 180° and latitude). Distances are measured as Spatial Join's **Closest** measures them: on the ground, between the features' nearest points, and zero for features that intersect, whose point is one they share. The point is found as the distance is, in a plane centred between the features, so on a long edge it can sit off the geodesic by up to about 0.01 % of the distance (5 m at 50 km). From 063de663f3129f9598911b54be16a936fb8b243b Mon Sep 17 00:00:00 2001 From: mavrick Date: Tue, 27 Oct 2026 10:06:59 -0500 Subject: [PATCH 6/8] fix(analysis): format Near's point wrap; say what the edge sweep prunes - Format the autofix commit's line, which the lint check refused. - The shared-point edge search's comment says what it prunes and that two long touching borders still cost about the product of their edges, once per output row. - The overlay's candidate zones are `meets`, not a local named like the Near tool; a test pins that Identity and Union copy an input feature no zone cuts unchanged, parts and coordinates. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/tools.rs | 42 +++++++++++++++++++++++++++--------- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index acf6e5852..b734829e6 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -919,7 +919,7 @@ fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result = index + let meets: Vec<(usize, &Zone)> = index .candidates(&ext) .into_iter() .filter_map(|zi| Some((zi, zones.get(zi)?))) @@ -927,9 +927,9 @@ fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result = Vec::new(); - for (zi, z) in near { + for (zi, z) in meets { 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); @@ -1954,8 +1954,11 @@ fn shared_point(a: &Geometry, a_geo: &geo::Geometry, b: &Geometry, b_geo: & if let Some(c) = overlap.as_ref().and_then(crate::join::first_coord) { return Some(c); } - // Lines, or polygons that only touch: edges by their west end, so each of `a`'s is tried only - // against those of `b` whose boxes may meet it. + // Lines, or polygons that only touch: the first point where an edge of `a` meets one of `b`. `b`'s + // edges are sorted by their west end, so the search for each of `a`'s stops at the first that + // starts east of it, and a pair whose boxes don't meet is skipped before the exact test; pairs + // within that window are still each tried, so two long touching borders cost about the product + // of their edges. This runs once per output row, for the nearest feature only. let edges = |g: &geo::Geometry| -> Vec> { let mut v: Vec> = match g { geo::Geometry::MultiLineString(m) => m.lines_iter().collect(), @@ -2042,11 +2045,7 @@ fn near(inp: &dyn Inputs, p: &Json) -> Result { let (t_ext, t_geo) = (g.extent(), to_geo(g)); let (idx, d) = nearby.closest(g, &t_ext, &t_geo, limit, same.then_some(r))?; let j = feats.get(idx)?; - let at = if d > 0.0 { - crate::join::nearest_points(g, j.geom).map(|(_, q)| q) - } else { - shared_point(g, &t_geo, j.geom, &j.geo) - }; + let at = if d > 0.0 { crate::join::nearest_points(g, j.geom).map(|(_, q)| q) } else { shared_point(g, &t_geo, j.geom, &j.geo) }; // The plane's longitudes run on past ±180°, and so can a shared point's: back within // [−180°, 180°). let at = at.map(|c| Coord::new((c.x + 180.0).rem_euclid(360.0) - 180.0, c.y)); @@ -3676,6 +3675,29 @@ mod tests { assert!(near(spans[1].0, 2.0) && near(spans[1].1, 3.0) && *spans[1].2 == Value::Null, "{spans:?}"); } + #[test] + fn identity_and_union_copy_an_uncut_input_feature_unchanged() { + // Two overlapping squares as one parcel, far from the zone; and a triangle whose extent + // meets the zone's but whose shape doesn't reach it. + let ring = |x: f64| { + vec![Coord::new(x, 10.0), Coord::new(x + 2.0, 10.0), Coord::new(x + 2.0, 12.0), Coord::new(x, 12.0), Coord::new(x, 10.0)] + }; + let overlapping = Geometry::Polygon(vec![vec![ring(0.0)], vec![ring(1.0)]]); + let triangle = + Geometry::Polygon(vec![vec![vec![Coord::new(0.0, 0.0), Coord::new(4.0, 0.0), Coord::new(0.0, 4.0), Coord::new(0.0, 0.0)]]]); + let mut parcels = Dataset::new("parcels", GeometryType::Polygon, vec![]); + parcels.push(Some(overlapping.clone()), vec![]); + parcels.push(Some(triangle.clone()), vec![]); + let mut zones = Dataset::new("zones", GeometryType::Polygon, vec![]); + zones.push(Some(sq(3.0, 3.0, 1.0)), vec![]); + let m = Map(HashMap::from([("parcels".to_string(), parcels), ("zones".to_string(), zones)])); + for tool in ["identity", "union"] { + let out = layer(crate::run(tool, &m, &serde_json::json!({"input": "parcels", "overlay": "zones"})).unwrap()); + let got: Vec> = out.features.iter().take(2).map(|f| f.geometry.as_ref()).collect(); + assert_eq!(got, [Some(&overlapping), Some(&triangle)], "{tool}: copied as they are, parts and coordinates"); + } + } + #[test] fn identity_keeps_points_outside_every_polygon() { let m = inputs(); From 5a61c82ea33de00f90cb9053f779a006a773b886 Mon Sep 17 00:00:00 2001 From: mavrick Date: Tue, 27 Oct 2026 10:12:49 -0500 Subject: [PATCH 7/8] docs(analysis): rewrap the Near paragraph Co-Authored-By: Claude Opus 5.5 --- docs/user-guide/analysis.md | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 907f7ce70..f01b2eb71 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -185,17 +185,17 @@ overlap each other, an input feature's part in both comes out once for each, and overlay polygon's own remainder, so the overlay's overlap outside the input comes out once for each too. -**Near** (`near`) copies `input` and adds, for each feature, the nearest feature of `near`: its -OID as the attribute table shows it (`near_fid`), the distance to it (`near_dist`, in `units`, +**Near** (`near`) copies `input` and adds, for each feature, the nearest feature of `near`: its OID +as the attribute table shows it (`near_fid`), the distance to it (`near_dist`, in `units`, Kilometers by default) and the point on it nearest the input feature (`near_x` and `near_y`, in -longitude from −180° up to 180° and latitude). Distances are measured as Spatial Join's **Closest** measures them: on the -ground, between the features' nearest points, and zero for features that intersect, whose point is -one they share. The point is found as the distance is, in a plane centred between the features, -so on a long edge it can sit off the geodesic by up to about 0.01 % of the distance (5 m at 50 km). -A `distance` limits how far it looks; features with nothing that near, or without -geometry, get empty values. Of equally near features, the first in the near layer wins. With the -same layer as both, each feature's nearest is another feature, not itself. A field name the input -already has gets a suffix (`near_dist_2`). Selections in both layers are honoured: +longitude from −180° up to 180° and latitude). Distances are measured as Spatial Join's **Closest** +measures them: on the ground, between the features' nearest points, and zero for features that +intersect, whose point is one they share. The point is found as the distance is, in a plane centred +between the features, so on a long edge it can sit off the geodesic by up to about 0.01 % of the +distance (5 m at 50 km). A `distance` limits how far it looks; features with nothing that near, or +without geometry, get empty values. Of equally near features, the first in the near layer wins. +With the same layer as both, each feature's nearest is another feature, not itself. A field name +the input already has gets a suffix (`near_dist_2`). Selections in both layers are honoured: ```text gp.run={"tool":"near","params":{"input":"addresses","near":"hydrants","units":"Feet"}} From fed9888743ae61e8f039a3c0b5a8bcfa2ef0fe0a Mon Sep 17 00:00:00 2001 From: mavrick Date: Tue, 27 Oct 2026 14:37:59 -0500 Subject: [PATCH 8/8] fix(analysis): Near gives a point where no shared one is found MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two features that touch at a place written 180° on one side and -180° on the other are 0 m apart, but their Web Mercator forms don't meet, so shared_point found nothing and Near wrote empty coordinates. Near now falls back to the nearest points' pair whenever no shared point is found. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/tools.rs | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index b734829e6..c9be28f7e 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -2045,7 +2045,10 @@ fn near(inp: &dyn Inputs, p: &Json) -> Result { let (t_ext, t_geo) = (g.extent(), to_geo(g)); let (idx, d) = nearby.closest(g, &t_ext, &t_geo, limit, same.then_some(r))?; let j = feats.get(idx)?; - let at = if d > 0.0 { crate::join::nearest_points(g, j.geom).map(|(_, q)| q) } else { shared_point(g, &t_geo, j.geom, &j.geo) }; + // Where they touch, a point they share; else, or when none is found in the map's plane (the + // same place at −180° and 180°), the nearest points' pair. + let shared = if d > 0.0 { None } else { shared_point(g, &t_geo, j.geom, &j.geo) }; + let at = shared.or_else(|| crate::join::nearest_points(g, j.geom).map(|(_, q)| q)); // The plane's longitudes run on past ±180°, and so can a shared point's: back within // [−180°, 180°). let at = at.map(|c| Coord::new((c.x + 180.0).rem_euclid(360.0) - 180.0, c.y)); @@ -3884,6 +3887,16 @@ mod tests { assert!((x.unwrap() + 160.0).abs() < 1e-7 && y.unwrap().abs() < 1e-7, "{x:?} {y:?}"); } + #[test] + fn near_gives_a_point_for_the_same_place_on_both_sides_of_the_antimeridian() { + let m = Map(HashMap::from([("a".to_string(), points("a", &[(180.0, 0.0)])), ("b".to_string(), points("b", &[(-180.0, 0.0)]))])); + let out = layer(crate::run("near", &m, &serde_json::json!({"input": "a", "near": "b", "units": "Meters"})).unwrap()); + let (fid, d, x, y) = near_fields(&out, 0); + assert_eq!(fid, Some(1)); + assert!(d.is_some_and(|d| d.abs() < 1e-6), "{d:?}"); + assert!(x.is_some_and(|x| (x.abs() - 180.0).abs() < 1e-7) && y.is_some_and(|y| y.abs() < 1e-7), "{x:?} {y:?}"); + } + #[test] fn near_finds_the_nearest_point_along_a_line_within_the_search_distance() { let line = |a: (f64, f64), b: (f64, f64)| Geometry::Line(vec![vec![Coord::new(a.0, a.1), Coord::new(b.0, b.1)]]);