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..037922ee8 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 overlay_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) = overlay_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,98 @@ 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(); + } + // Only the part of `ext` over the grid: a feature far larger than the indexed layer would + // otherwise reach cells up to the clamp, past the grid's own. + let t = &self.total; + let (min_x, min_y, max_x, max_y) = (ext.min_x.max(t.min_x), ext.min_y.max(t.min_y), ext.max_x.min(t.max_x), ext.max_y.min(t.max_y)); + let ((x0, y0), (x1, y1)) = (self.key(min_x, min_y), self.key(max_x, 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 +1746,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 +1835,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 +3136,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 +3149,105 @@ 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 extent_index_looks_only_at_its_own_cells() { + // 20 small extents in a 1° row: a grid of about 4 cells a side. + let index = ExtentIndex::new((0..20).map(|i| Extent::new(f64::from(i) * 0.05, 0.0, f64::from(i) * 0.05 + 0.01, 0.01))); + assert!(index.cell_size().is_some()); + // A query far larger than the layer finds every feature once, in order. + assert_eq!(index.candidates(&Extent::new(-170.0, -80.0, 170.0, 80.0)), (0..20).collect::>()); + // A small one finds the features of its cell, not the whole row. + let near_start = index.candidates(&Extent::new(0.0, 0.0, 0.02, 0.005)); + assert!(near_start.first() == Some(&0) && near_start.len() < 20, "{near_start:?}"); + assert!(index.candidates(&Extent::new(5.0, 5.0, 6.0, 6.0)).is_empty()); + } + + #[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..647864c2c 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:?}")); @@ -175,11 +190,17 @@ fn check(case: &Json, layers: &Layers) -> Vec { let want: Vec = ids.iter().filter_map(json_text).collect(); let got: Vec = match result.output { Output::Selection { layer, rows } => match layers.0.get(&layer) { - Some(ds) => rows.iter().filter_map(|r| text(ds, *r, "id")).collect(), + Some(ds) => { + let key = key_fields(expect); + rows.iter().map(|r| key_text(ds, *r, key.as_deref())).collect() + } 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..90031c5b6 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 22 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 | @@ -101,6 +102,7 @@ Use these ids in `gp.run`'s `tool` parameter (`gp.tools` lists the available too | Summarize Within | `countPoints` | Analysis › Summarize | | Summary Statistics | `summaryStatistics` | Analysis › Statistics | | Export Features | `exportFeatures` | Conversion | +| Zonal Statistics | `zonalStats` | Analysis › Raster | | IDW Interpolation | `idwInterpolation` | Analysis › Interpolation | | Ordinary Kriging | `kriging` | Analysis › Interpolation | @@ -152,6 +154,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 crossed by two flood zones 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: