test(analysis): check overlay, selection and join tools against GEOS

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

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<String, Dataset>);
+
+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<String> {
+    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<String> {
+    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<String, Vec<Geometry>> {
+    let mut out: BTreeMap<String, Vec<Geometry>> = 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<f64> {
+    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<f64> {
+    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<String> {
+    // 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<Dataset> {
+    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<String> {
+    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<String> = ids.iter().filter_map(json_text).collect();
+        let got: Vec<String> = 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<String> = (0..ds.features.len()).map(|row| text(&ds, row, "id").unwrap_or_default()).collect();
+        let mut want_ids: Vec<String> = 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<String, Vec<(Option<String>, Option<f64>)>>;
+        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<String>, Option<f64>)]| v.iter().map(|p| p.0.clone()).collect::<Vec<_>>();
+            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"));
+}