diff --git a/ATTRIBUTION.md b/ATTRIBUTION.md index 0f72db35e..1f4313fe9 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, 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 new file mode 100644 index 000000000..5b28b96f8 --- /dev/null +++ b/crates/analysis/tests/reference.rs @@ -0,0 +1,279 @@ +//! 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 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`. +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()) +} + +/// 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 + // 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() { + // 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:?}")); + } + } + } + } 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 } => 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()], + }, + }; + 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(); + // 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)); + 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) { + 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:?}")); + } + _ => {} + } + } + } + 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 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() + .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().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 "))); + } + } + 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..868fce829 --- /dev/null +++ b/crates/analysis/tests/reference/cases.json @@ -0,0 +1,616 @@ +{ + "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", + "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", + 9617.286477596233 + ], + [ + "e", + "P3", + 0.0 + ], + [ + "f", + "P5", + 191772.8154582311 + ], + [ + "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", + 20563.88085989103 + ], + [ + "s1", + "r1", + 27947.80972153128 + ], + [ + "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", + 14015.746772213406 + ], + [ + "a", + "r2", + 28983.146568103246 + ], + [ + "b", + "r2", + 0.0 + ], + [ + "b", + "r1", + 14084.263048419794 + ], + [ + "c", + "r2", + 6947.529600616946 + ], + [ + "c", + "r1", + 27599.770334571844 + ], + [ + "d", + "r3", + 5.612923586871125 + ], + [ + "e", + "r3", + 7906.554524361339 + ], + [ + "f", + null, + null + ], + [ + "g", + null, + null + ], + [ + "h", + null, + null + ], + [ + "i", + "r4", + 27790.741782081575 + ], + [ + "j", + "r6", + 6730.820161190676 + ] + ] + } + } + ] +} 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..20b574d17 --- /dev/null +++ b/crates/analysis/tests/reference/generate.py @@ -0,0 +1,332 @@ +# /// script +# requires-python = ">=3.11" +# dependencies = ["shapely==2.1.*", "pyproj==3.7.*"] +# /// +"""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 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 +GEOD = Geod(ellps="WGS84") + + +# --- 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 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 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): + 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 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's geodesic densely: a check on point_distance_m.""" + best = math.inf + for a, b in segments(g): + 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 + + +# --- 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, + "pyproj": pyproj.__version__, + "proj": pyproj.proj_version_str, + "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 + ] + ] + ] + } + } + ] +}