From cdc1aaf8ae65fff4890f291c55fdf0ed388a3696 Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Sun, 25 Oct 2026 14:41:51 -0300 Subject: [PATCH 1/3] Analysis: buffers, Tissot's circles and Spatial Join distances on the WGS 84 ellipsoid (#888) Aeqd becomes Karney's ellipsoidal azimuthal equidistant projection: forward is the inverse geodesic problem from the centre, inverse the direct one, both solved by geographiclib-rs (already in the tree through geo). A point's buffer and geodesic_destination (Tissot's circles) are exact on WGS 84; lines and polygons keep #429's planes, now on the ellipsoid. Spatial Join measures the nearest pair with wgs84_distance_m, and its pre-filter uses WGS 84's smallest radius of curvature (b^2/a), so the ellipsoid's shorter distances near the equator are never ruled out. Tests compare with GeographicLib's values (pyproj.Geod(ellps="WGS84")). Lengths and areas stay spherical (#894). Co-Authored-By: Claude Opus 5.5 --- Cargo.lock | 1 + Cargo.toml | 2 + crates/analysis/src/join.rs | 45 ++++++++---- crates/analysis/src/tissot.rs | 7 +- crates/analysis/src/tools.rs | 65 +++++++++++++---- crates/geom/Cargo.toml | 1 + crates/geom/src/lib.rs | 132 +++++++++++++++++++++------------- docs/user-guide/analysis.md | 14 ++-- 8 files changed, 187 insertions(+), 80 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index f630c165f..5d768fe07 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -4784,6 +4784,7 @@ dependencies = [ name = "rustgis-geom" version = "0.4.0" dependencies = [ + "geographiclib-rs", "serde", ] diff --git a/Cargo.toml b/Cargo.toml index 0097c4f9b..bf853d14f 100644 --- a/Cargo.toml +++ b/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"] } diff --git a/crates/analysis/src/join.rs b/crates/analysis/src/join.rs index 1128d739d..14bc2eac6 100644 --- a/crates/analysis/src/join.rs +++ b/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 { 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) { // 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, 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}"); } } diff --git a/crates/analysis/src/tissot.rs b/crates/analysis/src/tissot.rs index 5c92dfb57..eefa97874 100644 --- a/crates/analysis/src/tissot.rs +++ b/crates/analysis/src/tissot.rs @@ -208,7 +208,8 @@ pub fn run(inp: &dyn Inputs, p: &Json) -> Result { #[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:?}"); } } } diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index a13b398cf..f6ea18e03 100644 --- a/crates/analysis/src/tools.rs +++ b/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, ¶ms).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()); diff --git a/crates/geom/Cargo.toml b/crates/geom/Cargo.toml index 0ac472405..db9fb25c6 100644 --- a/crates/geom/Cargo.toml +++ b/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 } diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index da958ee76..acb5bbe5e 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -11,10 +11,17 @@ pub mod triangulate; pub use body::{BODIES, Body}; pub use sinusoidal::Sinusoidal; +use std::sync::LazyLock; + +use geographiclib_rs::{DirectGeodesic, Geodesic, InverseGeodesic}; use serde::{Deserialize, Serialize}; /// Mean Earth radius used by the haversine formula (metres). pub const EARTH_RADIUS_M: f64 = 6_371_008.8; +/// The smallest radius of curvature of the WGS 84 ellipsoid, b²/a, along the meridian at the +/// equator (metres). Every path on the ellipsoid is at least as long as on a sphere of this radius +/// through the same longitudes and latitudes, so a reach worked out on that sphere is never short. +pub const WGS84_MIN_RADIUS_M: f64 = 6_335_439.327_293; /// Web Mercator sphere radius (metres). pub const MERCATOR_R: f64 = 6_378_137.0; /// Latitude limit of Web Mercator. @@ -226,60 +233,54 @@ pub fn spherical_centre(pts: &[Coord]) -> Option { from_vector(s) } -/// The point `dist_m` metres from `c` along the great circle leaving it at `azimuth_deg` -/// (clockwise from north): the direct geodesic problem on the sphere. +/// The point `dist_m` metres from `c` along the geodesic leaving it at `azimuth_deg` (clockwise +/// from north): the direct geodesic problem on the WGS 84 ellipsoid. pub fn geodesic_destination(c: Coord, azimuth_deg: f64, dist_m: f64) -> Coord { let (s, co) = azimuth_deg.to_radians().sin_cos(); Aeqd::new(c).inverse(Coord::new(dist_m * s, dist_m * co)) } -/// The azimuthal equidistant projection on the sphere, centred on a lon/lat point: distances -/// and directions from the centre are true, so a buffer drawn in it is geodesic near the centre. -/// Coordinates are metres east and north of the centre. +/// The geodesic distance between two lon/lat points on the WGS 84 ellipsoid (metres): the inverse +/// geodesic problem, solved by GeographicLib. +pub fn wgs84_distance_m(a: Coord, b: Coord) -> f64 { + let (s12, _): (f64, f64) = WGS84.inverse(a.y.clamp(-90.0, 90.0), a.x, b.y.clamp(-90.0, 90.0), b.x); + s12 +} + +/// The WGS 84 ellipsoid's geodesic solver. +static WGS84: LazyLock = LazyLock::new(Geodesic::wgs84); + +/// The azimuthal equidistant projection on the WGS 84 ellipsoid, centred on a lon/lat point: +/// distances and directions from the centre are true, so a buffer drawn in it is geodesic near +/// the centre. Coordinates are metres east and north of the centre. +/// +/// It is Karney's ellipsoidal form (Karney 2013, "Algorithms for geodesics"): a point's +/// distance and azimuth from the centre are the inverse geodesic problem, and back again the +/// direct one, both solved by GeographicLib to a few nanometres. #[derive(Clone, Copy, Debug)] pub struct Aeqd { lon0: f64, - sin0: f64, - cos0: f64, + lat0: f64, } impl Aeqd { pub fn new(centre: Coord) -> Self { - let lat0 = centre.y.clamp(-90.0, 90.0).to_radians(); - Aeqd { lon0: centre.x, sin0: lat0.sin(), cos0: lat0.cos() } + Aeqd { lon0: centre.x, lat0: centre.y.clamp(-90.0, 90.0) } } /// Lon/lat (degrees) → metres from the centre. pub fn forward(&self, c: Coord) -> Coord { - let lat = c.y.clamp(-90.0, 90.0).to_radians(); - let dlon = (c.x - self.lon0).to_radians(); - let (sl, cl) = lat.sin_cos(); - let east = cl * dlon.sin(); - let north = self.cos0 * sl - self.sin0 * cl * dlon.cos(); - let cos_c = self.sin0 * sl + self.cos0 * cl * dlon.cos(); - let sin_c = east.hypot(north); - if sin_c < 1e-15 { - return Coord::new(0.0, 0.0); - } - let k = EARTH_RADIUS_M * sin_c.atan2(cos_c) / sin_c; - Coord::new(k * east, k * north) + let (s12, azi1, _, _): (f64, f64, f64, f64) = WGS84.inverse(self.lat0, self.lon0, c.y.clamp(-90.0, 90.0), c.x); + let (sa, ca) = azi1.to_radians().sin_cos(); + Coord::new(s12 * sa, s12 * ca) } /// Metres from the centre → lon/lat (degrees). The longitude is within 180° of the centre's, /// not folded into −180..180, so a shape across the antimeridian stays in one piece. pub fn inverse(&self, p: Coord) -> Coord { - let rho = p.x.hypot(p.y); - let c = rho / EARTH_RADIUS_M; - let (sc, cc) = c.sin_cos(); - let (sa, ca) = if rho > 0.0 { (p.x / rho, p.y / rho) } else { (0.0, 1.0) }; - let lat = (self.sin0 * cc + self.cos0 * sc * ca).clamp(-1.0, 1.0).asin(); - // At a pole every direction is south (or north): the azimuth is the meridian. - let dlon = if self.cos0 < 1e-9 { - if self.sin0 > 0.0 { p.x.atan2(-p.y) } else { p.x.atan2(p.y) } - } else { - (sa * sc * self.cos0).atan2(cc - self.sin0 * lat.sin()) - }; - Coord::new(self.lon0 + dlon.to_degrees(), lat.to_degrees()) + // At a pole the azimuth is taken from the meridian `lon0`, as `forward` gives it. + let (lat, lon): (f64, f64) = WGS84.direct(self.lat0, self.lon0, p.x.atan2(p.y).to_degrees(), p.x.hypot(p.y)); + Coord::new(self.lon0 + wrap_lon(lon - self.lon0), lat) } } @@ -377,28 +378,61 @@ mod tests { #[test] fn azimuthal_equidistant_keeps_distances_from_its_centre() { - let centre = Coord::new(10.0, 60.0); - let p = Aeqd::new(centre); - for (x, y) in [(10.0, 70.0), (-170.0, 80.0), (40.0, 20.0), (10.0, 60.0), (179.5, -30.0)] { - let c = Coord::new(x, y); + // Distances and azimuths on WGS 84 from GeographicLib (pyproj.Geod(ellps="WGS84").inv). + for (centre, (x, y), s12, azi1) in [ + ((10.0, 60.0), (10.0, 70.0), 1_114_907.908_280, 0.0), + ((10.0, 60.0), (-170.0, 80.0), 4_464_718.767_198, 0.0), + ((10.0, 60.0), (40.0, 20.0), 5_037_117.406_299, 138.557_618_275_686), + ((10.0, 60.0), (179.5, -30.0), 16_578_966.423_682, 17.827_025_777_178), + ((2.3522, 48.8566), (-0.1276, 51.5072), 343_896.891_267, -30.048_593_801_903), + ] { + let (centre, c) = (Coord::new(centre.0, centre.1), Coord::new(x, y)); + let p = Aeqd::new(centre); let m = p.forward(c); - assert!((m.x.hypot(m.y) - haversine_m(centre, c)).abs() < 1e-3, "{c:?}"); + assert!((m.x.hypot(m.y) - s12).abs() < 1e-5, "{c:?}: {m:?}"); + assert!((m.x.atan2(m.y).to_degrees() - azi1).abs() < 1e-9, "{c:?}: {m:?}"); let back = p.inverse(m); assert!((back.x - x).abs() < 1e-9 && (back.y - y).abs() < 1e-9, "{c:?} → {back:?}"); } + assert_eq!(Aeqd::new(Coord::new(10.0, 60.0)).forward(Coord::new(10.0, 60.0)), Coord::new(0.0, 0.0)); + assert!((wgs84_distance_m(Coord::new(2.3522, 48.8566), Coord::new(-0.1276, 51.5072)) - 343_896.891_267).abs() < 1e-5); + // A tenth of a degree north from the equator is shorter than on the mean sphere, but not + // than on the sphere of the smallest radius of curvature. + let (a, b) = (Coord::new(0.0, 0.0), Coord::new(0.0, 0.1)); + assert!(wgs84_distance_m(a, b) < haversine_m(a, b)); + assert!(wgs84_distance_m(a, b) >= haversine_m(a, b) * WGS84_MIN_RADIUS_M / EARTH_RADIUS_M); // The inverse stays next to the centre's longitude rather than folding into −180..180. - assert!((Aeqd::new(Coord::new(179.0, 0.0)).inverse(Coord::new(300_000.0, 0.0)).x - 181.698).abs() < 1e-3); + assert!((Aeqd::new(Coord::new(179.0, 0.0)).inverse(Coord::new(300_000.0, 0.0)).x - 181.694_945_852_359).abs() < 1e-9); + // Centred on a pole, the azimuth is taken from the centre's meridian, both ways. + for pole in [90.0, -90.0] { + let p = Aeqd::new(Coord::new(25.0, pole)); + for c in [Coord::new(25.0, pole * 0.9), Coord::new(-100.0, pole * 0.8), Coord::new(170.0, pole * 0.95)] { + let back = p.inverse(p.forward(c)); + assert!((wrap_lon(back.x - c.x)).abs() < 1e-9 && (back.y - c.y).abs() < 1e-9, "{c:?} → {back:?}"); + } + } } #[test] fn geodesic_destination_and_interpolation() { - let start = Coord::new(10.0, 60.0); - for az in [0.0, 45.0, 90.0, 180.0, 300.0] { - let d = geodesic_destination(start, az, 1_000_000.0); - assert!((haversine_m(start, d) - 1_000_000.0).abs() < 1e-3, "{az}: {d:?}"); + // Destinations on WGS 84 from GeographicLib (pyproj.Geod(ellps="WGS84").fwd): 10 km from + // the equator and from 46.5° and 60°N (#888), 1,000 km from 60°N, and from either pole. + for (lon, lat, az, d, lon2, lat2) in [ + (0.0, 0.0, 0.0, 10_000.0, 0.0, 0.090_436_946_951), + (0.0, 0.0, 90.0, 10_000.0, 0.089_831_528_412, 0.0), + (10.5, 46.5, 0.0, 10_000.0, 10.5, 46.589_958_829_711), + (10.5, 46.5, 90.0, 10_000.0, 10.630_271_666_436, 46.499_925_816_591), + (10.5, 46.5, 225.0, 10_000.0, 10.407_991_242_765, 46.436_351_607_689), + (30.0, 60.0, 90.0, 10_000.0, 30.179_211_026_165, 59.999_878_434_646), + (10.0, 60.0, 0.0, 1_000_000.0, 10.0, 68.969_947_795_748), + (10.0, 60.0, 300.0, 1_000_000.0, -7.552_217_029_942, 63.435_692_055_327), + (0.0, 90.0, 0.0, 1_000_000.0, 180.0, 81.046_232_815_951), + (0.0, -90.0, 90.0, 1_000_000.0, 90.0, -81.046_232_815_951), + ] { + let p = geodesic_destination(Coord::new(lon, lat), az, d); + // 1e-11° is about a micrometre. + assert!(wrap_lon(p.x - lon2).abs() < 1e-11 && (p.y - lat2).abs() < 1e-11, "{lon} {lat} {az}: {p:?}"); } - let north = geodesic_destination(start, 0.0, 1_000_000.0); - assert!((north.x - 10.0).abs() < 1e-9 && north.y > 68.9, "{north:?}"); let (a, b) = (Coord::new(2.3522, 48.8566), Coord::new(-0.1276, 51.5072)); let mid = geodesic_interpolate(a, b, 0.5); assert!((haversine_m(a, mid) - haversine_m(mid, b)).abs() < 1e-3); @@ -433,10 +467,12 @@ mod tests { let from = Coord::new(10.0, 60.0); for az in [0.0, 45.0, 90.0, 180.0, 270.0, 333.0] { let to = geodesic_destination(from, az, 500_000.0); - assert!((haversine_m(from, to) - 500_000.0).abs() < 1e-3, "{az}: {to:?}"); + let m = Aeqd::new(from).forward(to); + assert!((m.x.hypot(m.y) - 500_000.0).abs() < 1e-6, "{az}: {to:?}"); + assert!(wrap_lon(m.x.atan2(m.y).to_degrees() - az).abs() < 1e-9, "{az}: {to:?}"); } - // North and east go where they say. - assert!((geodesic_destination(from, 0.0, 111_195.0).y - 61.0).abs() < 1e-3); + // North and east go where they say: 60° to 61°N is 111,420.728 m on WGS 84 (GeographicLib). + assert!((geodesic_destination(from, 0.0, 111_420.727_870).y - 61.0).abs() < 1e-9); assert!(geodesic_destination(from, 90.0, 1000.0).x > 10.0); // Longitudes stay continuous across the antimeridian, so the points around one centre // form one shape (Tissot's indicatrix draws its circles with them). diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 3364ebf2b..9d2fc996e 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -111,8 +111,9 @@ Calculate Geometry Attributes updates **all** rows. A built-in tool's `output` i the project, never a file path. **Buffer** distances, and the distance of **Select Layer By Location**, are measured on the Earth -(a sphere of radius 6,371 km), not in the map's display projection. Each feature is buffered in -an azimuthal equidistant projection centred on it; long lines and large polygon outlines are +(the WGS 84 ellipsoid, as PROJ and GeographicLib measure it), not in the map's display +projection. Each feature is buffered in an azimuthal equidistant projection on the ellipsoid +centred on it, so a point's buffer is exact; long lines and large polygon outlines are buffered in stretches of a few hundred kilometres to limit projection error. Negative Buffer distances shrink polygons. Buffer and Select Layer By Location limit distances to 5,000 km. @@ -150,10 +151,11 @@ stored in longitude/latitude before the join. `match` is one of: one, every target gets its nearest feature. Ties keep the join layer's order. - **Within a distance**: every join feature within `distance` (required), nearest first. -Distances are on the ground, between the features' nearest points (edges and vertices alike), -and are zero for features that intersect. `distanceField` names a field added for the distance -to the joined feature, in `units`; it applies to **Closest** and **Within a distance** only, as -`distance` does, and the Tools pane shows them only for those. +Distances are on the ground (the WGS 84 ellipsoid, as for Buffer), between the features' nearest +points (edges and vertices alike), and are zero for features that intersect. `distanceField` +names a field added for the distance to the joined feature, in `units`; it applies to +**Closest** and **Within a distance** only, as `distance` does, and the Tools pane shows them +only for those. **Join one to one** (the default `operation`) copies the first matching feature's attributes; **Join one to many** makes a row for each match. Both add `join_count`, the total matches for From 620cc337ca03b72ccd0e87d95eefb41e65956830 Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Sun, 25 Oct 2026 16:18:05 -0300 Subject: [PATCH 2/3] Spatial Join: a Within a distance test through the grid index at the ellipsoid's 10 km edge; docs note lengths and areas stay spherical (#888) Review follow-up for #905: 20 stops so the join indexes them; the stop 0.0904369 deg north of the equator (9,999.995 m on WGS 84, 10,056 m on the mean sphere) is joined within 10,000 m and the one at 0.0905 deg is not. The test fails with the mean radius back in the join's pre-filter. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/tools.rs | 31 +++++++++++++++++++++++++++++++ docs/user-guide/analysis.md | 3 +++ 2 files changed, 34 insertions(+) diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index f6ea18e03..0b5e0da44 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -3142,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)]); diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 9d2fc996e..a1d3527fd 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -116,6 +116,9 @@ projection. Each feature is buffered in an azimuthal equidistant projection on t centred on it, so a point's buffer is exact; long lines and large polygon outlines are buffered in stretches of a few hundred kilometres to limit projection error. Negative Buffer distances shrink polygons. Buffer and Select Layer By Location limit distances to 5,000 km. +Lengths and areas (Calculate Geometry Attributes, and Spatial Join's **Largest overlap**) are +still measured on a sphere of radius 6,371 km, which differs from the ellipsoid by up to about +0.5 %. **IDW Interpolation** works like the QGIS tool of the same name: every point weighs on every cell by 1/distanceᴾ, where P is the distance coefficient (2 by default). Lines and polygons From 91e58efcca0537e193793ca59f3eafe30bd55d99 Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Sun, 25 Oct 2026 16:42:23 -0300 Subject: [PATCH 3/3] Aeqd: metres that aren't finite give no place; docs give separate length and area bounds (#888) Review follow-up for #905: Aeqd::inverse returned a finite longitude with a NaN latitude for NaN or infinite metres; it now returns NaN for both. The user guide's note on spherical lengths and areas now says up to about 0.6 % for lengths and 0.9 % for areas (a small area near a pole), not 0.5 %. Co-Authored-By: Claude Opus 5.5 --- crates/geom/src/lib.rs | 11 ++++++++++- docs/user-guide/analysis.md | 4 ++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index acb5bbe5e..b0646dbf4 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -276,8 +276,12 @@ impl Aeqd { } /// Metres from the centre → lon/lat (degrees). The longitude is within 180° of the centre's, - /// not folded into −180..180, so a shape across the antimeridian stays in one piece. + /// not folded into −180..180, so a shape across the antimeridian stays in one piece. Metres + /// that aren't finite give no place (both coordinates NaN). pub fn inverse(&self, p: Coord) -> Coord { + if !(p.x.is_finite() && p.y.is_finite()) { + return Coord::new(f64::NAN, f64::NAN); + } // At a pole the azimuth is taken from the meridian `lon0`, as `forward` gives it. let (lat, lon): (f64, f64) = WGS84.direct(self.lat0, self.lon0, p.x.atan2(p.y).to_degrees(), p.x.hypot(p.y)); Coord::new(self.lon0 + wrap_lon(lon - self.lon0), lat) @@ -401,6 +405,11 @@ mod tests { let (a, b) = (Coord::new(0.0, 0.0), Coord::new(0.0, 0.1)); assert!(wgs84_distance_m(a, b) < haversine_m(a, b)); assert!(wgs84_distance_m(a, b) >= haversine_m(a, b) * WGS84_MIN_RADIUS_M / EARTH_RADIUS_M); + // Metres that aren't finite are no place, not half of one. + for m in [Coord::new(f64::NAN, 0.0), Coord::new(f64::INFINITY, 0.0), Coord::new(0.0, f64::NEG_INFINITY)] { + let c = Aeqd::new(Coord::new(10.0, 45.0)).inverse(m); + assert!(c.x.is_nan() && c.y.is_nan(), "{m:?} → {c:?}"); + } // The inverse stays next to the centre's longitude rather than folding into −180..180. assert!((Aeqd::new(Coord::new(179.0, 0.0)).inverse(Coord::new(300_000.0, 0.0)).x - 181.694_945_852_359).abs() < 1e-9); // Centred on a pole, the azimuth is taken from the centre's meridian, both ways. diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index a1d3527fd..0ec792cca 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -117,8 +117,8 @@ centred on it, so a point's buffer is exact; long lines and large polygon outlin buffered in stretches of a few hundred kilometres to limit projection error. Negative Buffer distances shrink polygons. Buffer and Select Layer By Location limit distances to 5,000 km. Lengths and areas (Calculate Geometry Attributes, and Spatial Join's **Largest overlap**) are -still measured on a sphere of radius 6,371 km, which differs from the ellipsoid by up to about -0.5 %. +still measured on a sphere of radius 6,371 km: on the ellipsoid, lengths differ by up to about +0.6 % and areas by up to about 0.9 % (a small area near a pole). **IDW Interpolation** works like the QGIS tool of the same name: every point weighs on every cell by 1/distanceᴾ, where P is the distance coefficient (2 by default). Lines and polygons