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..0b5e0da44 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()); @@ -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)]); 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..b0646dbf4 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,58 @@ 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. + /// 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 { - 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()) + 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) } } @@ -377,28 +382,66 @@ 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); + // 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.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 +476,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..0ec792cca 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -111,10 +111,14 @@ 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. +Lengths and areas (Calculate Geometry Attributes, and Spatial Join's **Largest overlap**) are +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 @@ -150,10 +154,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