Analysis: buffers, Tissot's circles and Spatial Join distances on the WGS 84 ellipsoid (#888)

Cargo.toml

@@ -35,6 +35,8 @@ web-time = "1"
 
 # geo / formats / rendering
 geo = "0.33"
+# Geodesics on the WGS 84 ellipsoid (already in the tree through geo; opengeos/RustGIS#888)
+geographiclib-rs = { version = "0.2", default-features = false }
 # A5, the equal-area pentagonal DGGS (opengeos/RustGIS#551)
 a5 = "1"
 shapefile = { version = "0.9", features = ["encoding_rs", "yore"] }

crates/analysis/src/join.rs

@@ -2,12 +2,12 @@
 //! rules out features too far away to matter, the overlap a feature has with a polygon, and the
 //! rules that roll the values of several matches up into one.
 //!
-//! Distances are on the sphere buffers use. The nearest points of two features are found in an
-//! azimuthal equidistant plane, centred first between the features and then between the points
-//! found, and measured along the great circle between them.
+//! Distances are on the WGS 84 ellipsoid, as buffers' are. The nearest points of two features are
+//! found in an azimuthal equidistant plane, centred first between the features and then between
+//! the points found, and measured along the geodesic between them.
 
 use rustgis_data::{Coord, Dataset, FieldType, Geometry, Value};
-use rustgis_geom::{Aeqd, EARTH_RADIUS_M, Extent, geodesic_interpolate, geodesic_length_m, geodesic_ring_area_m2, haversine_m};
+use rustgis_geom::{Aeqd, Extent, WGS84_MIN_RADIUS_M, geodesic_interpolate, geodesic_length_m, geodesic_ring_area_m2, wgs84_distance_m};
 use serde_json::Value as Json;
 
 /// Whether a coordinate is a place on the Earth in lon/lat: finite, with a latitude of ±90° at most
@@ -39,7 +39,7 @@ pub fn distance_m(a: &Geometry, b: &Geometry) -> Option<f64> {
     let (fa, fb) = (first_coord(a)?, first_coord(b)?);
     let (p, q) = nearest_pair(a, b, &Aeqd::new(geodesic_interpolate(fa, fb, 0.5)), Some((fa, fb)))?;
     let (p, q) = nearest_pair(a, b, &Aeqd::new(geodesic_interpolate(p, q, 0.5)), Some((p, q))).unwrap_or((p, q));
-    let d = haversine_m(p, q);
+    let d = wgs84_distance_m(p, q);
     d.is_finite().then_some(d)
 }
 
@@ -131,10 +131,11 @@ fn closest_on(p: Coord, a: Coord, b: Coord) -> Coord {
 /// other can be when one is no further than `lat` degrees from the equator. The longitude is
 /// `None` when they can reach a pole, so any longitudes can be that close. Two points whose
 /// longitudes differ by Δλ, both no further from the equator than φ, are at least
-/// 2R·asin(cos φ·sin(Δλ/2)) apart.
+/// 2R·asin(cos φ·sin(Δλ/2)) apart, with R the ellipsoid's smallest radius of curvature, so the
+/// ellipsoid's shorter distances near the equator are never ruled out.
 fn spread(lat: f64, metres: f64) -> (f64, Option<f64>) {
     // A little slack, so a feature exactly at the distance isn't ruled out by rounding.
-    let r = (metres.max(0.0) * (1.0 + 1e-9) + 1e-3) / EARTH_RADIUS_M;
+    let r = (metres.max(0.0) * (1.0 + 1e-9) + 1e-3) / WGS84_MIN_RADIUS_M;
     let dlat = r.to_degrees();
     let phi = (lat.abs() + dlat).to_radians();
     let reach = (r / 2.0).min(std::f64::consts::FRAC_PI_2).sin() / phi.cos();
@@ -346,12 +347,32 @@ pub fn parse_merge(join: &Dataset, v: Option<&Json>) -> Result<Vec<Merge>, Strin
 mod tests {
     use super::*;
 
+    use rustgis_geom::haversine_m;
+
     /// The ground distance from `p` to a line, by sampling each edge finely along its great circle.
+    /// Only the samples within 1.2% of the nearest on the sphere are measured on the ellipsoid,
+    /// whose distances are within 0.6% of the sphere's.
     fn brute(p: Coord, line: &[Coord]) -> f64 {
-        line.windows(2)
+        let samples: Vec<(Coord, f64)> = line
+            .windows(2)
             .flat_map(|w| (0..=400).map(move |k| geodesic_interpolate(w[0], w[1], f64::from(k) / 400.0)))
-            .map(|c| haversine_m(p, c))
-            .fold(f64::INFINITY, f64::min)
+            .map(|c| (c, haversine_m(p, c)))
+            .collect();
+        let nearest = samples.iter().map(|(_, d)| *d).fold(f64::INFINITY, f64::min);
+        samples.iter().filter(|(_, d)| *d <= nearest * 1.012 + 1.0).map(|(c, _)| wgs84_distance_m(p, *c)).fold(f64::INFINITY, f64::min)
+    }
+
+    #[test]
+    fn distances_are_on_the_wgs84_ellipsoid() {
+        // 0.0904369° north of the equator is 10,000 m on WGS 84 (GeographicLib) but 10,056 m on the
+        // mean sphere: within 10 km in the measure, and not ruled out before it (#888).
+        let (a, b) = (Coord::new(0.0, 0.0), Coord::new(0.0, 0.090_436_9));
+        assert!(haversine_m(a, b) > 10_050.0);
+        let d = distance_m(&Geometry::Point(a), &Geometry::Point(b)).unwrap();
+        assert!((d - 9_999.994_808).abs() < 1e-5, "{d}");
+        let (ea, eb) = (Extent::new(a.x, a.y, a.x, a.y), Extent::new(b.x, b.y, b.x, b.y));
+        assert!(may_be_within(&ea, &eb, 10_000.0));
+        assert!(search_box(&ea, 10_000.0).is_some_and(|e| e.max_y >= b.y));
     }
 
     #[test]
@@ -370,7 +391,7 @@ mod tests {
         // between the two real vertices, past the point between them.
         let line = Geometry::Line(vec![vec![Coord::new(0.0, 0.0), Coord::new(600.0, 80.0), Coord::new(10.0, 0.0)]]);
         let d = distance_m(&Geometry::Point(Coord::new(5.0, 0.0)), &line).unwrap();
-        let truth = haversine_m(Coord::new(5.0, 0.0), Coord::new(0.0, 0.0));
+        let truth = wgs84_distance_m(Coord::new(5.0, 0.0), Coord::new(0.0, 0.0));
         assert!((d - truth).abs() < 1e-6 * truth, "{d} vs {truth}");
         // Nor does a ring close across the gap.
         let ring = Geometry::Polygon(vec![vec![vec![
@@ -393,7 +414,7 @@ mod tests {
         let a = Geometry::MultiPoint(vec![Coord::new(f64::MAX, 0.0), origin]);
         for pts in [vec![Coord::new(-f64::MAX, 0.0), near], vec![near, Coord::new(-f64::MAX, 0.0)]] {
             let d = distance_m(&a, &Geometry::MultiPoint(pts.clone()));
-            let truth = haversine_m(origin, near);
+            let truth = wgs84_distance_m(origin, near);
             assert!(d.is_some_and(|d| (d - truth).abs() < 1e-6 * truth), "{pts:?}: {d:?} vs {truth}");
         }
     }

crates/analysis/src/tissot.rs

@@ -208,7 +208,8 @@ pub fn run(inp: &dyn Inputs, p: &Json) -> Result<ToolResult, String> {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use rustgis_geom::{haversine_m, lonlat_to_merc};
+    use geo::{Distance, Geodesic, Point};
+    use rustgis_geom::lonlat_to_merc;
 
     /// No map: no projection, no view.
     struct NoMap;
@@ -271,7 +272,9 @@ mod tests {
             assert_eq!(ring.len(), VERTICES + 1);
             assert_eq!(ring.first(), ring.last());
             for v in ring {
-                assert!((haversine_m(centre, *v) - 500_000.0).abs() < 1.0, "{centre:?} → {v:?}");
+                // 500 km on the WGS 84 ellipsoid, as GeographicLib measures it through `geo`.
+                let d = Geodesic.distance(Point::new(centre.x, centre.y), Point::new(v.x, v.y));
+                assert!((d - 500_000.0).abs() < 1e-3, "{centre:?} → {v:?}");
             }
         }
     }

crates/analysis/src/tools.rs

@@ -2124,8 +2124,8 @@ mod tests {
         let Geometry::Polygon(polys) = ds.features[0].geometry.as_ref().unwrap() else { panic!("expected buffer polygon") };
         // 500 international feet is exactly 152.4 metres, at every vertex.
         for vertex in polys.iter().flatten().flatten() {
-            let radius = rustgis_geom::haversine_m(center, *vertex);
-            assert!((radius - 152.4).abs() < 0.05, "radius: {radius}");
+            let radius = on_wgs84_m(center, *vertex);
+            assert!((radius - 152.4).abs() < 1e-6, "radius: {radius}");
         }
     }
 
@@ -2349,15 +2349,31 @@ mod tests {
         }
     }
 
-    /// The great-circle distance (km) from `c` to a line, sampled every 500 m along it.
-    fn km_to(c: Coord, line: &[Coord]) -> f64 {
-        line.windows(2)
+    /// The geodesic distance on the WGS 84 ellipsoid (metres), from GeographicLib through `geo`.
+    fn on_wgs84_m(a: Coord, b: Coord) -> f64 {
+        use geo::{Distance, Geodesic, Point};
+        Geodesic.distance(Point::new(a.x, a.y), Point::new(b.x, b.y))
+    }
+
+    /// The distance on WGS 84 (km) from `c` to a line, sampled every `step` metres along it.
+    /// Only the samples within 1.2% of the nearest on the sphere are measured on the ellipsoid,
+    /// whose distances are within 0.6% of the sphere's.
+    fn km_to_every(c: Coord, line: &[Coord], step: f64) -> f64 {
+        let samples: Vec<(Coord, f64)> = line
+            .windows(2)
             .flat_map(|w| {
-                let n = (rustgis_geom::haversine_m(w[0], w[1]) / 500.0).ceil().max(1.0) as usize;
+                let n = (rustgis_geom::haversine_m(w[0], w[1]) / step).ceil().max(1.0) as usize;
                 (0..=n).map(move |i| rustgis_geom::geodesic_interpolate(w[0], w[1], i as f64 / n as f64))
             })
-            .map(|p| rustgis_geom::haversine_m(c, p) / 1000.0)
-            .fold(f64::INFINITY, f64::min)
+            .map(|p| (p, rustgis_geom::haversine_m(c, p)))
+            .collect();
+        let nearest = samples.iter().map(|(_, d)| *d).fold(f64::INFINITY, f64::min);
+        samples.iter().filter(|(_, d)| *d <= nearest * 1.012 + 1.0).map(|(p, _)| on_wgs84_m(c, *p) / 1000.0).fold(f64::INFINITY, f64::min)
+    }
+
+    /// The distance on WGS 84 (km) from `c` to a line, sampled every 500 m along it.
+    fn km_to(c: Coord, line: &[Coord]) -> f64 {
+        km_to_every(c, line, 500.0)
     }
 
     /// Every vertex of every ring is `km` from `line`, within 0.2%. Where the arcs of two pieces
@@ -2371,6 +2387,28 @@ mod tests {
         }
     }
 
+    /// Buffers are geodesic on the WGS 84 ellipsoid, not on a sphere (#888): on a sphere of the
+    /// mean radius a 10 km buffer reached 9,944 to 10,011 m at the equator and 9,997 to 10,029 m
+    /// at 46.5°N, and a 1 km buffer of a 2 km line 999.7 to 1,002.9 m.
+    #[test]
+    fn buffers_keep_their_distance_on_the_ellipsoid() {
+        for centre in [Coord::new(0.0, 0.0), Coord::new(10.5, 46.5), Coord::new(30.0, 60.0), Coord::new(-70.0, -33.4)] {
+            for c in buffered(Geometry::Point(centre), 10.0).iter().flatten().flatten() {
+                let d = on_wgs84_m(centre, *c);
+                assert!((d - 10_000.0).abs() < 1e-3, "{centre:?}: {c:?} is {d} m away");
+            }
+        }
+        let a = Coord::new(10.5, 46.5);
+        for azimuth in [90.0, 0.0, 45.0] {
+            let line = [a, rustgis_geom::geodesic_destination(a, azimuth, 2_000.0)];
+            for c in buffered(Geometry::Line(vec![line.to_vec()]), 1.0).iter().flatten().flatten() {
+                // Sampled every 0.5 m along the line, the nearest sample is at most 0.03 mm further.
+                let d = km_to_every(*c, &line, 0.5) * 1000.0;
+                assert!((d - 1_000.0).abs() < 1e-3, "{azimuth}: {c:?} is {d} m away");
+            }
+        }
+    }
+
     /// Buffers are geodesic (#394): Web Mercator's scale changes with latitude, so one scale for
     /// the whole feature made a 1000 km circle at 60°N reach 872 km north and 1142 km south, and a
     /// 50 km buffer of a line from 5°S to 33°S 59 km wide at one end and 50 km at the other.
@@ -2381,7 +2419,9 @@ mod tests {
         all_at(&circle, &[centre, centre], 1000.0);
         let north = circle.iter().flatten().flatten().map(|c| c.y).fold(f64::MIN, f64::max);
         let south = circle.iter().flatten().flatten().map(|c| c.y).fold(f64::MAX, f64::min);
-        assert!((north - 60.0 - (60.0 - south)).abs() < 0.01, "as far north as south: {north} {south}");
+        // As far north as south, where GeographicLib puts 1000 km on WGS 84 (a degree of
+        // latitude grows towards the pole, so 8.970° north and 8.982° south).
+        assert!((north - 68.969_947_795_748).abs() < 1e-9 && (south - 51.017_861_794_779).abs() < 1e-9, "{north} {south}");
 
         let line = [Coord::new(-45.0, -5.0), Coord::new(-45.0, -33.0)];
         all_at(&buffered(Geometry::Line(vec![line.to_vec()]), 50.0), &line, 50.0);
@@ -2972,7 +3012,8 @@ mod tests {
         ]));
         let run = |params: serde_json::Value| layer(crate::run("spatialJoin", &m, &params).unwrap_or_else(|e| panic!("{e}")));
         let num_at = |ds: &Dataset, row: usize, f: &str| ds.value(row, ds.field_index(f).unwrap()).as_f64();
-        let hav = |a: (f64, f64), b: (f64, f64)| rustgis_geom::haversine_m(Coord::new(a.0, a.1), Coord::new(b.0, b.1));
+        // Distances on WGS 84, as the join measures them (#888).
+        let hav = |a: (f64, f64), b: (f64, f64)| rustgis_geom::wgs84_distance_m(Coord::new(a.0, a.1), Coord::new(b.0, b.1));
 
         // Closest, unlimited: the nearest stop, across the antimeridian and over the pole too.
         let closest =
@@ -2996,7 +3037,7 @@ mod tests {
         );
         assert!((num_at(&feet, 0, "ft").unwrap() - s0 / 0.3048).abs() < 1e-3);
         // A line's nearest point can be between its vertices: north of the home, 2° of longitude
-        // from the vertices. The edge is a great circle, a little poleward of latitude 1 midway.
+        // from the vertices. The edge is a geodesic, a little poleward of latitude 1 midway.
         let to_road =
             run(serde_json::json!({"target": "homes", "join": "road", "match": "Closest", "units": "Meters", "distanceField": "d"}));
         let (north, bulge) = (hav((0.004, 0.0), (0.004, 1.0)), hav((0.004, 0.0), (0.004, 1.000_153)));
@@ -3046,7 +3087,7 @@ mod tests {
         }
         let brute = |x: f64, y: f64| {
             let d = |f: &rustgis_data::Feature| match &f.geometry {
-                Some(Geometry::Point(c)) => rustgis_geom::haversine_m(Coord::new(x, y), *c),
+                Some(Geometry::Point(c)) => rustgis_geom::wgs84_distance_m(Coord::new(x, y), *c),
                 _ => f64::INFINITY,
             };
             let best = many.features.iter().min_by(|a, b| d(a).total_cmp(&d(b))).map(|f| f.attrs[0].clone());
@@ -3101,6 +3142,37 @@ mod tests {
         assert!(crate::run("spatialJoin", &m, &serde_json::json!({"target": "homes", "join": "stops", "distance": " "})).is_ok());
     }
 
+    /// Within a distance through the join's grid index (more than 16 join features), at the edge
+    /// of 10 km due north of the equator (#888): 0.0904369° is 9,999.995 m on WGS 84 but 10,056 m
+    /// on the mean sphere, so a search worked out on that sphere left it out; 0.0905° is
+    /// 10,006.972 m and stays out.
+    #[test]
+    fn spatial_join_within_a_distance_at_the_ellipsoids_edge() {
+        let point = |x: f64, y: f64| Some(Geometry::Point(Coord::new(x, y)));
+        let mut stops = Dataset::new("stops", GeometryType::Point, vec![Field::new("stop", FieldType::Text)]);
+        stops.push(point(0.0, 0.090_436_9), vec![Value::Text("edge".into())]);
+        stops.push(point(0.0, 0.0905), vec![Value::Text("past".into())]);
+        // More stops, all over 124 km away, so the join indexes them.
+        for k in 0..18 {
+            stops.push(point(1.0 + f64::from(k) * 0.1, 0.5 + f64::from(k) * 0.02), vec![Value::Text(format!("far {k}"))]);
+        }
+        let mut homes = Dataset::new("homes", GeometryType::Point, vec![]);
+        homes.push(point(0.0, 0.0), vec![]);
+        let m = Map(HashMap::from([("stops".to_string(), stops), ("homes".to_string(), homes)]));
+        let within = layer(
+            crate::run(
+                "spatialJoin",
+                &m,
+                &serde_json::json!({"target": "homes", "join": "stops", "match": "Within a distance", "distance": 10_000, "units": "Meters", "operation": "Join one to many", "distanceField": "d"}),
+            )
+            .unwrap(),
+        );
+        let names: Vec<_> = (0..within.features.len()).filter_map(|r| text_at(&within, r, "join_stop")).collect();
+        assert_eq!(names, ["edge"]);
+        let d = within.value(0, within.field_index("d").unwrap()).as_f64().unwrap();
+        assert!((d - 9_999.994_808).abs() < 1e-5, "{d}");
+    }
+
     #[test]
     fn spatial_join_merge_rules() {
         let mut districts = Dataset::new("districts", GeometryType::Polygon, vec![Field::new("d", FieldType::Text)]);

crates/geom/Cargo.toml

@@ -9,4 +9,5 @@ description = "RustGIS L0: coordinates, extents, projections, geodesic measureme
 workspace = true
 
 [dependencies]
+geographiclib-rs = { workspace = true }
 serde = { workspace = true }