From 0294428245a573da7df858b7717323063cf60afa Mon Sep 17 00:00:00 2001 From: mavrick Date: Sun, 25 Oct 2026 19:18:20 -0500 Subject: [PATCH 1/2] test(analysis): check overlay, selection and join tools against GEOS A reference test runs the analysis tools on small made-up layers and compares their results with what GEOS (the engine of QGIS and PostGIS) gives for them: 21 cases across Clip, Erase, Dissolve, Minimum Bounding Geometry, Select Layer By Location, Summarize Within and Spatial Join. The layers put points on shared edges, corners and hole boundaries, and polygons that only touch, where engines tend to disagree. generate.py writes the layers and the expected results; GEOS runs only then and is no dependency of RustGIS. Co-Authored-By: Claude Opus 5.5 --- ATTRIBUTION.md | 1 + crates/analysis/tests/reference.rs | 241 +++++++ crates/analysis/tests/reference/cases.json | 614 ++++++++++++++++++ .../reference/expected/clip_parcels.geojson | 227 +++++++ .../reference/expected/clip_roads.geojson | 81 +++ .../reference/expected/dissolve_all.geojson | 230 +++++++ .../reference/expected/dissolve_code.geojson | 259 ++++++++ .../reference/expected/erase_parcels.geojson | 325 +++++++++ .../tests/reference/expected/hull_all.geojson | 46 ++ .../reference/expected/hull_each.geojson | 269 ++++++++ crates/analysis/tests/reference/generate.py | 337 ++++++++++ .../tests/reference/layers/parcels.geojson | 329 ++++++++++ .../tests/reference/layers/points.geojson | 135 ++++ .../tests/reference/layers/roads.geojson | 119 ++++ .../tests/reference/layers/stations.geojson | 31 + .../tests/reference/layers/zones.geojson | 203 ++++++ 16 files changed, 3447 insertions(+) create mode 100644 crates/analysis/tests/reference.rs create mode 100644 crates/analysis/tests/reference/cases.json create mode 100644 crates/analysis/tests/reference/expected/clip_parcels.geojson create mode 100644 crates/analysis/tests/reference/expected/clip_roads.geojson create mode 100644 crates/analysis/tests/reference/expected/dissolve_all.geojson create mode 100644 crates/analysis/tests/reference/expected/dissolve_code.geojson create mode 100644 crates/analysis/tests/reference/expected/erase_parcels.geojson create mode 100644 crates/analysis/tests/reference/expected/hull_all.geojson create mode 100644 crates/analysis/tests/reference/expected/hull_each.geojson create mode 100644 crates/analysis/tests/reference/generate.py create mode 100644 crates/analysis/tests/reference/layers/parcels.geojson create mode 100644 crates/analysis/tests/reference/layers/points.geojson create mode 100644 crates/analysis/tests/reference/layers/roads.geojson create mode 100644 crates/analysis/tests/reference/layers/stations.geojson create mode 100644 crates/analysis/tests/reference/layers/zones.geojson diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 0f72db35e..d444337d3 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -13,6 +13,7 @@ which every package ships next to LICENSE. They are generated from `Cargo.lock` | Sample data (`countries`, `countries_chn`, `cities`, `rivers`, `lakes` GeoJSON), embedded at build time from [rustgis-assets](https://github.com/opengeos/rustgis-assets/tree/main/data) | Natural Earth 1:110m, attributes trimmed; `countries_chn` is Natural Earth's China point of view (1:10m, simplified) | Public domain | | QGIS screenshot (`qgis-world-atlas.webp` in [rustgis-assets](https://github.com/opengeos/rustgis-assets/tree/main/images)), on the [RustGIS and QGIS](docs/comparison/qgis.md) page | Taken by the RustGIS project: [QGIS](https://qgis.org) 4.2 showing the Natural Earth sample data above | QGIS is GPL-2.0-or-later; the data is public domain | | QGIS test corpus (`crates/map/testdata/qgis/generated/`: styles, projects and their sample data) | Dominik Strebel ([@pick2510](https://github.com/pick2510)) for the RustGIS project: written by [QGIS](https://qgis.org) 4.2.3 through its API, from `generate.py` there; the data is made up | Apache-2.0 | +| Analysis reference corpus (`crates/analysis/tests/reference/`: layers and expected results) | Written for the RustGIS project from `generate.py` there: the layers are made up; the expected results were computed by [GEOS](https://libgeos.org) 3.13.1 through [Shapely](https://shapely.readthedocs.io) 2.1.2, run only to write them (GEOS is LGPL-2.1, Shapely BSD-3-Clause; neither is included in or linked into RustGIS) | Apache-2.0 | Raster coordinate systems are looked up by EPSG code in the PROJ definitions of the [crs-definitions](https://crates.io/crates/crs-definitions) crate (CC0-1.0), compiled into the diff --git a/crates/analysis/tests/reference.rs b/crates/analysis/tests/reference.rs new file mode 100644 index 000000000..4ef2e5506 --- /dev/null +++ b/crates/analysis/tests/reference.rs @@ -0,0 +1,241 @@ +//! Runs the analysis tools on the layers in `reference/layers` and compares what they give with the +//! answers in `reference/cases.json`, worked out by GEOS (the geometry engine of QGIS and PostGIS) +//! with `reference/generate.py`. Every case runs and every difference is reported together. +//! +//! When a case fails, fix the tool with RustGIS's own code. Don't translate GEOS's source to make +//! it pass: GEOS is LGPL-2.1, and RustGIS is Apache-2.0 with no third-party code. + +use std::collections::{BTreeMap, HashMap}; +use std::path::{Path, PathBuf}; + +use geo::{Area, BooleanOps, Euclidean, HausdorffDistance, Length, MultiLineString, MultiPolygon}; +use rustgis_analysis::convert::{to_multiline, to_multipolygon}; +use rustgis_analysis::{Inputs, Output}; +use rustgis_data::{Dataset, Geometry, Value}; +use serde_json::Value as Json; + +/// Polygons may differ by this fraction of their area (0.1 ppm: 160 m² of a 40 km square). +const AREA_TOLERANCE: f64 = 1e-7; +/// Lines may differ in length, and sit, this far (Web Mercator metres) from the reference. +const POSITION_TOLERANCE_M: f64 = 1e-3; +/// Distances may differ by this fraction, plus a millimetre. The reference measures along great +/// circles; the join finds the nearest points in an azimuthal equidistant plane (see `join.rs`), +/// whose straight edges are within this of the great circles, as its own tests require. +const DISTANCE_TOLERANCE: f64 = 1e-4; + +/// The reference corpus: the layers, the expected results and `cases.json`. +fn dir() -> PathBuf { + Path::new(env!("CARGO_MANIFEST_DIR")).join("tests/reference") +} + +/// Reads a GeoJSON file as a layer named `name`, the name the cases' parameters refer to it by. +fn read_layer(path: &Path, name: &str) -> Dataset { + let text = std::fs::read_to_string(path).unwrap_or_else(|e| panic!("{}: {e}", path.display())); + rustgis_data::io::read_geojson(&text, name).unwrap_or_else(|e| panic!("{}: {e}", path.display())).0 +} + +/// The reference layers by name, none with a selection: every tool works on all their features. +struct Layers(HashMap); + +impl Inputs for Layers { + fn dataset(&self, layer: &str) -> Option<&Dataset> { + self.0.get(layer) + } + fn selection(&self, _: &str) -> &[usize] { + &[] + } +} + +/// A field's value as text, for matching features by key (`None` for a missing field or null). +fn text(ds: &Dataset, row: usize, field: &str) -> Option { + let fi = ds.field_index(field)?; + match ds.value(row, fi) { + Value::Null => None, + Value::Text(t) => Some(t.clone()), + v => Some(v.to_string()), + } +} + +/// A key from `cases.json` as text, so it compares with [`text`] whether written as a string or a number. +fn json_text(v: &Json) -> Option { + match v { + Json::Null => None, + Json::String(s) => Some(s.clone()), + v => Some(v.to_string()), + } +} + +/// 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> { + 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 +} + +/// The polygons among `gs`, together, in Web Mercator metres. +fn polygons(gs: &[Geometry]) -> MultiPolygon { + MultiPolygon::new(gs.iter().filter_map(to_multipolygon).flat_map(|m| m.0).collect()) +} + +/// The lines among `gs`, together, in Web Mercator metres. +fn lines(gs: &[Geometry]) -> MultiLineString { + MultiLineString::new(gs.iter().filter_map(to_multiline).flat_map(|m| m.0).collect()) +} + +/// How one geometry differs from the reference, if it does beyond the tolerances. +fn compare_geometry(got: &[Geometry], want: &[Geometry]) -> Option { + if want.iter().any(|g| matches!(g, Geometry::Polygon(_))) { + let (g, w) = (polygons(got), polygons(want)); + // The area the two don't share, not vertex positions: GEOS keeps vertices partway along a + // straight edge that other engines drop, with the same shape. + let (area, off) = (w.unsigned_area(), g.xor(&w).unsigned_area()); + (off > AREA_TOLERANCE * area) + .then(|| format!("area {:.1} m² vs {area:.1} m²; {off:.1} m² ({:.2e} of it) not shared", g.unsigned_area(), off / area)) + } else { + let (g, w) = (lines(got), lines(want)); + let (gl, wl) = (Euclidean.length(&g), Euclidean.length(&w)); + let far = if g.0.is_empty() { f64::INFINITY } else { g.hausdorff_distance(&w) }; + ((gl - wl).abs() > POSITION_TOLERANCE_M || far > POSITION_TOLERANCE_M) + .then(|| format!("length {gl:.3} m vs {wl:.3} m; up to {far:.3} m apart")) + } +} + +/// The dataset a tool's output carries, whether a new layer, an updated one or a table. +fn layer_of(output: Output) -> Option { + match output { + Output::Layer(ds) | Output::Update { dataset: ds, .. } | Output::Table(ds) => Some(ds), + _ => None, + } +} + +/// Checks one case; the differences found, empty when it matches. +fn check(case: &Json, layers: &Layers) -> Vec { + let tool = case["tool"].as_str().unwrap_or_default(); + let expect = &case["expect"]; + let result = match rustgis_analysis::run(tool, layers, &case["params"]) { + Ok(r) => r, + Err(e) => return vec![format!("the tool failed: {e}")], + }; + let mut problems = Vec::new(); + + 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); + for (k, w) in &want { + match got.get(k) { + None => problems.push(format!("{k}: missing (the reference has it)")), + Some(g) => problems.extend(compare_geometry(g, w).map(|d| format!("{k}: {d}"))), + } + } + for k in got.keys().filter(|k| !want.contains_key(*k)) { + problems.push(format!("{k}: present (the reference has nothing)")); + } + if let Some(counts) = expect["counts"].as_object() { + for (row, _) in ds.features.iter().enumerate() { + let k = key.and_then(|key| text(&ds, row, key)).unwrap_or_default(); + 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:?}")); + } + } + } + } else if let Some(ids) = expect["ids"].as_array() { + let want: Vec = ids.iter().filter_map(json_text).collect(); + let got: Vec = match result.output { + Output::Selection { layer, rows } => { + let ds = &layers.0[&layer]; + rows.iter().filter_map(|r| text(ds, *r, "id")).collect() + } + other => match layer_of(other) { + Some(ds) => (0..ds.features.len()).filter_map(|r| text(&ds, r, "id")).collect(), + None => return vec!["neither a selection nor a layer came out".into()], + }, + }; + if got != want { + problems.push(format!("gave {got:?}, the reference {want:?}")); + } + } else if let Some(values) = expect["values"].as_object() { + let Some(ds) = layer_of(result.output) else { return vec!["no layer came out".into()] }; + let field = expect["field"].as_str().unwrap_or_default(); + for row in 0..ds.features.len() { + let id = text(&ds, row, "id").unwrap_or_default(); + let (got, want) = (text(&ds, row, field), values.get(&id).and_then(json_text)); + if got != want { + problems.push(format!("{id}: {field} {got:?}, the reference {want:?}")); + } + } + } else if let Some(pairs) = expect["pairs"].as_array() { + let Some(ds) = layer_of(result.output) else { return vec!["no layer came out".into()] }; + // Per target, its joined features in order, each with its distance when one is written. + type Joined = BTreeMap, Option)>>; + let mut want = Joined::new(); + for p in pairs { + let d = p.get(2).and_then(Json::as_f64); + want.entry(json_text(&p[0]).unwrap_or_default()).or_default().push((json_text(&p[1]), d)); + } + let mut got = Joined::new(); + let dist = ds.field_index("dist"); + for row in 0..ds.features.len() { + let d = dist.and_then(|fi| ds.value(row, fi).as_f64()); + got.entry(text(&ds, row, "id").unwrap_or_default()).or_default().push((text(&ds, row, "join_id"), d)); + } + // Order matters only for Within a distance (nearest first). + let ordered = case["params"]["match"] == "Within a distance"; + for (t, w) in &want { + let mut g = got.get(t).cloned().unwrap_or_default(); + let mut w = w.clone(); + if !ordered { + g.sort_by(|a, b| a.0.cmp(&b.0)); + w.sort_by(|a, b| a.0.cmp(&b.0)); + } + let ids = |v: &[(Option, Option)]| v.iter().map(|p| p.0.clone()).collect::>(); + if ids(&g) != ids(&w) { + problems.push(format!("{t}: joined {:?}, the reference {:?}", ids(&g), ids(&w))); + continue; + } + for ((j, gd), (_, wd)) in g.iter().zip(&w) { + if let (Some(gd), Some(wd)) = (gd, wd) + && (gd - wd).abs() > DISTANCE_TOLERANCE * wd + 1e-3 + { + problems.push(format!("{t}→{}: {gd:.4} m, the reference {wd:.4} m", j.clone().unwrap_or_default())); + } else if gd.is_some() != wd.is_some() { + problems.push(format!("{t}→{}: distance {gd:?}, the reference {wd:?}", j.clone().unwrap_or_default())); + } + } + } + for t in got.keys().filter(|t| !want.contains_key(*t)) { + problems.push(format!("{t}: present (the reference has no such target)")); + } + } else { + problems.push("the case says nothing to expect".into()); + } + problems +} + +#[test] +fn tools_agree_with_the_reference() { + let cases: Json = serde_json::from_str(&std::fs::read_to_string(dir().join("cases.json")).unwrap()).unwrap(); + let layers = Layers( + ["parcels", "zones", "points", "roads", "stations"] + .into_iter() + .map(|n| (n.to_string(), read_layer(&dir().join(format!("layers/{n}.geojson")), n))) + .collect(), + ); + let mut report = Vec::new(); + for case in cases["cases"].as_array().unwrap() { + let problems = check(case, &layers); + if !problems.is_empty() { + report.push(format!("{}:\n {}", case["name"].as_str().unwrap_or("?"), problems.join("\n "))); + } + } + let n = cases["cases"].as_array().map_or(0, Vec::len); + assert!(report.is_empty(), "{} of {n} cases differ from the reference:\n{}", report.len(), report.join("\n")); +} diff --git a/crates/analysis/tests/reference/cases.json b/crates/analysis/tests/reference/cases.json new file mode 100644 index 000000000..9783c70da --- /dev/null +++ b/crates/analysis/tests/reference/cases.json @@ -0,0 +1,614 @@ +{ + "about": "Answers from GEOS; regenerate with `uv run generate.py` (see its header).", + "shapely": "2.1.2", + "geos": "3.13.1", + "cases": [ + { + "name": "clip polygons", + "tool": "clip", + "params": { + "input": "parcels", + "clip": "zones" + }, + "expect": { + "features": "expected/clip_parcels.geojson" + } + }, + { + "name": "clip lines", + "tool": "clip", + "params": { + "input": "roads", + "clip": "zones" + }, + "expect": { + "features": "expected/clip_roads.geojson" + } + }, + { + "name": "clip points", + "tool": "clip", + "params": { + "input": "points", + "clip": "zones" + }, + "expect": { + "ids": [ + "b", + "d", + "e", + "j" + ] + } + }, + { + "name": "erase", + "tool": "erase", + "params": { + "input": "parcels", + "erase": "zones" + }, + "expect": { + "features": "expected/erase_parcels.geojson" + } + }, + { + "name": "dissolve by field", + "tool": "dissolve", + "params": { + "input": "parcels", + "field": "code" + }, + "expect": { + "features": "expected/dissolve_code.geojson", + "key": "code", + "counts": { + "A": 3, + "B": 2, + "C": 2, + "D": 1 + } + } + }, + { + "name": "dissolve all", + "tool": "dissolve", + "params": { + "input": "parcels" + }, + "expect": { + "features": "expected/dissolve_all.geojson", + "key": null + } + }, + { + "name": "convex hull of each", + "tool": "convexHull", + "params": { + "input": "parcels", + "group": "None" + }, + "expect": { + "features": "expected/hull_each.geojson" + } + }, + { + "name": "convex hull of all", + "tool": "convexHull", + "params": { + "input": "parcels", + "group": "All" + }, + "expect": { + "features": "expected/hull_all.geojson", + "key": null + } + }, + { + "name": "select parcels by zones", + "tool": "selectByLocation", + "params": { + "input": "parcels", + "select": "zones" + }, + "expect": { + "ids": [ + "P1", + "P2", + "P3", + "P4", + "P6", + "P8" + ] + } + }, + { + "name": "select parcels by zones, inverted", + "tool": "selectByLocation", + "params": { + "input": "parcels", + "select": "zones", + "invert": "Yes" + }, + "expect": { + "ids": [ + "P5", + "P7" + ] + } + }, + { + "name": "select roads by parcels", + "tool": "selectByLocation", + "params": { + "input": "roads", + "select": "parcels" + }, + "expect": { + "ids": [ + "r1", + "r2", + "r4", + "r6" + ] + } + }, + { + "name": "select points by parcels", + "tool": "selectByLocation", + "params": { + "input": "points", + "select": "parcels" + }, + "expect": { + "ids": [ + "a", + "b", + "c", + "e", + "g", + "h", + "i", + "j" + ] + } + }, + { + "name": "select roads by zones", + "tool": "selectByLocation", + "params": { + "input": "roads", + "select": "zones" + }, + "expect": { + "ids": [ + "r1", + "r2", + "r3", + "r6" + ] + } + }, + { + "name": "summarize within", + "tool": "countPoints", + "params": { + "polygons": "parcels", + "points": "points" + }, + "expect": { + "field": "point_count", + "values": { + "P1": 3, + "P2": 2, + "P3": 2, + "P4": 0, + "P5": 2, + "P6": 1, + "P7": 1, + "P8": 1 + } + } + }, + { + "name": "join intersects", + "tool": "spatialJoin", + "params": { + "target": "parcels", + "join": "zones", + "match": "Intersects", + "operation": "Join one to many" + }, + "expect": { + "pairs": [ + [ + "P1", + "Z1" + ], + [ + "P2", + "Z1" + ], + [ + "P2", + "Z6" + ], + [ + "P3", + "Z1" + ], + [ + "P3", + "Z4" + ], + [ + "P3", + "Z6" + ], + [ + "P4", + "Z2" + ], + [ + "P5", + null + ], + [ + "P6", + "Z1" + ], + [ + "P7", + null + ], + [ + "P8", + "Z5" + ] + ] + } + }, + { + "name": "join contains", + "tool": "spatialJoin", + "params": { + "target": "parcels", + "join": "points", + "match": "Contains", + "operation": "Join one to many" + }, + "expect": { + "pairs": [ + [ + "P1", + "a" + ], + [ + "P2", + null + ], + [ + "P3", + null + ], + [ + "P4", + null + ], + [ + "P5", + "g" + ], + [ + "P5", + "h" + ], + [ + "P6", + "c" + ], + [ + "P7", + null + ], + [ + "P8", + "j" + ] + ] + } + }, + { + "name": "join within", + "tool": "spatialJoin", + "params": { + "target": "points", + "join": "parcels", + "match": "Within", + "operation": "Join one to many" + }, + "expect": { + "pairs": [ + [ + "a", + "P1" + ], + [ + "b", + null + ], + [ + "c", + "P6" + ], + [ + "d", + null + ], + [ + "e", + null + ], + [ + "f", + null + ], + [ + "g", + "P5" + ], + [ + "h", + "P5" + ], + [ + "i", + null + ], + [ + "j", + "P8" + ] + ] + } + }, + { + "name": "join roads intersect parcels", + "tool": "spatialJoin", + "params": { + "target": "roads", + "join": "parcels", + "match": "Intersects", + "operation": "Join one to many" + }, + "expect": { + "pairs": [ + [ + "r1", + "P1" + ], + [ + "r1", + "P2" + ], + [ + "r1", + "P6" + ], + [ + "r2", + "P1" + ], + [ + "r2", + "P2" + ], + [ + "r2", + "P6" + ], + [ + "r3", + null + ], + [ + "r4", + "P3" + ], + [ + "r4", + "P5" + ], + [ + "r5", + null + ], + [ + "r6", + "P8" + ] + ] + } + }, + { + "name": "join closest", + "tool": "spatialJoin", + "params": { + "target": "points", + "join": "parcels", + "match": "Closest", + "units": "Meters", + "distanceField": "dist", + "operation": "Join one to many" + }, + "expect": { + "pairs": [ + [ + "a", + "P1", + 0.0 + ], + [ + "b", + "P1", + 0.0 + ], + [ + "c", + "P1", + 0.0 + ], + [ + "d", + "P3", + 9589.67458520921 + ], + [ + "e", + "P3", + 0.0 + ], + [ + "f", + "P5", + 191521.6332492596 + ], + [ + "g", + "P5", + 0.0 + ], + [ + "h", + "P5", + 0.0 + ], + [ + "i", + "P3", + 0.0 + ], + [ + "j", + "P8", + 0.0 + ] + ] + } + }, + { + "name": "join within a distance", + "tool": "spatialJoin", + "params": { + "target": "stations", + "join": "roads", + "match": "Within a distance", + "distance": 30000.0, + "units": "Meters", + "distanceField": "dist", + "operation": "Join one to many" + }, + "expect": { + "pairs": [ + [ + "s1", + "r2", + 20512.8663375754 + ], + [ + "s1", + "r1", + 27957.68284808907 + ], + [ + "s2", + null, + null + ] + ] + } + }, + { + "name": "join points within a distance", + "tool": "spatialJoin", + "params": { + "target": "points", + "join": "roads", + "match": "Within a distance", + "distance": 30000.0, + "units": "Meters", + "distanceField": "dist", + "operation": "Join one to many" + }, + "expect": { + "pairs": [ + [ + "a", + "r1", + 14020.411061039082 + ], + [ + "a", + "r2", + 28900.35711457829 + ], + [ + "b", + "r2", + 0.0 + ], + [ + "b", + "r1", + 14088.731617278729 + ], + [ + "c", + "r2", + 6949.6925145955665 + ], + [ + "c", + "r1", + 27609.423470296257 + ], + [ + "d", + "r3", + 5.591934318966577 + ], + [ + "e", + "r3", + 7892.016604980594 + ], + [ + "f", + null, + null + ], + [ + "g", + null, + null + ], + [ + "h", + null, + null + ], + [ + "i", + "r4", + 27798.630819739807 + ], + [ + "j", + "r6", + 6721.086468365455 + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/expected/clip_parcels.geojson b/crates/analysis/tests/reference/expected/clip_parcels.geojson new file mode 100644 index 000000000..19d71fd94 --- /dev/null +++ b/crates/analysis/tests/reference/expected/clip_parcels.geojson @@ -0,0 +1,227 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "P1" + }, + "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" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 46.37499999999999 + ], + [ + 10.875, + 46.37499999999999 + ], + [ + 10.875, + 46.50000000000001 + ], + [ + 11.0, + 46.50000000000001 + ], + [ + 11.0, + 45.999999999999986 + ], + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.500000000000002, + 46.37499999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 46.625 + ], + [ + 11.375, + 46.625 + ], + [ + 11.375, + 46.37499999999999 + ], + [ + 11.500000000000002, + 46.37499999999999 + ], + [ + 11.500000000000002, + 46.25 + ], + [ + 11.500000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P4" + }, + "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" + }, + "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" + }, + "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/clip_roads.geojson b/crates/analysis/tests/reference/expected/clip_roads.geojson new file mode 100644 index 000000000..e3ba16c2c --- /dev/null +++ b/crates/analysis/tests/reference/expected/clip_roads.geojson @@ -0,0 +1,81 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "r1" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.250000000000002, + 46.25 + ], + [ + 11.125, + 46.25 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r2" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.500000000000002, + 46.0625 + ], + [ + 10.500000000000002, + 46.37499999999999 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r3" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 11.5625, + 46.31250000000001 + ], + [ + 11.6875, + 46.4375 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r6" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 11.206763423774959, + 47.12499999999999 + ], + [ + 11.500000000000002, + 47.34499451815567 + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/expected/dissolve_all.geojson b/crates/analysis/tests/reference/expected/dissolve_all.geojson new file mode 100644 index 000000000..2f7d2a84b --- /dev/null +++ b/crates/analysis/tests/reference/expected/dissolve_all.geojson @@ -0,0 +1,230 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "all" + }, + "geometry": { + "type": "MultiPolygon", + "coordinates": [ + [ + [ + [ + 10.250000000000002, + 46.50000000000001 + ], + [ + 10.250000000000002, + 46.74999999999999 + ], + [ + 10.75, + 46.74999999999999 + ], + [ + 10.75, + 46.50000000000001 + ], + [ + 11.0, + 46.50000000000001 + ], + [ + 11.0, + 45.999999999999986 + ], + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.0, + 45.999999999999986 + ], + [ + 10.0, + 46.50000000000001 + ], + [ + 10.250000000000002, + 46.50000000000001 + ] + ] + ], + [ + [ + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 46.74999999999999 + ] + ], + [ + [ + 11.500000000000002, + 46.25 + ], + [ + 11.749999999999998, + 46.25 + ], + [ + 11.749999999999998, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.25 + ] + ] + ], + [ + [ + [ + 12.749999999999998, + 46.25 + ], + [ + 13.250000000000002, + 46.25 + ], + [ + 13.250000000000002, + 45.999999999999986 + ], + [ + 12.499999999999998, + 45.999999999999986 + ], + [ + 12.499999999999998, + 46.74999999999999 + ], + [ + 12.749999999999998, + 46.74999999999999 + ], + [ + 12.749999999999998, + 46.25 + ] + ] + ], + [ + [ + [ + 12.250000000000002, + 47.0 + ], + [ + 12.250000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 47.0 + ], + [ + 12.250000000000002, + 47.0 + ] + ] + ], + [ + [ + [ + 10.250000000000002, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.0 + ], + [ + 10.0, + 47.0 + ], + [ + 10.0, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.250000000000014 + ] + ] + ], + [ + [ + [ + 10.75, + 47.0 + ], + [ + 10.500000000000002, + 47.0 + ], + [ + 10.500000000000002, + 47.250000000000014 + ], + [ + 10.75, + 47.250000000000014 + ], + [ + 10.75, + 47.0 + ] + ] + ], + [ + [ + [ + 11.250000000000002, + 47.5 + ], + [ + 11.749999999999998, + 47.0 + ], + [ + 11.0, + 47.0 + ], + [ + 11.250000000000002, + 47.5 + ] + ] + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/expected/dissolve_code.geojson b/crates/analysis/tests/reference/expected/dissolve_code.geojson new file mode 100644 index 000000000..49fd11265 --- /dev/null +++ b/crates/analysis/tests/reference/expected/dissolve_code.geojson @@ -0,0 +1,259 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "A" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.0, + 45.999999999999986 + ], + [ + 10.0, + 46.50000000000001 + ], + [ + 10.250000000000002, + 46.50000000000001 + ], + [ + 10.250000000000002, + 46.74999999999999 + ], + [ + 10.75, + 46.74999999999999 + ], + [ + 10.75, + 46.50000000000001 + ], + [ + 11.0, + 46.50000000000001 + ], + [ + 11.0, + 45.999999999999986 + ], + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.0, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "B" + }, + "geometry": { + "type": "MultiPolygon", + "coordinates": [ + [ + [ + [ + 11.250000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 45.999999999999986 + ] + ], + [ + [ + 11.749999999999998, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.50000000000001 + ], + [ + 11.500000000000002, + 46.25 + ], + [ + 11.749999999999998, + 46.25 + ], + [ + 11.749999999999998, + 46.50000000000001 + ] + ] + ], + [ + [ + [ + 10.0, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.250000000000014 + ], + [ + 10.250000000000002, + 47.0 + ], + [ + 10.0, + 47.0 + ], + [ + 10.0, + 47.250000000000014 + ] + ] + ], + [ + [ + [ + 10.500000000000002, + 47.250000000000014 + ], + [ + 10.75, + 47.250000000000014 + ], + [ + 10.75, + 47.0 + ], + [ + 10.500000000000002, + 47.0 + ], + [ + 10.500000000000002, + 47.250000000000014 + ] + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "C" + }, + "geometry": { + "type": "MultiPolygon", + "coordinates": [ + [ + [ + [ + 12.499999999999998, + 46.74999999999999 + ], + [ + 12.749999999999998, + 46.74999999999999 + ], + [ + 12.749999999999998, + 46.25 + ], + [ + 13.250000000000002, + 46.25 + ], + [ + 13.250000000000002, + 45.999999999999986 + ], + [ + 12.499999999999998, + 45.999999999999986 + ], + [ + 12.499999999999998, + 46.74999999999999 + ] + ] + ], + [ + [ + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 47.0 + ], + [ + 12.250000000000002, + 47.0 + ], + [ + 12.250000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 46.74999999999999 + ] + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "D" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.0, + 47.0 + ], + [ + 11.749999999999998, + 47.0 + ], + [ + 11.250000000000002, + 47.5 + ], + [ + 11.0, + 47.0 + ] + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/expected/erase_parcels.geojson b/crates/analysis/tests/reference/expected/erase_parcels.geojson new file mode 100644 index 000000000..a030906dc --- /dev/null +++ b/crates/analysis/tests/reference/expected/erase_parcels.geojson @@ -0,0 +1,325 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "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" + }, + "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" + }, + "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" + }, + "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, + 46.74999999999999 + ], + [ + 12.749999999999998, + 46.74999999999999 + ], + [ + 12.749999999999998, + 46.25 + ], + [ + 13.250000000000002, + 46.25 + ], + [ + 13.250000000000002, + 45.999999999999986 + ], + [ + 12.499999999999998, + 45.999999999999986 + ], + [ + 12.499999999999998, + 46.74999999999999 + ] + ] + ] + } + }, + { + "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.000000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 47.0 + ], + [ + 12.250000000000002, + 47.0 + ], + [ + 12.250000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 46.74999999999999 + ] + ] + ] + } + }, + { + "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/hull_all.geojson b/crates/analysis/tests/reference/expected/hull_all.geojson new file mode 100644 index 000000000..9a868127d --- /dev/null +++ b/crates/analysis/tests/reference/expected/hull_all.geojson @@ -0,0 +1,46 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "all" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.0, + 45.999999999999986 + ], + [ + 10.0, + 47.250000000000014 + ], + [ + 11.250000000000002, + 47.5 + ], + [ + 12.749999999999998, + 46.74999999999999 + ], + [ + 13.250000000000002, + 46.25 + ], + [ + 13.250000000000002, + 45.999999999999986 + ], + [ + 10.0, + 45.999999999999986 + ] + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/expected/hull_each.geojson b/crates/analysis/tests/reference/expected/hull_each.geojson new file mode 100644 index 000000000..fdafcbda6 --- /dev/null +++ b/crates/analysis/tests/reference/expected/hull_each.geojson @@ -0,0 +1,269 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "P1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.0, + 45.999999999999986 + ], + [ + 10.0, + 46.50000000000001 + ], + [ + 10.500000000000002, + 46.50000000000001 + ], + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.0, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.500000000000002, + 45.999999999999986 + ], + [ + 10.500000000000002, + 46.50000000000001 + ], + [ + 11.0, + 46.50000000000001 + ], + [ + 11.0, + 45.999999999999986 + ], + [ + 10.500000000000002, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.250000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 45.999999999999986 + ], + [ + 11.250000000000002, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P4" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.0, + 47.0 + ], + [ + 10.0, + 47.250000000000014 + ], + [ + 10.75, + 47.250000000000014 + ], + [ + 10.75, + 47.0 + ], + [ + 10.0, + 47.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P5" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 12.499999999999998, + 45.999999999999986 + ], + [ + 12.499999999999998, + 46.74999999999999 + ], + [ + 12.749999999999998, + 46.74999999999999 + ], + [ + 13.250000000000002, + 46.25 + ], + [ + 13.250000000000002, + 45.999999999999986 + ], + [ + 12.499999999999998, + 45.999999999999986 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P6" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.250000000000002, + 46.25 + ], + [ + 10.250000000000002, + 46.74999999999999 + ], + [ + 10.75, + 46.74999999999999 + ], + [ + 10.75, + 46.25 + ], + [ + 10.250000000000002, + 46.25 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P7" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 12.000000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 47.0 + ], + [ + 12.250000000000002, + 47.0 + ], + [ + 12.250000000000002, + 46.74999999999999 + ], + [ + 12.000000000000002, + 46.74999999999999 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P8" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.0, + 47.0 + ], + [ + 11.250000000000002, + 47.5 + ], + [ + 11.749999999999998, + 47.0 + ], + [ + 11.0, + 47.0 + ] + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/generate.py b/crates/analysis/tests/reference/generate.py new file mode 100644 index 000000000..95f908d61 --- /dev/null +++ b/crates/analysis/tests/reference/generate.py @@ -0,0 +1,337 @@ +# /// script +# requires-python = ">=3.11" +# dependencies = ["shapely==2.1.*"] +# /// +"""Writes the reference layers and the answers RustGIS's tools must give for them. + +The answers come from GEOS (through Shapely), the geometry engine of QGIS and PostGIS. Run it +from this directory with `uv run generate.py`; the Rust test `reference.rs` reads what it writes. + +Overlays and spatial predicates are worked out in Web Mercator, the plane RustGIS's tools use, so +the comparison tests the geometry engine rather than the choice of plane. Distances are on a +sphere of radius 6,371,008.8 m along great circles, as RustGIS measures them, worked out exactly +from unit vectors (and checked against dense sampling). +""" + +import json +import math +from pathlib import Path + +import shapely +from shapely.geometry import LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, box, mapping +from shapely.ops import transform, unary_union + +HERE = Path(__file__).parent +MERC_R = 6378137.0 +MAX_LAT = 85.051_128_779_806_59 +EARTH_R = 6_371_008.8 + + +# --- the layers ------------------------------------------------------------------------------- +# Coordinates are multiples of 1/16°, so "on the edge" is exactly on the edge in either plane. + +def ring(*pts): + return list(pts) + + +PARCELS = [ + ("P1", "A", box(10, 46, 10.5, 46.5)), + # Clockwise, unlike the others: tools must not depend on the winding. + ("P2", "A", Polygon([(10.5, 46), (10.5, 46.5), (11, 46.5), (11, 46), (10.5, 46)])), + ("P3", "B", Polygon(box(11.25, 46, 12, 46.75).exterior.coords, [box(11.5, 46.25, 11.75, 46.5).exterior.coords])), + ("P4", "B", MultiPolygon([box(10, 47, 10.25, 47.25), box(10.5, 47, 10.75, 47.25)])), + ("P5", "C", Polygon([(12.5, 46), (13.25, 46), (13.25, 46.25), (12.75, 46.25), (12.75, 46.75), (12.5, 46.75)])), + ("P6", "A", box(10.25, 46.25, 10.75, 46.75)), + ("P7", "C", box(12, 46.75, 12.25, 47)), + ("P8", "D", Polygon([(11, 47), (11.75, 47), (11.25, 47.5)])), +] +ZONES = [ + ("Z1", box(10.25, 45.75, 11.5, 46.375)), + ("Z2", box(9.75, 46.875, 10.5, 47.5)), + ("Z3", box(13.5, 46, 14, 46.5)), + ("Z4", box(11.5, 46.25, 11.75, 46.5)), + ("Z5", box(11.125, 47.125, 11.5, 47.75)), + ("Z6", box(10.875, 46.125, 11.375, 46.625)), +] +POINTS = [ + ("a", Point(10.125, 46.125)), + ("b", Point(10.5, 46.125)), + ("c", Point(10.5, 46.5)), + ("d", Point(11.625, 46.375)), + ("e", Point(11.5, 46.375)), + ("f", Point(14.5, 48)), + ("g", Point(12.625, 46.125)), + ("h", Point(12.625, 46.125)), + ("i", Point(12, 46.75)), + ("j", Point(11.25, 47.25)), +] +ROADS = [ + ("r1", LineString([(9.875, 46.25), (11.125, 46.25)])), + ("r2", LineString([(10.5, 46.0625), (10.5, 46.4375)])), + ("r3", LineString([(11.5625, 46.3125), (11.6875, 46.4375)])), + ("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)])), +] +STATIONS = [ + ("s1", Point(10.75, 46.0)), + ("s2", Point(12.25, 47.25)), +] + + +def write_layer(name, rows, fields): + feats = [ + {"type": "Feature", "properties": dict(zip(fields, row[:-1])), "geometry": mapping(row[-1])} + for row in rows + ] + (HERE / "layers").mkdir(exist_ok=True) + text = json.dumps({"type": "FeatureCollection", "features": feats}, indent=1) + (HERE / "layers" / f"{name}.geojson").write_text(text + "\n") + + +# --- planes and distances --------------------------------------------------------------------- + +def to_merc(x, y, z=None): + lat = min(max(y, -MAX_LAT), MAX_LAT) + return MERC_R * math.radians(x), MERC_R * math.log(math.tan(math.pi / 4 + math.radians(lat) / 2)) + + +def to_lonlat(x, y, z=None): + return math.degrees(x / MERC_R), math.degrees(2 * math.atan(math.exp(y / MERC_R)) - math.pi / 2) + + +def m(g): + return transform(to_merc, g) + + +def ll(g): + return transform(to_lonlat, g) + + +def unit(lon, lat): + lo, la = math.radians(lon), math.radians(lat) + return (math.cos(la) * math.cos(lo), math.cos(la) * math.sin(lo), math.sin(la)) + + +def cross(a, b): + return (a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]) + + +def dot(a, b): + return sum(x * y for x, y in zip(a, b)) + + +def norm(a): + return math.sqrt(dot(a, a)) + + +def angle(a, b): + return math.atan2(norm(cross(a, b)), dot(a, b)) + + +def point_arc_m(p, a, b): + """Distance on the sphere from point p to the great-circle arc a→b (lon/lat pairs).""" + P, A, B = unit(*p), unit(*a), unit(*b) + n = cross(A, B) + if norm(n) < 1e-15: + return EARTH_R * angle(P, A) + n = tuple(x / norm(n) for x in n) + s = dot(P, n) + C = tuple(pi - s * ni for pi, ni in zip(P, n)) + if norm(C) > 1e-15 and dot(cross(A, C), n) >= 0 and dot(cross(C, B), n) >= 0: + return EARTH_R * angle(P, C) + return EARTH_R * min(angle(P, A), angle(P, B)) + + +def segments(g): + if isinstance(g, (Polygon, MultiPolygon)): + polys = [g] if isinstance(g, Polygon) else list(g.geoms) + for poly in polys: + for r in [poly.exterior, *poly.interiors]: + cs = list(r.coords) + yield from zip(cs, cs[1:]) + elif isinstance(g, (LineString, MultiLineString)): + lines = [g] if isinstance(g, LineString) else list(g.geoms) + for line in lines: + cs = list(line.coords) + yield from zip(cs, cs[1:]) + + +def point_distance_m(p, g): + """Ground distance from point p to g: zero when they intersect (in the tools' plane).""" + if isinstance(g, Point): + return EARTH_R * angle(unit(p.x, p.y), unit(g.x, g.y)) + if m(g).intersects(m(p)): + return 0.0 + return min(point_arc_m((p.x, p.y), a, b) for a, b in segments(g)) + + +def sampled_distance_m(p, g, step_m=2.0): + """The same distance by sampling every edge densely: a check on point_distance_m.""" + best = math.inf + P = unit(p.x, p.y) + for a, b in segments(g): + A, B = unit(*a), unit(*b) + w = angle(A, B) + n = max(1, math.ceil(w * EARTH_R / step_m)) + for k in range(n + 1): + t = k / n + s0, s1 = math.sin((1 - t) * w), math.sin(t * w) + q = tuple((s0 * x + s1 * y) / math.sin(w) for x, y in zip(A, B)) if w > 0 else A + best = min(best, EARTH_R * angle(P, q)) + return best + + +# --- the answers ------------------------------------------------------------------------------ + +def polygonal(g): + parts = [p for p in shapely.get_parts(g) if isinstance(p, (Polygon, MultiPolygon)) and p.area > 0] + return unary_union(parts) if parts else None + + +def linear(g): + parts = [p for p in shapely.get_parts(shapely.line_merge(g) if g.geom_type.endswith("LineString") else g) + if isinstance(p, (LineString, MultiLineString)) and p.length > 0] + return shapely.line_merge(MultiLineString([q for p in parts for q in shapely.get_parts(p)])) if parts else None + + +def write_expected(name, rows): + feats = [{"type": "Feature", "properties": {"id": i}, "geometry": mapping(g)} for i, g in rows] + (HERE / "expected").mkdir(exist_ok=True) + text = json.dumps({"type": "FeatureCollection", "features": feats}, indent=1) + (HERE / "expected" / f"{name}.geojson").write_text(text + "\n") + return f"expected/{name}.geojson" + + +def overlay(name, rows, op): + out = [] + for i, g in rows: + r = op(m(g)) + if r is not None and not r.is_empty: + out.append((i, ll(r))) + return write_expected(name, out) + + +def main(): + write_layer("parcels", PARCELS, ["id", "code"]) + write_layer("zones", ZONES, ["id"]) + write_layer("points", POINTS, ["id"]) + write_layer("roads", ROADS, ["id"]) + write_layer("stations", STATIONS, ["id"]) + + zones_m = unary_union([m(g) for _, g in ZONES]) + parcels = [(i, g) for i, _, g in PARCELS] + cases = [] + + def case(name, tool, params, expect): + cases.append({"name": name, "tool": tool, "params": params, "expect": expect}) + + # Clip and Erase: each input feature's part inside (outside) the union of the zones; features + # left with nothing (or only a shared edge or corner) are dropped. + case("clip polygons", "clip", {"input": "parcels", "clip": "zones"}, + {"features": overlay("clip_parcels", parcels, lambda g: polygonal(g.intersection(zones_m)))}) + case("clip lines", "clip", {"input": "roads", "clip": "zones"}, + {"features": overlay("clip_roads", ROADS, lambda g: linear(g.intersection(zones_m)))}) + case("clip points", "clip", {"input": "points", "clip": "zones"}, + {"ids": [i for i, g in POINTS if zones_m.intersects(m(g))]}) + case("erase", "erase", {"input": "parcels", "erase": "zones"}, + {"features": overlay("erase_parcels", parcels, lambda g: polygonal(g.difference(zones_m)))}) + + # 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] + case("dissolve by field", "dissolve", {"input": "parcels", "field": "code"}, + {"features": write_expected("dissolve_code", rows), "key": "code", + "counts": {c: sum(1 for _, k, _ in PARCELS if k == c) for c in codes}}) + case("dissolve all", "dissolve", {"input": "parcels"}, + {"features": write_expected("dissolve_all", [("all", ll(unary_union([m(g) for _, g in parcels])))]), + "key": None}) + + # Minimum Bounding Geometry: convex hulls in the tools' plane. + case("convex hull of each", "convexHull", {"input": "parcels", "group": "None"}, + {"features": overlay("hull_each", parcels, lambda g: g.convex_hull)}) + case("convex hull of all", "convexHull", {"input": "parcels", "group": "All"}, + {"features": write_expected("hull_all", [("all", ll(MultiPolygon([m(g) for _, g in parcels if g.geom_type == "Polygon"] + + [p for _, g in parcels if g.geom_type == "MultiPolygon" for p in m(g).geoms]).convex_hull))]), + "key": None}) + + # Select Layer By Location: features that intersect any selecting feature (edges count). + layers = {"parcels": parcels, "zones": ZONES, "points": POINTS, "roads": ROADS} + for inp, sel in [("parcels", "zones"), ("roads", "parcels"), ("points", "parcels"), ("roads", "zones")]: + sm = [m(g) for _, g in layers[sel]] + hit = [i for i, g in layers[inp] if any(s.intersects(m(g)) for s in sm)] + case(f"select {inp} by {sel}", "selectByLocation", {"input": inp, "select": sel}, {"ids": hit}) + if inp == "parcels": + miss = [i for i, _ in layers[inp] if i not in hit] + case(f"select {inp} by {sel}, inverted", "selectByLocation", {"input": inp, "select": sel, "invert": "Yes"}, + {"ids": miss}) + + # Summarize Within: points that intersect each polygon (edges count; duplicates count twice). + case("summarize within", "countPoints", {"polygons": "parcels", "points": "points"}, + {"field": "point_count", "values": {i: sum(1 for _, p in POINTS if m(g).intersects(m(p))) for i, g in parcels}}) + + # Spatial Join, one to many, keeping every target: the set of (target, join) pairs. + def pairs(targets, joins, rel): + out = [] + for ti, tg in targets: + hits = [ji for ji, jg in joins if rel(m(tg), m(jg))] + out += [[ti, ji] for ji in hits] or [[ti, None]] + return out + + many = {"operation": "Join one to many"} + case("join intersects", "spatialJoin", {"target": "parcels", "join": "zones", "match": "Intersects", **many}, + {"pairs": pairs(parcels, ZONES, lambda t, j: t.intersects(j))}) + case("join contains", "spatialJoin", {"target": "parcels", "join": "points", "match": "Contains", **many}, + {"pairs": pairs(parcels, POINTS, lambda t, j: t.contains(j))}) + case("join within", "spatialJoin", {"target": "points", "join": "parcels", "match": "Within", **many}, + {"pairs": pairs(POINTS, parcels, lambda t, j: t.within(j))}) + case("join roads intersect parcels", "spatialJoin", {"target": "roads", "join": "parcels", "match": "Intersects", **many}, + {"pairs": pairs(ROADS, parcels, lambda t, j: t.intersects(j))}) + + # Closest and Within a distance: ground distances, checked against dense sampling. + def dist(p, g): + d = point_distance_m(p, g) + if d > 0: + # Sampling every 2 m overestimates by at most about step² / 8d. + s = sampled_distance_m(p, g) + assert -1e-6 < s - d < 1e-3 + 0.5 / d, (p, d, s) + return d + + rows = [] + for ti, tp in POINTS: + ds = [(dist(tp, g), ji) for ji, g in parcels] + best = min(d for d, _ in ds) + ji = next(j for d, j in ds if d == best) # ties keep the join layer's order + rows.append([ti, ji, best]) + case("join closest", "spatialJoin", + {"target": "points", "join": "parcels", "match": "Closest", "units": "Meters", "distanceField": "dist", **many}, + {"pairs": rows}) + + limit = 30_000.0 + rows = [] + for ti, tp in STATIONS + POINTS: + ds = sorted((dist(tp, g), ji) for ji, g in ROADS) + assert all(abs(d - limit) > 100 for d, _ in ds), "a distance too near the limit" + hits = [[ti, ji, d] for d, ji in ds if d <= limit] + rows += hits or [[ti, None, None]] + case("join within a distance", "spatialJoin", + {"target": "stations", "join": "roads", "match": "Within a distance", "distance": limit, "units": "Meters", + "distanceField": "dist", **many}, + {"pairs": [r for r in rows if r[0].startswith("s")]}) + case("join points within a distance", "spatialJoin", + {"target": "points", "join": "roads", "match": "Within a distance", "distance": limit, "units": "Meters", + "distanceField": "dist", **many}, + {"pairs": [r for r in rows if not r[0].startswith("s")]}) + + meta = { + "about": "Answers from GEOS; regenerate with `uv run generate.py` (see its header).", + "shapely": shapely.__version__, + "geos": shapely.geos_version_string, + "cases": cases, + } + (HERE / "cases.json").write_text(json.dumps(meta, indent=1) + "\n") + print(f"{len(cases)} cases") + + +if __name__ == "__main__": + main() diff --git a/crates/analysis/tests/reference/layers/parcels.geojson b/crates/analysis/tests/reference/layers/parcels.geojson new file mode 100644 index 000000000..3d8a04fdf --- /dev/null +++ b/crates/analysis/tests/reference/layers/parcels.geojson @@ -0,0 +1,329 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "P1", + "code": "A" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.5, + 46.0 + ], + [ + 10.5, + 46.5 + ], + [ + 10.0, + 46.5 + ], + [ + 10.0, + 46.0 + ], + [ + 10.5, + 46.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P2", + "code": "A" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.5, + 46.0 + ], + [ + 10.5, + 46.5 + ], + [ + 11.0, + 46.5 + ], + [ + 11.0, + 46.0 + ], + [ + 10.5, + 46.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P3", + "code": "B" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 12.0, + 46.0 + ], + [ + 12.0, + 46.75 + ], + [ + 11.25, + 46.75 + ], + [ + 11.25, + 46.0 + ], + [ + 12.0, + 46.0 + ] + ], + [ + [ + 11.75, + 46.25 + ], + [ + 11.75, + 46.5 + ], + [ + 11.5, + 46.5 + ], + [ + 11.5, + 46.25 + ], + [ + 11.75, + 46.25 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P4", + "code": "B" + }, + "geometry": { + "type": "MultiPolygon", + "coordinates": [ + [ + [ + [ + 10.25, + 47.0 + ], + [ + 10.25, + 47.25 + ], + [ + 10.0, + 47.25 + ], + [ + 10.0, + 47.0 + ], + [ + 10.25, + 47.0 + ] + ] + ], + [ + [ + [ + 10.75, + 47.0 + ], + [ + 10.75, + 47.25 + ], + [ + 10.5, + 47.25 + ], + [ + 10.5, + 47.0 + ], + [ + 10.75, + 47.0 + ] + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P5", + "code": "C" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 12.5, + 46.0 + ], + [ + 13.25, + 46.0 + ], + [ + 13.25, + 46.25 + ], + [ + 12.75, + 46.25 + ], + [ + 12.75, + 46.75 + ], + [ + 12.5, + 46.75 + ], + [ + 12.5, + 46.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P6", + "code": "A" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.75, + 46.25 + ], + [ + 10.75, + 46.75 + ], + [ + 10.25, + 46.75 + ], + [ + 10.25, + 46.25 + ], + [ + 10.75, + 46.25 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P7", + "code": "C" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 12.25, + 46.75 + ], + [ + 12.25, + 47.0 + ], + [ + 12.0, + 47.0 + ], + [ + 12.0, + 46.75 + ], + [ + 12.25, + 46.75 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "P8", + "code": "D" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.0, + 47.0 + ], + [ + 11.75, + 47.0 + ], + [ + 11.25, + 47.5 + ], + [ + 11.0, + 47.0 + ] + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/layers/points.geojson b/crates/analysis/tests/reference/layers/points.geojson new file mode 100644 index 000000000..61b565d3e --- /dev/null +++ b/crates/analysis/tests/reference/layers/points.geojson @@ -0,0 +1,135 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "a" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 10.125, + 46.125 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "b" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 10.5, + 46.125 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "c" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 10.5, + 46.5 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "d" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 11.625, + 46.375 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "e" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 11.5, + 46.375 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "f" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 14.5, + 48.0 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "g" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 12.625, + 46.125 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "h" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 12.625, + 46.125 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "i" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 12.0, + 46.75 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "j" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 11.25, + 47.25 + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/layers/roads.geojson b/crates/analysis/tests/reference/layers/roads.geojson new file mode 100644 index 000000000..c7d2d2056 --- /dev/null +++ b/crates/analysis/tests/reference/layers/roads.geojson @@ -0,0 +1,119 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "r1" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 9.875, + 46.25 + ], + [ + 11.125, + 46.25 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r2" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.5, + 46.0625 + ], + [ + 10.5, + 46.4375 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r3" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 11.5625, + 46.3125 + ], + [ + 11.6875, + 46.4375 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r4" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 12.0, + 46.5 + ], + [ + 12.5, + 46.5 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r5" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 15.0, + 49.0 + ], + [ + 15.5, + 49.5 + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "r6" + }, + "geometry": { + "type": "LineString", + "coordinates": [ + [ + 10.875, + 46.875 + ], + [ + 11.875, + 47.625 + ] + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/layers/stations.geojson b/crates/analysis/tests/reference/layers/stations.geojson new file mode 100644 index 000000000..7e3cceef0 --- /dev/null +++ b/crates/analysis/tests/reference/layers/stations.geojson @@ -0,0 +1,31 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "s1" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 10.75, + 46.0 + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "s2" + }, + "geometry": { + "type": "Point", + "coordinates": [ + 12.25, + 47.25 + ] + } + } + ] +} diff --git a/crates/analysis/tests/reference/layers/zones.geojson b/crates/analysis/tests/reference/layers/zones.geojson new file mode 100644 index 000000000..0c5d6bea0 --- /dev/null +++ b/crates/analysis/tests/reference/layers/zones.geojson @@ -0,0 +1,203 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": { + "id": "Z1" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.5, + 45.75 + ], + [ + 11.5, + 46.375 + ], + [ + 10.25, + 46.375 + ], + [ + 10.25, + 45.75 + ], + [ + 11.5, + 45.75 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "Z2" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 10.5, + 46.875 + ], + [ + 10.5, + 47.5 + ], + [ + 9.75, + 47.5 + ], + [ + 9.75, + 46.875 + ], + [ + 10.5, + 46.875 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "Z3" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 14.0, + 46.0 + ], + [ + 14.0, + 46.5 + ], + [ + 13.5, + 46.5 + ], + [ + 13.5, + 46.0 + ], + [ + 14.0, + 46.0 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "Z4" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.75, + 46.25 + ], + [ + 11.75, + 46.5 + ], + [ + 11.5, + 46.5 + ], + [ + 11.5, + 46.25 + ], + [ + 11.75, + 46.25 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "Z5" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.5, + 47.125 + ], + [ + 11.5, + 47.75 + ], + [ + 11.125, + 47.75 + ], + [ + 11.125, + 47.125 + ], + [ + 11.5, + 47.125 + ] + ] + ] + } + }, + { + "type": "Feature", + "properties": { + "id": "Z6" + }, + "geometry": { + "type": "Polygon", + "coordinates": [ + [ + [ + 11.375, + 46.125 + ], + [ + 11.375, + 46.625 + ], + [ + 10.875, + 46.625 + ], + [ + 10.875, + 46.125 + ], + [ + 11.375, + 46.125 + ] + ] + ] + } + } + ] +} From e726a5b3605a0ee116f0d3e71c572827e38a53bb Mon Sep 17 00:00:00 2001 From: mavrick Date: Sun, 25 Oct 2026 21:17:27 -0500 Subject: [PATCH 2/2] test(analysis): measure reference distances on WGS 84; review follow-ups #905 moved Spatial Join distances to the WGS 84 ellipsoid, and the first CI run caught it: distances up to 0.3 % from the spherical reference. The reference now measures with GeographicLib (pyproj): point to point directly, and to an edge as the nearest point of its geodesic, found by searching along it and checked against dense sampling. Review follow-ups: Summarize Within's features are compared once each (missing, extra or repeated ones fail); an output with another kind of geometry than the reference fails; dissolve counts without a key match by_key; a selection in an unknown layer is reported, not a panic; the distance checks are a single match. Co-Authored-By: Claude Opus 5.5 --- ATTRIBUTION.md | 2 +- crates/analysis/tests/reference.rs | 70 +++++++++++++---- crates/analysis/tests/reference/cases.json | 28 ++++--- crates/analysis/tests/reference/generate.py | 87 ++++++++++----------- 4 files changed, 111 insertions(+), 76 deletions(-) diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index d444337d3..1f4313fe9 100644 --- a/ATTRIBUTION.md +++ b/ATTRIBUTION.md @@ -13,7 +13,7 @@ which every package ships next to LICENSE. They are generated from `Cargo.lock` | Sample data (`countries`, `countries_chn`, `cities`, `rivers`, `lakes` GeoJSON), embedded at build time from [rustgis-assets](https://github.com/opengeos/rustgis-assets/tree/main/data) | Natural Earth 1:110m, attributes trimmed; `countries_chn` is Natural Earth's China point of view (1:10m, simplified) | Public domain | | QGIS screenshot (`qgis-world-atlas.webp` in [rustgis-assets](https://github.com/opengeos/rustgis-assets/tree/main/images)), on the [RustGIS and QGIS](docs/comparison/qgis.md) page | Taken by the RustGIS project: [QGIS](https://qgis.org) 4.2 showing the Natural Earth sample data above | QGIS is GPL-2.0-or-later; the data is public domain | | QGIS test corpus (`crates/map/testdata/qgis/generated/`: styles, projects and their sample data) | Dominik Strebel ([@pick2510](https://github.com/pick2510)) for the RustGIS project: written by [QGIS](https://qgis.org) 4.2.3 through its API, from `generate.py` there; the data is made up | Apache-2.0 | -| Analysis reference corpus (`crates/analysis/tests/reference/`: layers and expected results) | Written for the RustGIS project from `generate.py` there: the layers are made up; the expected results were computed by [GEOS](https://libgeos.org) 3.13.1 through [Shapely](https://shapely.readthedocs.io) 2.1.2, run only to write them (GEOS is LGPL-2.1, Shapely BSD-3-Clause; neither is included in or linked into RustGIS) | Apache-2.0 | +| Analysis reference corpus (`crates/analysis/tests/reference/`: layers and expected results) | Written for the RustGIS project from `generate.py` there: the layers are made up; the expected results were computed by [GEOS](https://libgeos.org) 3.13.1 through [Shapely](https://shapely.readthedocs.io) 2.1.2, and distances by GeographicLib in [PROJ](https://proj.org) 9.5.1 through [pyproj](https://pyproj4.github.io/pyproj/) 3.7.2, run only to write them (GEOS is LGPL-2.1, Shapely BSD-3-Clause, PROJ and pyproj MIT; none is included in or linked into RustGIS) | Apache-2.0 | Raster coordinate systems are looked up by EPSG code in the PROJ definitions of the [crs-definitions](https://crates.io/crates/crs-definitions) crate (CC0-1.0), compiled into the diff --git a/crates/analysis/tests/reference.rs b/crates/analysis/tests/reference.rs index 4ef2e5506..5b28b96f8 100644 --- a/crates/analysis/tests/reference.rs +++ b/crates/analysis/tests/reference.rs @@ -18,9 +18,9 @@ use serde_json::Value as Json; const AREA_TOLERANCE: f64 = 1e-7; /// Lines may differ in length, and sit, this far (Web Mercator metres) from the reference. const POSITION_TOLERANCE_M: f64 = 1e-3; -/// Distances may differ by this fraction, plus a millimetre. The reference measures along great -/// circles; the join finds the nearest points in an azimuthal equidistant plane (see `join.rs`), -/// whose straight edges are within this of the great circles, as its own tests require. +/// Distances may differ by this fraction, plus a millimetre. The reference measures to the WGS 84 +/// geodesic between an edge's vertices; the join finds the nearest points in an azimuthal +/// equidistant plane (see `join.rs`), whose straight edges are within this of the geodesics. const DISTANCE_TOLERANCE: f64 = 1e-4; /// The reference corpus: the layers, the expected results and `cases.json`. @@ -88,8 +88,28 @@ fn lines(gs: &[Geometry]) -> MultiLineString { MultiLineString::new(gs.iter().filter_map(to_multiline).flat_map(|m| m.0).collect()) } +/// The kinds of geometry among `gs`, by name, in order and without repeats. +fn kinds(gs: &[Geometry]) -> Vec<&'static str> { + let mut out: Vec<&'static str> = gs + .iter() + .map(|g| match g { + Geometry::Point(_) | Geometry::MultiPoint(_) => "point", + Geometry::Line(_) => "line", + Geometry::Polygon(_) => "polygon", + }) + .collect(); + out.sort_unstable(); + out.dedup(); + out +} + /// How one geometry differs from the reference, if it does beyond the tolerances. fn compare_geometry(got: &[Geometry], want: &[Geometry]) -> Option { + // The comparisons below look only at the reference's kind of geometry, so another kind in the + // output would go unseen. + if kinds(got) != kinds(want) { + return Some(format!("{:?} geometry, the reference {:?}", kinds(got), kinds(want))); + } if want.iter().any(|g| matches!(g, Geometry::Polygon(_))) { let (g, w) = (polygons(got), polygons(want)); // The area the two don't share, not vertex positions: GEOS keeps vertices partway along a @@ -140,7 +160,11 @@ fn check(case: &Json, layers: &Layers) -> Vec { } if let Some(counts) = expect["counts"].as_object() { for (row, _) in ds.features.iter().enumerate() { - let k = key.and_then(|key| text(&ds, row, key)).unwrap_or_default(); + // 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 (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:?}")); @@ -150,10 +174,10 @@ fn check(case: &Json, layers: &Layers) -> Vec { } else if let Some(ids) = expect["ids"].as_array() { let want: Vec = ids.iter().filter_map(json_text).collect(); let got: Vec = match result.output { - Output::Selection { layer, rows } => { - let ds = &layers.0[&layer]; - rows.iter().filter_map(|r| text(ds, *r, "id")).collect() - } + Output::Selection { layer, rows } => match layers.0.get(&layer) { + Some(ds) => rows.iter().filter_map(|r| text(ds, *r, "id")).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(), None => return vec!["neither a selection nor a layer came out".into()], @@ -165,6 +189,14 @@ fn check(case: &Json, layers: &Layers) -> Vec { } else if let Some(values) = expect["values"].as_object() { let Some(ds) = layer_of(result.output) else { return vec!["no layer came out".into()] }; let field = expect["field"].as_str().unwrap_or_default(); + // Every feature once: a missing, extra or repeated one fails even when the counts agree. + let mut ids: Vec = (0..ds.features.len()).map(|row| text(&ds, row, "id").unwrap_or_default()).collect(); + let mut want_ids: Vec = values.keys().cloned().collect(); + ids.sort(); + want_ids.sort(); + if ids != want_ids { + problems.push(format!("features {ids:?}, the reference {want_ids:?}")); + } for row in 0..ds.features.len() { let id = text(&ds, row, "id").unwrap_or_default(); let (got, want) = (text(&ds, row, field), values.get(&id).and_then(json_text)); @@ -202,12 +234,16 @@ fn check(case: &Json, layers: &Layers) -> Vec { continue; } for ((j, gd), (_, wd)) in g.iter().zip(&w) { - if let (Some(gd), Some(wd)) = (gd, wd) - && (gd - wd).abs() > DISTANCE_TOLERANCE * wd + 1e-3 - { - problems.push(format!("{t}→{}: {gd:.4} m, the reference {wd:.4} m", j.clone().unwrap_or_default())); - } else if gd.is_some() != wd.is_some() { - problems.push(format!("{t}→{}: distance {gd:?}, the reference {wd:?}", j.clone().unwrap_or_default())); + let j = j.clone().unwrap_or_default(); + match (gd, wd) { + (Some(gd), Some(wd)) if (gd - wd).abs() > DISTANCE_TOLERANCE * wd + 1e-3 => { + problems.push(format!("{t}→{j}: {gd:.4} m, the reference {wd:.4} m")); + } + // A distance where the reference has none, or none where it has one. + (Some(_), None) | (None, Some(_)) => { + problems.push(format!("{t}→{j}: distance {gd:?}, the reference {wd:?}")); + } + _ => {} } } } @@ -222,7 +258,9 @@ fn check(case: &Json, layers: &Layers) -> Vec { #[test] fn tools_agree_with_the_reference() { - let cases: Json = serde_json::from_str(&std::fs::read_to_string(dir().join("cases.json")).unwrap()).unwrap(); + let path = dir().join("cases.json"); + let text = std::fs::read_to_string(&path).unwrap_or_else(|e| panic!("{}: {e}", path.display())); + let cases: Json = serde_json::from_str(&text).unwrap_or_else(|e| panic!("{}: {e}", path.display())); let layers = Layers( ["parcels", "zones", "points", "roads", "stations"] .into_iter() @@ -230,7 +268,7 @@ fn tools_agree_with_the_reference() { .collect(), ); let mut report = Vec::new(); - for case in cases["cases"].as_array().unwrap() { + for case in cases["cases"].as_array().expect("cases.json has a \"cases\" list") { let problems = check(case, &layers); if !problems.is_empty() { report.push(format!("{}:\n {}", case["name"].as_str().unwrap_or("?"), problems.join("\n "))); diff --git a/crates/analysis/tests/reference/cases.json b/crates/analysis/tests/reference/cases.json index 9783c70da..868fce829 100644 --- a/crates/analysis/tests/reference/cases.json +++ b/crates/analysis/tests/reference/cases.json @@ -2,6 +2,8 @@ "about": "Answers from GEOS; regenerate with `uv run generate.py` (see its header).", "shapely": "2.1.2", "geos": "3.13.1", + "pyproj": "3.7.2", + "proj": "9.5.1", "cases": [ { "name": "clip polygons", @@ -461,7 +463,7 @@ [ "d", "P3", - 9589.67458520921 + 9617.286477596233 ], [ "e", @@ -471,7 +473,7 @@ [ "f", "P5", - 191521.6332492596 + 191772.8154582311 ], [ "g", @@ -513,12 +515,12 @@ [ "s1", "r2", - 20512.8663375754 + 20563.88085989103 ], [ "s1", "r1", - 27957.68284808907 + 27947.80972153128 ], [ "s2", @@ -545,12 +547,12 @@ [ "a", "r1", - 14020.411061039082 + 14015.746772213406 ], [ "a", "r2", - 28900.35711457829 + 28983.146568103246 ], [ "b", @@ -560,27 +562,27 @@ [ "b", "r1", - 14088.731617278729 + 14084.263048419794 ], [ "c", "r2", - 6949.6925145955665 + 6947.529600616946 ], [ "c", "r1", - 27609.423470296257 + 27599.770334571844 ], [ "d", "r3", - 5.591934318966577 + 5.612923586871125 ], [ "e", "r3", - 7892.016604980594 + 7906.554524361339 ], [ "f", @@ -600,12 +602,12 @@ [ "i", "r4", - 27798.630819739807 + 27790.741782081575 ], [ "j", "r6", - 6721.086468365455 + 6730.820161190676 ] ] } diff --git a/crates/analysis/tests/reference/generate.py b/crates/analysis/tests/reference/generate.py index 95f908d61..20b574d17 100644 --- a/crates/analysis/tests/reference/generate.py +++ b/crates/analysis/tests/reference/generate.py @@ -1,6 +1,6 @@ # /// script # requires-python = ">=3.11" -# dependencies = ["shapely==2.1.*"] +# dependencies = ["shapely==2.1.*", "pyproj==3.7.*"] # /// """Writes the reference layers and the answers RustGIS's tools must give for them. @@ -8,23 +8,26 @@ from this directory with `uv run generate.py`; the Rust test `reference.rs` reads what it writes. Overlays and spatial predicates are worked out in Web Mercator, the plane RustGIS's tools use, so -the comparison tests the geometry engine rather than the choice of plane. Distances are on a -sphere of radius 6,371,008.8 m along great circles, as RustGIS measures them, worked out exactly -from unit vectors (and checked against dense sampling). +the comparison tests the geometry engine rather than the choice of plane. Distances are on the +WGS 84 ellipsoid, as RustGIS measures them: GeographicLib (through pyproj) gives the distance +between two points, and the nearest point of an edge, the geodesic between its vertices, is found +by searching along it (and checked against dense sampling). """ import json import math from pathlib import Path +import pyproj import shapely from shapely.geometry import LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, box, mapping +from pyproj import Geod from shapely.ops import transform, unary_union HERE = Path(__file__).parent MERC_R = 6378137.0 MAX_LAT = 85.051_128_779_806_59 -EARTH_R = 6_371_008.8 +GEOD = Geod(ellps="WGS84") # --- the layers ------------------------------------------------------------------------------- @@ -108,39 +111,34 @@ def ll(g): return transform(to_lonlat, g) -def unit(lon, lat): - lo, la = math.radians(lon), math.radians(lat) - return (math.cos(la) * math.cos(lo), math.cos(la) * math.sin(lo), math.sin(la)) - - -def cross(a, b): - return (a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]) - - -def dot(a, b): - return sum(x * y for x, y in zip(a, b)) - - -def norm(a): - return math.sqrt(dot(a, a)) - - -def angle(a, b): - return math.atan2(norm(cross(a, b)), dot(a, b)) +def geodesic_m(p, q): + """Distance on WGS 84 between two lon/lat points.""" + return GEOD.inv(p[0], p[1], q[0], q[1], return_back_azimuth=True)[2] def point_arc_m(p, a, b): - """Distance on the sphere from point p to the great-circle arc a→b (lon/lat pairs).""" - P, A, B = unit(*p), unit(*a), unit(*b) - n = cross(A, B) - if norm(n) < 1e-15: - return EARTH_R * angle(P, A) - n = tuple(x / norm(n) for x in n) - s = dot(P, n) - C = tuple(pi - s * ni for pi, ni in zip(P, n)) - if norm(C) > 1e-15 and dot(cross(A, C), n) >= 0 and dot(cross(C, B), n) >= 0: - return EARTH_R * angle(P, C) - return EARTH_R * min(angle(P, A), angle(P, B)) + """Distance on WGS 84 from point p to the geodesic a→b (lon/lat pairs): its nearest point is + bracketed by sampling along the geodesic, then found by golden-section search.""" + az, _, length = GEOD.inv(a[0], a[1], b[0], b[1], return_back_azimuth=True) + if length == 0: + return geodesic_m(p, a) + + def at(s): + lon, lat, _ = GEOD.fwd(a[0], a[1], az, s) + return geodesic_m(p, (lon, lat)) + + n = 400 + ds = [at(length * k / n) for k in range(n + 1)] + k = min(range(n + 1), key=ds.__getitem__) + lo, hi = length * max(k - 1, 0) / n, length * min(k + 1, n) / n + g = (math.sqrt(5) - 1) / 2 + while hi - lo > 1e-6: + x1, x2 = hi - g * (hi - lo), lo + g * (hi - lo) + if at(x1) <= at(x2): + hi = x2 + else: + lo = x1 + return min(at((lo + hi) / 2), ds[0], ds[-1]) def segments(g): @@ -160,25 +158,20 @@ def segments(g): def point_distance_m(p, g): """Ground distance from point p to g: zero when they intersect (in the tools' plane).""" if isinstance(g, Point): - return EARTH_R * angle(unit(p.x, p.y), unit(g.x, g.y)) + return geodesic_m((p.x, p.y), (g.x, g.y)) if m(g).intersects(m(p)): return 0.0 return min(point_arc_m((p.x, p.y), a, b) for a, b in segments(g)) def sampled_distance_m(p, g, step_m=2.0): - """The same distance by sampling every edge densely: a check on point_distance_m.""" + """The same distance by sampling every edge's geodesic densely: a check on point_distance_m.""" best = math.inf - P = unit(p.x, p.y) for a, b in segments(g): - A, B = unit(*a), unit(*b) - w = angle(A, B) - n = max(1, math.ceil(w * EARTH_R / step_m)) - for k in range(n + 1): - t = k / n - s0, s1 = math.sin((1 - t) * w), math.sin(t * w) - q = tuple((s0 * x + s1 * y) / math.sin(w) for x, y in zip(A, B)) if w > 0 else A - best = min(best, EARTH_R * angle(P, q)) + pts = GEOD.inv_intermediate(a[0], a[1], b[0], b[1], del_s=step_m, initial_idx=0, terminus_idx=0, return_back_azimuth=True) + lons, lats = list(pts.lons), list(pts.lats) + _, _, d = GEOD.inv([p.x] * len(lons), [p.y] * len(lats), lons, lats, return_back_azimuth=True) + best = min(best, *d) return best @@ -327,6 +320,8 @@ def dist(p, g): "about": "Answers from GEOS; regenerate with `uv run generate.py` (see its header).", "shapely": shapely.__version__, "geos": shapely.geos_version_string, + "pyproj": pyproj.__version__, + "proj": pyproj.proj_version_str, "cases": cases, } (HERE / "cases.json").write_text(json.dumps(meta, indent=1) + "\n")