diff --git a/crates/analysis/src/convert.rs b/crates/analysis/src/convert.rs index ef7a2e7c9..b71885103 100644 --- a/crates/analysis/src/convert.rs +++ b/crates/analysis/src/convert.rs @@ -1,18 +1,26 @@ -//! Conversions between RustGIS geometry and `geo` types (in Web Mercator metres). +//! Conversions between RustGIS geometry and `geo` types, in the analysis plane +//! ([`rustgis_geom::lonlat_to_plane`]): Web Mercator meters within ±85.05°, carried on to the poles +//! beyond, so polar geometry isn't clamped onto 85.05° (#1033). use geo::{Coord as GC, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon}; use rustgis_data::Geometry; -use rustgis_geom::{Coord, lonlat_to_merc, merc_to_lonlat}; +use rustgis_geom::{Coord, lonlat_to_plane, plane_to_lonlat}; /// Projects without folding the longitude: a ring kept unwrapped across the antimeridian (as a /// geodesic buffer is) stays one piece, rather than folding its far vertex across the world. fn gc(c: Coord) -> GC { - let m = lonlat_to_merc(c); + let m = lonlat_to_plane(c); GC { x: m.x, y: m.y } } fn back(c: &GC) -> Coord { - merc_to_lonlat(Coord::new(c.x, c.y)) + plane_to_lonlat(Coord::new(c.x, c.y)) +} + +/// The poles something spanning latitudes `south` to `north` reaches (to within +/// [`rustgis_geom::POLE_TOLERANCE`]): the North Pole, then the South Pole. +pub(crate) fn reaches_poles(south: f64, north: f64) -> (bool, bool) { + (north >= 90.0 - rustgis_geom::POLE_TOLERANCE, south <= -90.0 + rustgis_geom::POLE_TOLERANCE) } pub fn ring(r: &[Coord]) -> LineString { @@ -96,7 +104,7 @@ mod tests { let xs: Vec = mp.0.iter().flat_map(|p| p.exterior().0.iter().map(|c| c.x)).collect(); let width = xs.iter().copied().fold(f64::NEG_INFINITY, f64::max) - xs.iter().copied().fold(f64::INFINITY, f64::min); assert!(width > 20_000.0 && width < 25_000.0, "{width}"); - let back: Vec = mp.0.iter().flat_map(|p| p.exterior().0.iter().map(|c| merc_to_lonlat(Coord::new(c.x, c.y)))).collect(); + let back: Vec = mp.0.iter().flat_map(|p| p.exterior().0.iter().map(|c| plane_to_lonlat(Coord::new(c.x, c.y)))).collect(); for (a, b) in ring.iter().zip(back.iter()) { assert!((a.x - b.x).abs() < 1e-9 && (a.y - b.y).abs() < 1e-9, "{a:?} → {b:?}"); } diff --git a/crates/analysis/src/geodesic.rs b/crates/analysis/src/geodesic.rs index ff03d781c..7e1316324 100644 --- a/crates/analysis/src/geodesic.rs +++ b/crates/analysis/src/geodesic.rs @@ -8,7 +8,9 @@ //! //! The outline keeps its distance between its vertices too, not only at them (#1159): an edge //! drawn straight in the plane is split before it comes back to lon/lat wherever a straight edge -//! there would stray from it ([`OUTLINE_TOLERANCE`]), so up to 85° of latitude the sides are +//! there would stray from it ([`OUTLINE_TOLERANCE`]), straight as the tools draw edges (in the +//! analysis plane: Web Mercator within ±85.05°, straight in lon/lat beyond, #1033), so at any +//! latitude, as far as [`MAX_SPLITS`] halvings reach, the sides are //! within 0.012% of the distance (or 1 cm) wherever the feature is within 150 km of its //! projection's centre: the projection's own error, under 0.0093% there, and the edges' 0.002%. //! A winding line or outline buffered in stretches can have vertices up to [`STRETCH_M`] from a @@ -26,7 +28,7 @@ use geo::{ orient::Direction, }; use rustgis_data::{Coord, Geometry}; -use rustgis_geom::{Aeqd, geodesic_interpolate, haversine_m, lonlat_to_merc, merc_to_lonlat, spherical_centre}; +use rustgis_geom::{Aeqd, geodesic_interpolate, haversine_m, lonlat_to_plane, plane_to_lonlat, spherical_centre}; /// A feature within this distance of its centre is buffered in one projection (the error is /// under 0.01%). @@ -240,9 +242,9 @@ impl Plane { } /// The points to add between `a` and `b`, an edge straight in the plane, so that back in lon/lat - /// straight edges between them stay within `tolerance` metres of it. The edge's quarter points - /// and midpoint are compared with the same points of the straight edge in Web Mercator, the - /// plane the tools draw edges in (one point alone misses a stray that changes side, as along a + /// straight edges between them stay within `tolerance` meters of it. The edge's quarter points + /// and midpoint are compared with the same points of the straight edge in the analysis plane, + /// the plane the tools draw edges in (one point alone misses a stray that changes side, as along a /// line crossing the equator); an edge that strays is halved, and each half checked again, up /// to [`MAX_SPLITS`] times, while the ring's `budget` of added points lasts. Each point, the /// ends and those added to `out`, comes with its place in lon/lat. @@ -273,27 +275,19 @@ impl Plane { } /// Whether the edge from `a` to `b` (straight in the plane, `la` and `lb` in lon/lat) strays - /// from the same edge straight in Web Mercator by more than `tolerance` metres at its quarter - /// points or midpoint. + /// from the same edge straight in the analysis plane ([`lonlat_to_plane`]: Web Mercator within + /// ±85.05°, straight in lon/lat beyond) by more than `tolerance` meters at its quarter points + /// or midpoint. fn strays(&self, a: GC, b: GC, la: Coord, lb: Coord, tolerance: f64) -> bool { let at = |t: f64| self.0.inverse(Coord::new(a.x + t * (b.x - a.x), a.y + t * (b.y - a.y))); - let beyond = |c: Coord| c.y.abs() >= rustgis_geom::MAX_LAT; - match (beyond(la), beyond(lb)) { - // Beyond Web Mercator's 85.05° there's no straight edge to compare with. - (true, true) => false, - // Across that edge: halved until each piece is on one side of it or too short to stray. - (true, false) | (false, true) => true, - (false, false) => { - // `b` on `a`'s side of the antimeridian, so the straight edge is the short one. - let near_b = Coord::new(la.x + wrap(lb.x - la.x), lb.y); - let (ma, mb) = (lonlat_to_merc(la), lonlat_to_merc(near_b)); - [0.25, 0.5, 0.75].into_iter().any(|t| { - let p = at(t); - let chord = merc_to_lonlat(Coord::new(ma.x + t * (mb.x - ma.x), ma.y + t * (mb.y - ma.y))); - haversine_m(p, Coord::new(p.x + wrap(chord.x - p.x), chord.y)) > tolerance - }) - } - } + // `b` on `a`'s side of the antimeridian, so the straight edge is the short one. + let near_b = Coord::new(la.x + wrap(lb.x - la.x), lb.y); + let (ma, mb) = (lonlat_to_plane(la), lonlat_to_plane(near_b)); + [0.25, 0.5, 0.75].into_iter().any(|t| { + let p = at(t); + let chord = plane_to_lonlat(Coord::new(ma.x + t * (mb.x - ma.x), ma.y + t * (mb.y - ma.y))); + haversine_m(p, Coord::new(p.x + wrap(chord.x - p.x), chord.y)) > tolerance + }) } /// A ring round a pole, which goes once round the globe in longitude, as a lon/lat ring: cut at @@ -333,8 +327,11 @@ impl Plane { } /// Whether the edge from `a` to `b` in the plane (`la` and `lb` in lon/lat) can stray as far as -/// `tolerance` from the same edge straight in Web Mercator, so that [`Plane::strays`] must look. -/// Straight in Web Mercator is a rhumb line, which runs off the geodesic by about L² tan φ / 8R +/// `tolerance` from the same edge straight in the analysis plane, so that [`Plane::strays`] must +/// look. Beyond 85.05° the analysis plane's straight edges are straight in lon/lat: spirals about a +/// pole, not rhumb lines, bounded as such; an edge reaching across 85.05°, or with its ends beyond +/// it about both poles, is always looked at. Within ±85.05° straight in +/// the analysis plane is straight in Web Mercator, a rhumb line, which runs off the geodesic by about L² tan φ / 8R /// at latitude φ, the heading of a geodesic turning by L tan φ / R along it: 2.0e-8 per square /// metre at 46°, 1.9e-7 at 84°, as measured for [`STRAY_PER_M2`], plus [`FLAT_ALLOWANCE`] for the /// ellipsoid and near the equator. Straight in the plane is the @@ -345,18 +342,45 @@ impl Plane { /// the look that most of a buffer's time went to. fn may_stray(a: GC, b: GC, la: Coord, lb: Coord, tolerance: f64) -> bool { let length2 = (b.x - a.x).powi(2) + (b.y - a.y).powi(2); + let length = length2.sqrt(); + let rho = a.x.hypot(a.y).max(b.x.hypot(b.y)) + length; + let (beyond_a, beyond_b) = (la.y.abs() > rustgis_geom::MAX_LAT, lb.y.abs() > rustgis_geom::MAX_LAT); + if beyond_a || beyond_b { + // [`STRAY_PER_M2`] holds for rhumb lines, not here, where even a short edge near a pole + // strays meters: an edge with an end within ±85.05° (where the plane bends) or ends about + // both poles is always looked at. + if !(beyond_a && beyond_b) || la.y * lb.y <= 0.0 { + return true; + } + // Straight in lon/lat is, about the pole, a spiral ([`spiral_stray`]). Its chord there, + // straight about the pole, turns from the geodesic as this plane's straight edges do, at + // its own distance from its center. + let from_pole = |c: Coord| (90.0 - c.y.abs()).to_radians() * rustgis_geom::WGS84_POLAR_RADIUS_M; + let (pa, pb) = (from_pole(la), from_pole(lb)); + let spiral = spiral_stray(pa.max(pb), (pa - pb).abs(), wrap(lb.x - la.x).to_radians().abs()); + let turn = (FLAT_ALLOWANCE + PLANE_TURN * (rho + pa.max(pb)) / EARTH_R_M) / (8.0 * EARTH_R_M) * length2; + let stray = STRAY_MARGIN * (spiral + turn); + return !stray.is_finite() || stray > tolerance; + } if STRAY_PER_M2 * length2 <= tolerance { return false; } - let length = length2.sqrt(); let lat = la.y.abs().max(lb.y.abs()).to_radians() + length / EARTH_R_M; - let rho = a.x.hypot(a.y).max(b.x.hypot(b.y)) + length; if !(lat.is_finite() && rho.is_finite()) || lat >= 89f64.to_radians() { return true; } STRAY_MARGIN * (lat.tan() + FLAT_ALLOWANCE + PLANE_TURN * rho / EARTH_R_M) / (8.0 * EARTH_R_M) * length2 > tolerance } +/// How far, at most, an edge straight in longitude and latitude strays from its chord in a plane +/// about its pole (meters): through `dl` radians of longitude, its distance from the pole changing +/// by `drho` meters, `rho` being its end farther from the pole, it's a spiral, which strays by at +/// most Δλ(ρΔλ + 2Δρ)/8 (the most its second derivative reaches, over 8). For the buffers' edges +/// and the dangle check's pieces beyond 85.05°. +pub(crate) fn spiral_stray(rho: f64, drho: f64, dl: f64) -> f64 { + dl * (rho * dl + 2.0 * drho) / 8.0 +} + /// A longitude difference folded into −180..180. fn wrap(d: f64) -> f64 { (d + 180.0).rem_euclid(360.0) - 180.0 @@ -536,16 +560,99 @@ mod tests { assert_eq!(ring.0.len(), n + MAX_ADDED as usize); } + #[test] + fn an_outline_beyond_85_degrees_follows_its_edges_as_the_tools_draw_them() { + // An octagon 100 km around 87°N in its projection, edges of about 77 km, all beyond 85.05°, + // where the tools draw an edge straight in lon/lat (#1033). Back in lon/lat, each edge of + // the outline, drawn as the tools draw it, stays within the tolerance of the octagon's: it + // gains the points that keep it there. + let plane = Plane::new(Coord::new(10.0, 87.0), 10.0); + let corner = |k: usize| { + let (s, c) = (k as f64 * std::f64::consts::FRAC_PI_4).sin_cos(); + GC { x: 100_000.0 * c, y: 100_000.0 * s } + }; + let octagon = LineString::new((0..=8).map(corner).collect()); + let tolerance = 100_000.0 * OUTLINE_TOLERANCE; + let ring = plane.ring(&octagon, tolerance); + assert!(ring.0.len() > octagon.0.len(), "{}", ring.0.len()); + let worst = farthest_drawn(&plane, &ring, &octagon); + assert!(worst <= 2.0 * tolerance, "{worst} m off the octagon, tolerance {tolerance} m"); + } + + #[test] + fn a_small_buffer_near_a_pole_follows_its_circle_as_the_tools_draw_it() { + // A point's buffer beyond 85.05° has short edges (49 m for 1 km), but near a pole an edge + // straight in lon/lat strays meters even so: each is looked at, not let pass for its length. + // Back in lon/lat, each edge of the outline, drawn as the tools draw it, stays within the + // tolerance of the circle drawn in the plane: across ±180° too, and round a pole (where the + // ring is closed by edges along ±180° and the pole, which aren't the circle's). + let cases = [ + (10.0, 85.5, 1_000.0), + (10.0, 88.0, 1_000.0), + (10.0, 89.9, 100.0), + (10.0, 89.9, 10_000.0), + (10.0, 89.99, 1_000.0), + (10.0, -89.99, 1_000.0), + (180.0, 89.9, 10_000.0), + (10.0, 89.95, 10_000.0), + ]; + for (lon, lat, meters) in cases { + let c = Coord::new(lon, lat); + let plane = Plane::new(c, c.x); + let circle = Point::new(0.0, 0.0).buffer_with_style(buffer_style(meters)); + let Some(outline) = circle.0.first().map(Polygon::exterior) else { panic!("no circle at {lat}°") }; + let tolerance = (meters * OUTLINE_TOLERANCE).max(MIN_TOLERANCE_M); + let ring = plane.ring(outline, tolerance); + // Only the last circle holds its pole. + assert_eq!(ring.0.iter().any(|p| p.y.abs() >= 90.0), (lon, lat) == (10.0, 89.95), "{lon}°, {lat}°"); + let worst = farthest_drawn(&plane, &ring, outline); + assert!(worst <= 2.0 * tolerance, "{lon}°, {lat}°, {meters} m: {worst} m off the circle, tolerance {tolerance} m"); + } + } + + /// How far, at most, `ring` (in lon/lat, as [`Plane::ring`] gives it) drawn as the tools draw + /// its edges (straight in the analysis plane) comes from `outline` in `plane`, at each edge's + /// quarter points and midpoint; edges to a pole, which close a ring round it, aren't measured. + fn farthest_drawn(plane: &Plane, ring: &LineString, outline: &LineString) -> f64 { + let mut worst = 0.0_f64; + for w in ring.0.windows(2) { + let [a, b] = w else { continue }; + if a.y.abs().max(b.y.abs()) >= 90.0 { + continue; + } + let (pa, pb) = (lonlat_to_plane(Coord::new(a.x, a.y)), lonlat_to_plane(Coord::new(b.x, b.y))); + for t in [0.25, 0.5, 0.75] { + let drawn = plane_to_lonlat(Coord::new(pa.x + t * (pb.x - pa.x), pa.y + t * (pb.y - pa.y))); + let p = plane.0.forward(drawn); + let off = outline + .0 + .windows(2) + .filter_map(|e| { + Some(rustgis_data::geometry::seg_dist2( + &[Coord::new(e.first()?.x, e.first()?.y), Coord::new(e.get(1)?.x, e.get(1)?.y)], + p, + )) + }) + .fold(f64::INFINITY, f64::min) + .sqrt(); + worst = worst.max(off); + } + } + worst + } + #[test] fn an_edge_let_pass_by_its_latitude_does_not_stray() { - // Wherever `may_stray` lets an edge pass without looking, looking agrees: planes from 88°S - // to 88°N (past Web Mercator's 85.05° too), edges of 100 m to 1,000 km in eight directions, - // up to 5,000 km from the centre to the side or at 45°, at the tolerances of buffers from - // 100 m to 5,000 km (the largest Buffer makes). - let mut passed = 0; - for lat in (-84..=84).step_by(12).chain([-88, -86, 86, 88]) { - let plane = Plane::new(Coord::new(10.0, f64::from(lat)), 10.0); - for len in [100.0, 500.0, 1_000.0, 2_000.0, 5_000.0, 10_000.0, 30_000.0, 100_000.0, 300_000.0, 1_000_000.0] { + // Wherever `may_stray` lets an edge pass without looking, looking agrees: planes from + // 89.99°S to 89.99°N (past Web Mercator's 85.05° too), edges of 10 m to 1,000 km in eight + // directions, up to 5,000 km from the center to the side or at 45°, at the tolerances of + // buffers from 100 m to 5,000 km (the largest Buffer makes). Beyond 85.05° every edge is + // tried, short ones too, as the spirals' bound lets them pass. + let (mut passed, mut polar) = (0, 0); + let lats = (-84..=84).step_by(12).map(f64::from).chain([-89.99, -89.9, -88.0, -86.0, 85.5, 86.0, 88.0, 89.0, 89.9, 89.99]); + for lat in lats { + let plane = Plane::new(Coord::new(10.0, lat), 10.0); + for len in [10.0, 50.0, 100.0, 500.0, 1_000.0, 2_000.0, 5_000.0, 10_000.0, 30_000.0, 100_000.0, 300_000.0, 1_000_000.0] { for k in 0..8 { let (sin, cos) = (f64::from(k) * std::f64::consts::PI / 8.0).sin_cos(); // The edge's centre off it to the side (where the plane turns it most) or obliquely. @@ -558,10 +665,12 @@ mod tests { let a = GC { x: cx - cos * len / 2.0, y: cy - sin * len / 2.0 }; let b = GC { x: cx + cos * len / 2.0, y: cy + sin * len / 2.0 }; let (la, lb) = (plane.to_lonlat(a), plane.to_lonlat(b)); + let beyond = la.y.abs().min(lb.y.abs()) > rustgis_geom::MAX_LAT; for metres in [100.0, 1_000.0, 5_000.0, 50_000.0, 500_000.0, 2_000_000.0, 5_000_000.0] { let tolerance = (metres * OUTLINE_TOLERANCE).max(MIN_TOLERANCE_M); - if STRAY_PER_M2 * len * len > tolerance && !may_stray(a, b, la, lb, tolerance) { + if (beyond || STRAY_PER_M2 * len * len > tolerance) && !may_stray(a, b, la, lb, tolerance) { passed += 1; + polar += usize::from(beyond); assert!(!plane.strays(a, b, la, lb, tolerance), "{len} m at {lat}°, {off} m off, {metres} m"); } } @@ -569,8 +678,8 @@ mod tests { } } } - // Most edges a buffer has at mid latitudes are let pass. - assert!(passed > 1_000, "{passed}"); + // Most edges a buffer has at mid latitudes are let pass, and many beyond 85.05°. + assert!(passed > 1_000 && polar > 1_000, "{passed}, {polar} beyond 85.05°"); } #[test] diff --git a/crates/analysis/src/lib.rs b/crates/analysis/src/lib.rs index 586abe97a..0b3fbe867 100644 --- a/crates/analysis/src/lib.rs +++ b/crates/analysis/src/lib.rs @@ -1,9 +1,13 @@ //! L2: geoprocessing tools. //! //! Each tool is a [`ToolSpec`]: an id, a toolbox, parameters for the Tools pane, and a -//! pure `run` function over datasets. Inputs refer to layers by name through [`Inputs`]. Work is -//! done in Web Mercator metres, except buffers, which are geodesic: drawn in a projection centred -//! on each feature (see [`geodesic`]) so a 100 km buffer is 100 km on the ground. +//! pure `run` function over datasets. Inputs refer to layers by name through [`Inputs`]. +//! The geometry tools and topology checks work in the analysis plane +//! ([`rustgis_geom::lonlat_to_plane`]: Web Mercator meters within ±85.05°, continued beyond to +//! the poles), radius searches in Web Mercator meters (on the ground when the center or the +//! feature's box reaches beyond 85.05°, see [`network::within`]), and buffers are geodesic: drawn +//! in a projection centered on each feature (see [`geodesic`]) so a 100 km buffer is 100 km on the +//! ground. //! //! Tools from external suites (the GeoLibre/Whitebox tools in `rustgis-geolibre`, and plugins //! in `rustgis-plugins`) join the built-in ones through [`register_external`]: they're listed by [`tools`], found by diff --git a/crates/analysis/src/network.rs b/crates/analysis/src/network.rs index 40dec9598..b51597e72 100644 --- a/crates/analysis/src/network.rs +++ b/crates/analysis/src/network.rs @@ -378,7 +378,8 @@ fn radius_coord_supported(point: Coord) -> bool { /// Local Mercator distance in metres, using the same centre-latitude correction as [`within`]. /// /// This is not a geodesic distance and does not include `within`'s legacy geographic -/// bounding-box prefilter. Inputs must be canonical lon/lat within the Mercator domain. +/// bounding-box prefilter. Inputs must be canonical lon/lat within the Mercator domain. ([`within`] +/// measures a shape reaching beyond 85.05° on the ground instead.) pub fn radius_distance_m(at: Coord, point: Coord) -> Option { if !radius_coord_supported(at) || !radius_coord_supported(point) { return None; @@ -395,7 +396,8 @@ pub fn radius_distance_m(at: Coord, point: Coord) -> Option { /// /// Radii must be positive and at most 5,000 km. Points outside the canonical Mercator /// domain are unsupported rather than wrapped or clamped. The outline is approximate: -/// [`within`] additionally preserves its historical geographic bounding-box prefilter. +/// [`within`] additionally preserves its historical geographic bounding-box prefilter, and it +/// measures a shape reaching beyond 85.05° on the ground, which may then count outside the outline. pub fn radius_boundary_point(at: Coord, radius_m: f64, angle_rad: f64) -> Option { if !radius_coord_supported(at) || !radius_m.is_finite() || radius_m <= 0.0 || radius_m > 5_000_000.0 || !angle_rad.is_finite() { return None; @@ -408,16 +410,60 @@ pub fn radius_boundary_point(at: Coord, radius_m: f64, angle_rad: f64) -> Option radius_coord_supported(point).then_some(point) } -/// Rows of `ds` whose geometry comes within `radius_m` of a point, with their distances. +/// Rows of `ds` whose geometry comes within `radius_m` of a point, with their distances. A shape +/// is measured in local Web Mercator about the point, as [`radius_distance_m`] is; where the point +/// lies beyond 85.05°, where Web Mercator stops, or the shape's box reaches beyond it, it's +/// measured on the ground instead (`polar_reach`) (#1033), to within its steps' sag (up to about +/// 125 m within 5,000 km; farther, more, but only ever nearer). Whether a shape is measured in Web +/// Mercator or on the ground doesn't depend on the radius, and a larger radius never loses a row, +/// to within that sag (though past the search's `REACH_STEPS_TOTAL` steps, which a larger radius +/// may spend sooner, a shape counts at its box's distance; [`within_counted`] says how many did). +/// A center that's no place (a latitude past ±90°, a coordinate missing or infinite), or a radius +/// missing or negative, reaches no row: callers check them first, to say why. pub fn within(ds: &Dataset, rows: &[usize], at: Coord, radius_m: f64) -> Vec<(usize, f64)> { + within_counted(ds, rows, at, radius_m).0 +} + +/// [`within`], and how many of its rows have a distance that's a bound rather than measured: a +/// point or edge as near as its box ([`least_distance_to_box`]), past the search's +/// `REACH_STEPS_TOTAL` steps or along an edge longer than [`REACH_STEPS_MAX`] steps, or a polygon +/// counted as holding the center (at 0) when its turns around the pole outrun the steps or are too +/// many to count. They may lie farther, even beyond the radius. +pub fn within_counted(ds: &Dataset, rows: &[usize], at: Coord, radius_m: f64) -> (Vec<(usize, f64)>, usize) { + use rustgis_geom::MAX_LAT; + // A center that's no place (a latitude past ±90°, a coordinate missing) reaches nothing. + if !(at.x.is_finite() && at.y.abs() <= 90.0 && radius_m >= 0.0) { + return (Vec::new(), 0); + } let m = rustgis_geom::lonlat_to_merc(at); let k = radius_scale(at); + let polar_point = at.y.abs() > MAX_LAT; + // On the ground, its longitude within ±180° (to within a rounding, however far around it's + // written), so every measure there sees the same place; Web Mercator takes it as written. + let lon = at.x.rem_euclid(360.0); + let center = Coord::new(if lon >= 180.0 { lon - 360.0 } else { lon }, at.y); + let ground = rustgis_geom::Aeqd::new(center); + // The steps along polar shapes' edges this search may take, shared by them all. + let mut steps = Steps::new(REACH_STEPS_TOTAL); + let mut bounded = 0; let mut out: Vec<(usize, f64)> = rows .iter() .filter_map(|r| { let g = ds.features.get(*r)?.geometry.as_ref()?; // Skip shapes whose box is farther than the radius before projecting them. let e = g.extent(); + if polar_point || e.min_y < -MAX_LAT || e.max_y > MAX_LAT { + if least_distance_to_box(center, (e.min_x, e.max_x), (e.min_y, e.max_y)) > radius_m { + return None; + } + steps.bound = f64::INFINITY; + let d = polar_reach(g, &e, center, &ground, radius_m, &mut steps); + // A shape with nothing to measure is no distance away, however far the radius. + let found = d <= radius_m && d.is_finite(); + // Its distance a bound: it may lie farther, even beyond the radius. + bounded += usize::from(found && steps.bound <= d); + return found.then_some((*r, d)); + } let dlat = radius_m / 111_320.0; let dlon = dlat * k; if at.x < e.min_x - dlon || at.x > e.max_x + dlon || at.y < e.min_y - dlat || at.y > e.max_y + dlat { @@ -429,7 +475,228 @@ pub fn within(ds: &Dataset, rows: &[usize], at: Coord, radius_m: f64) -> Vec<(us }) .collect(); out.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal)); - out + (out, bounded) +} + +/// The steps a search ([`within`]) has left to take along its shapes; and, for the shape being +/// measured, the least distance it counted at a bound rather than measured: a point or edge as +/// near as its box, once the steps it needed weren't there or were more than one edge takes. +struct Steps { + left: usize, + bound: f64, +} + +impl Steps { + fn new(left: usize) -> Self { + Self { left, bound: f64::INFINITY } + } + + /// Notes a distance counted at a bound. + fn bound(&mut self, d: f64) -> f64 { + self.bound = self.bound.min(d); + d + } + + /// Takes `n` steps if that many are left. + fn take(&mut self, n: usize) -> bool { + let enough = n <= self.left; + if enough { + self.left -= n; + } + enough + } +} + +/// The longest step, in degrees of longitude or latitude, that [`polar_reach`] follows an edge by. +const REACH_STEP_DEG: f64 = 0.5; +/// The most steps it takes along one edge; one needing more (over 512° long, as only a hostile +/// edge is) counts as near as its box. +const REACH_STEPS_MAX: f64 = 1024.0; +/// The most steps one search ([`within`]) takes along all its shapes' edges, a point measured +/// counting as one; past them, a point or edge counts as near as its box +/// ([`least_distance_to_box`]), so a shape may be counted in that's farther than the radius, never +/// left out for want of steps. +const REACH_STEPS_TOTAL: usize = 250_000; +/// More than a degree of longitude or latitude is long anywhere on the ground (meters): a step's +/// length is at most this times the degrees it spans in each. +const DEGREE_M_MAX: f64 = 111_700.0; + +/// The least distance on the ground (meters) from `p` to a point with longitude in `lon` (west, +/// then east) and latitude in `lat` (degrees; a longitude range 360° wide holds every longitude), +/// never more than the true one: the central angle (from its tangent, which rounds little at any +/// distance) on a sphere of WGS 84's least radius of curvature, less a centimeter. An edge +/// straight in the analysis plane stays within its ends' box, so this bounds it too. +fn least_distance_to_box(p: Coord, (lon0, lon1): (f64, f64), (lat0, lat1): (f64, f64)) -> f64 { + if !(p.is_finite() && lon0.is_finite() && lon1.is_finite() && lat0.is_finite() && lat1.is_finite()) { + return 0.0; + } + let (lat0, lat1) = (lat0.min(lat1).clamp(-90.0, 90.0), lat0.max(lat1).clamp(-90.0, 90.0)); + let phi = p.y.clamp(-90.0, 90.0); + // How far `p`'s longitude is from the range, the short way around (none within it). + let (width, past) = (lon1 - lon0, (p.x - lon0).rem_euclid(360.0)); + let dlon = if width >= 360.0 || past <= width { 0.0 } else { (past - width).min(360.0 - past) }; + let angle = if dlon <= 0.0 { + // Along `p`'s own meridian, to the nearer latitude. + (lat0 - phi).max(phi - lat1).max(0.0).to_radians() + } else { + // On the range's nearer edge meridian. Along it, from pole to pole, the distance from `p` + // has one turn: while Δλ is under 90°, a least one at atan(tan φ / cos Δλ), so the nearest + // in range is that latitude clamped; from 90° on, a greatest one, so the nearest in range + // is one of its ends. + let (s, c) = dlon.to_radians().sin_cos(); + let (sin, cos) = phi.to_radians().sin_cos(); + let angle_to = |q: f64| { + let (q_sin, q_cos) = q.to_radians().sin_cos(); + (q_cos * s).hypot(cos * q_sin - sin * q_cos * c).atan2(sin * q_sin + cos * q_cos * c) + }; + if c > 0.0 { angle_to((sin / cos / c).atan().to_degrees().clamp(lat0, lat1)) } else { angle_to(lat0).min(angle_to(lat1)) } + }; + (angle * rustgis_geom::WGS84_MIN_RADIUS_M - 0.01).max(0.0) +} + +/// How far `g` comes from `at` on the ground (meters), as [`within`] measures it near a pole: zero +/// when `at` (its longitude within ±180°, as [`within`] gives it, moved by whole turns into `g`'s +/// box `e`) lies in a polygon in the analysis plane, as the other tools see it; otherwise the +/// least distance to its points or edges, each edge followed as the plane draws it (straight in +/// the plane: in longitude and latitude beyond 85.05°) in steps of at most half a degree, in +/// `ground`, an azimuthal equidistant plane at `at`. Its vertices alone would draw a long edge as +/// a chord, or one along a parallel from −180° to 180° as a point. A point or edge whose box can't +/// come within `radius_m`, or nearer than one already measured, isn't measured; each point +/// measured and each step takes one of the search's `steps`, and once they run out, a point or +/// edge counts as near as its box ([`least_distance_to_box`]). +fn polar_reach(g: &Geometry, e: &rustgis_geom::Extent, at: Coord, ground: &rustgis_geom::Aeqd, radius_m: f64, steps: &mut Steps) -> f64 { + use rustgis_geom::{lonlat_to_plane, plane_to_lonlat}; + let from = |c: Coord| { + let p = ground.forward(c); + (p, p.x.hypot(p.y)) + }; + if let Geometry::Polygon(_) = g { + // The center at its longitude moved by whole turns into the polygon's box: a polygon + // written across ±180°, or from 0° to 360°, holds a center written a turn or more away; one + // whose box holds none of them holds no center. Every turn the box spans is tried: the + // first four (all an ordinary polygon has) freely, each further one taking as many steps as + // the polygon has vertices. When, past the first four, the steps left can't pay for the + // rest, or the turns are too many to count one by one (past 2⁵³), the polygon counts as + // holding the center, a bound, at once, taking none (as an edge too long to step along + // does), so one polygon written a billion degrees long doesn't spend the search's steps. + if e.min_y <= at.y && at.y <= e.max_y { + let (first, last) = (((e.min_x - at.x) / 360.0).ceil(), ((e.max_x - at.x) / 360.0).floor()); + let mut plane = None; + let mut tried = 0_u32; + let mut turn = first; + // The steps each turn past the first four takes. + let mut cost = 1; + while turn <= last { + if tried == 4 { + let mut n = 0_usize; + g.for_each_coord(|_| n += 1); + cost = n.max(1); + // (Finite, or infinite for a box billions of turns wide: never missing.) + if (last - turn + 1.0) * cost as f64 > steps.left as f64 { + return steps.bound(0.0); + } + } + if tried >= 4 && !steps.take(cost) { + return steps.bound(0.0); + } + let c = Coord::new(at.x + 360.0 * turn, at.y); + if plane.get_or_insert_with(|| g.map_coords(&lonlat_to_plane)).hit(lonlat_to_plane(c), 0.0) { + return 0.0; + } + tried = tried.saturating_add(1); + let next = first + f64::from(tried); + if next <= turn { + return steps.bound(0.0); + } + turn = next; + } + } + } + let mut best = f64::INFINITY; + let mut point = |c: Coord, best: &mut f64| { + let near = if c.is_finite() { least_distance_to_box(at, (c.x, c.x), (c.y, c.y)) } else { f64::INFINITY }; + if near <= radius_m && near < *best { + *best = if steps.take(1) { best.min(from(c).1) } else { steps.bound(near) }; + } + }; + // A line's part of one point, or a ring of one point (closed or not), has no edges: it's + // measured as that point. + match g { + Geometry::Point(_) | Geometry::MultiPoint(_) => { + g.for_each_coord(|c| point(c, &mut best)); + return best; + } + Geometry::Line(parts) => parts.iter().for_each(|p| { + if let [c] = p.as_slice() { + point(*c, &mut best); + } + }), + Geometry::Polygon(polys) => polys.iter().flatten().for_each(|r| match r.as_slice() { + [c] => point(*c, &mut best), + [c, d] if c == d => point(*c, &mut best), + _ => {} + }), + } + let origin = Coord::new(0.0, 0.0); + // The last vertex measured: the next edge starts there when they're joined. + let mut joined: Option<(Coord, (Coord, f64))> = None; + g.for_each_edge(|_, a, b| { + // A coordinate missing can't be measured. + if !(a.is_finite() && b.is_finite()) { + return; + } + let near = least_distance_to_box(at, (a.x.min(b.x), a.x.max(b.x)), (a.y.min(b.y), a.y.max(b.y))); + if near > radius_m || near >= best { + return; + } + // Its start, measured already when the edge before ended there, or one more step. + let measured = joined.filter(|(c, _)| *c == a).map(|(_, measured)| measured); + // An edge of no length, or along a pole (as rings cut at ±180° are closed: both ends at the + // North Pole, or both at the South), is one place on the ground, measured once, as a point + // is. + let tolerance = rustgis_geom::POLE_TOLERANCE; + let along = a.y.min(b.y) >= 90.0 - tolerance || a.y.max(b.y) <= -90.0 + tolerance; + if a == b || along { + let start = measured.or_else(|| steps.take(1).then(|| from(a))); + match start { + Some(start) => { + best = best.min(start.1); + joined = Some((b, start)); + } + None => best = best.min(steps.bound(near)), + } + return; + } + // Steps equal in the plane: at most half a degree of longitude, and of latitude where the + // plane stretches latitude least along the edge (at its lowest latitude, or 85.05° beyond). + let (pa, pb) = (lonlat_to_plane(a), lonlat_to_plane(b)); + let low = if a.y * b.y <= 0.0 { 0.0 } else { a.y.abs().min(b.y.abs()).min(rustgis_geom::MAX_LAT) }; + let dlat = ((pb.y - pa.y).abs() / rustgis_geom::MERCATOR_R * low.to_radians().cos()).to_degrees(); + let span = (b.x - a.x).abs().max(dlat) / REACH_STEP_DEG; + let n = if span > 1.0 { span.ceil().min(REACH_STEPS_MAX) as usize } else { 1 }; + let cost = n + usize::from(measured.is_none()); + // An edge needing more steps than one edge takes (over 512° long), or more than are left, + // counts as near as its box: longer steps' chords could pass farther than the edge does. + if span > REACH_STEPS_MAX || !steps.take(cost) { + best = best.min(steps.bound(near)); + return; + } + let start = measured.unwrap_or_else(|| from(a)); + let (mut last, (mut last_p, mut last_d)) = (a, start); + for k in 1..=n { + let t = k as f64 / n as f64; + let q = if k == n { b } else { plane_to_lonlat(Coord::new(pa.x + (pb.x - pa.x) * t, pa.y + (pb.y - pa.y) * t)) }; + let (q_p, q_d) = from(q); + // A step is no nearer than its nearer end less its length: near `at`'s antipode, where + // the plane tears, its chord can pass close to `at` though the step is far. + let chord = rustgis_data::geometry::seg_dist2(&[last_p, q_p], origin).sqrt(); + let length = ((q.x - last.x).abs() + (q.y - last.y).abs()) * DEGREE_M_MAX; + best = best.min(chord.max(last_d.min(q_d) - length)); + (last, last_p, last_d) = (q, q_p, q_d); + } + joined = Some((b, (last_p, last_d))); + }); + best } fn pg_dist(g: &Geometry, p: Coord) -> f64 { @@ -447,6 +714,201 @@ mod tests { use super::*; use rustgis_data::{Field, FieldType, GeometryType, Value}; + #[test] + fn polar_reach_spends_its_budget_then_bounds_edges_by_their_boxes() { + // A line around 89.95°N five times, 5.58 km from a center at 89.9°N: with room for two + // edges' 718 steps and the first's start (the second starts where the first ended), it + // follows two and counts the third as near as its box, 5.53 km, which is never more than + // the line really is. + let c = Coord::new; + let line = Geometry::Line(vec![(0..=5).map(|i| c(if i % 2 == 0 { -179.5 } else { 179.5 }, 89.95)).collect()]); + let at = c(0.0, 89.9); + let mut steps = Steps::new(2 * 718 + 1); + let d = polar_reach(&line, &line.extent(), at, &rustgis_geom::Aeqd::new(at), 10_000.0, &mut steps); + assert!(steps.left == 0 && steps.bound == d); + assert!((5_520.0..5_540.0).contains(&d), "{d}"); + // With room enough, every edge is followed: the line where it is. + let mut steps = Steps::new(10 * 718); + let d = polar_reach(&line, &line.extent(), at, &rustgis_geom::Aeqd::new(at), 10_000.0, &mut steps); + assert!((5_560.0..5_600.0).contains(&d) && steps.left > 0 && steps.bound.is_infinite(), "{d}"); + // A box bound is never more than the distance: the box of an edge along a parallel. + let bound = least_distance_to_box(at, (-179.5, 179.5), (89.95, 89.95)); + assert!(bound <= d && bound > 5_500.0, "{bound}"); + } + + #[test] + fn a_box_is_never_farther_than_its_nearest_point() { + // From 1°N, a box from 150°E to 160°E reaching the South Pole: its nearest point is the + // pole, 91° away (past 90° of longitude away, distances along a meridian peak in between); + // and one at 30°E from 10°N to 20°N, nearest at its corner at 10°N. + let c = Coord::new; + let p = c(0.0, 1.0); + for (lon, lat, nearest) in [((150.0, 160.0), (-90.0, -80.0), c(150.0, -90.0)), ((30.0, 40.0), (10.0, 20.0), c(30.0, 10.0))] { + let (bound, d) = (least_distance_to_box(p, lon, lat), rustgis_geom::wgs84_distance_m(p, nearest)); + assert!(bound <= d && bound > 0.99 * d, "{lon:?} {lat:?}: {bound} vs {d}"); + } + // Micrometers away, where an arccosine would round to centimeters. + let (p, q) = (c(-135.34187, -75.78973), c(-135.34187 + 4e-11, -75.78973)); + assert!(least_distance_to_box(p, (q.x, q.x), (q.y, q.y)) <= rustgis_geom::wgs84_distance_m(p, q)); + } + + #[test] + fn a_center_is_where_its_longitude_falls_however_it_is_written() { + // 1e300° is a whole number of turns: from 87°N there, as from 0°, a polygon from 175°E to + // 185°E is 556 km away over the pole, not around the center. + let c = Coord::new; + let mut ds = Dataset::new("caps", GeometryType::Polygon, vec![]); + ds.push( + Some(Geometry::Polygon(vec![vec![vec![c(175.0, 86.0), c(185.0, 86.0), c(185.0, 88.0), c(175.0, 88.0), c(175.0, 86.0)]]])), + vec![], + ); + let (far, near) = (within(&ds, &[0], c(1e300, 87.0), 600_000.0), within(&ds, &[0], c(0.0, 87.0), 600_000.0)); + assert!(far.len() == 1 && far == near && (550_000.0..560_000.0).contains(&far[0].1), "{far:?} {near:?}"); + assert!(within(&ds, &[0], c(1e300, 87.0), 1.0).is_empty()); + // A center at 180.1° lies in a polygon written from 179.6° to 180.2°, as on the map. + let mut across = Dataset::new("across", GeometryType::Polygon, vec![]); + let ring = vec![c(179.6, 10.0), c(180.2, 10.0), c(180.2, 10.4), c(179.6, 10.4), c(179.6, 10.0)]; + across.push(Some(Geometry::Polygon(vec![vec![ring]])), vec![]); + assert_eq!(within(&across, &[0], c(180.1, 10.2), 1_000.0), vec![(0, 0.0)]); + // Within ±85.05°, Web Mercator takes a center as written, as before: one at 180° finds a + // point at 179.9°, and one at 190° a point written 190.1°. + let mut pts = Dataset::new("pts", GeometryType::Point, vec![]); + pts.push(Some(Geometry::Point(c(179.9, 0.0))), vec![]); + pts.push(Some(Geometry::Point(c(190.1, 0.0))), vec![]); + let found = |at: Coord| within(&pts, &[0, 1], at, 50_000.0).into_iter().map(|(r, _)| r).collect::>(); + assert_eq!((found(c(180.0, 0.0)), found(c(190.0, 0.0))), (vec![0], vec![1])); + } + + #[test] + fn an_edge_along_a_pole_is_measured_once() { + // From 180° to −180° along the South Pole, as rings cut at ±180° are closed: one place, one + // step of the budget, 1° (111.7 km) from 89°S. + let c = Coord::new; + let (at, line) = (c(0.0, -89.0), Geometry::Line(vec![vec![c(180.0, -90.0), c(-180.0, -90.0)]])); + let mut steps = Steps::new(10); + let d = polar_reach(&line, &line.extent(), at, &rustgis_geom::Aeqd::new(at), 200_000.0, &mut steps); + assert!(steps.left == 9 && (d - rustgis_geom::wgs84_distance_m(at, c(0.0, -90.0))).abs() < 1e-3, "{} {d}", steps.left); + // Going on up 180°W from there, the next edge starts where that one ended: one step more. + let line = Geometry::Line(vec![vec![c(180.0, -90.0), c(-180.0, -90.0), c(-180.0, -89.75)]]); + let mut steps = Steps::new(10); + polar_reach(&line, &line.extent(), at, &rustgis_geom::Aeqd::new(at), 200_000.0, &mut steps); + assert_eq!(steps.left, 8); + // A ring whose vertices are all one point is that point, once, as is a vertex repeated + // beside a coordinate missing: 10 km off, it's found. + let p = c(0.0, -89.09); + let ring = Geometry::Polygon(vec![vec![vec![p, p, p, p]]]); + let mut steps = Steps::new(10); + let d = polar_reach(&ring, &ring.extent(), at, &rustgis_geom::Aeqd::new(at), 200_000.0, &mut steps); + assert!(steps.left == 9 && (d - rustgis_geom::wgs84_distance_m(at, p)).abs() < 1e-3, "{} {d}", steps.left); + let mut beside = Dataset::new("beside", GeometryType::Polyline, vec![]); + beside.push(Some(Geometry::Line(vec![vec![p, p, c(f64::NAN, -89.5)]])), vec![]); + let found = within(&beside, &[0], at, 20_000.0); + assert!(found.len() == 1 && found.iter().all(|(_, d)| (d - rustgis_geom::wgs84_distance_m(at, p)).abs() < 1e-3), "{found:?}"); + // A shape with nothing to measure isn't found, however far the radius reaches. + let mut empty = Dataset::new("empty", GeometryType::Polyline, vec![]); + empty.push(Some(Geometry::Line(vec![])), vec![]); + assert!(within(&empty, &[0], at, f64::INFINITY).is_empty()); + } + + #[test] + fn polar_reach_counts_points_against_its_budget() { + // Points at 89.95°N, 16.7, 12.5 and 5.6 km from a center at 89.9°N, nearest last: each + // measured takes one of the budget; with none left, a point counts as near as its box, + // never more than it is; and one out of reach takes none. + let c = Coord::new; + let at = c(0.0, 89.9); + let (ground, pts) = (rustgis_geom::Aeqd::new(at), Geometry::MultiPoint(vec![c(180.0, 89.95), c(90.0, 89.95), c(0.0, 89.95)])); + let mut steps = Steps::new(10); + let d = polar_reach(&pts, &pts.extent(), at, &ground, 50_000.0, &mut steps); + assert!(steps.left == 7 && steps.bound.is_infinite() && (5_570.0..5_600.0).contains(&d), "{} {d}", steps.left); + let mut steps = Steps::new(1); + let bound = polar_reach(&pts, &pts.extent(), at, &ground, 50_000.0, &mut steps); + assert!(steps.left == 0 && steps.bound == bound && bound <= d && bound > 5_500.0, "{bound}"); + let mut steps = Steps::new(10); + let out_of_reach = polar_reach(&pts, &pts.extent(), at, &ground, 1_000.0, &mut steps); + assert!(out_of_reach.is_infinite() && steps.left == 10 && steps.bound.is_infinite()); + } + + #[test] + fn a_center_is_tried_at_every_turn_a_polygon_spans() { + // A strip along 86°N from 0° to 1441°, rising into a block from 1430° to 1441°: a center at + // 86.5°N on 0° lies in the block, 1440° on, its fifth turn. A sliver a billion degrees long, + // rising from 86°N to 87°N along it, reaches 86.5°N only some 1.4 million turns on: it isn't + // tried at every turn for ever, but past the steps counts as holding the center, a bound, + // and says so. + let c = Coord::new; + let mut ds = Dataset::new("strips", GeometryType::Polygon, vec![]); + let block = vec![c(0.0, 86.0), c(1441.0, 86.0), c(1441.0, 87.0), c(1430.0, 87.0), c(1430.0, 86.1), c(0.0, 86.1), c(0.0, 86.0)]; + ds.push(Some(Geometry::Polygon(vec![vec![block]])), vec![]); + let long = vec![c(0.0, 86.0), c(1e9, 86.0), c(1e9, 87.0), c(0.0, 86.0)]; + ds.push(Some(Geometry::Polygon(vec![vec![long]])), vec![]); + assert_eq!(within_counted(&ds, &[0], c(0.0, 86.5), 1_000.0), (vec![(0, 0.0)], 0)); + assert_eq!(within_counted(&ds, &[1], c(0.0, 86.5), 50_000.0), (vec![(1, 0.0)], 1)); + // Past its first four turns, the sliver counts at its bound at once, taking no steps: a + // line after it, about 160 km off, is still measured, and left out of 50 km. + ds.push(Some(Geometry::Line(vec![vec![c(-170.0, 86.0), c(170.0, 89.9)]])), vec![]); + assert_eq!(within_counted(&ds, &[1, 2], c(0.0, 86.5), 50_000.0), (vec![(1, 0.0)], 1)); + // A center on the sliver's corner is found there, at its first turn, exactly. + assert_eq!(within_counted(&ds, &[1], c(0.0, 86.0), 50_000.0), (vec![(1, 0.0)], 0)); + // A wedge written 1e19° on, its turns some 2.8e16, past 2⁵³, where a turn and the next + // round to one number: the center isn't in it at the first, which can't be counted on from, + // so it's a bound, at once, taking no steps. + let wedge = Geometry::Polygon(vec![vec![vec![c(1e19, 86.0), c(1e19 + 1e6, 86.0), c(1e19 + 1e6, 86.6), c(1e19, 86.0)]]]); + let mut steps = Steps::new(REACH_STEPS_TOTAL); + let at = c(0.0, 86.5); + assert!((1e19 / 360.0f64).ceil() + 1.0 == (1e19 / 360.0f64).ceil()); + let d = polar_reach(&wedge, &wedge.extent(), at, &rustgis_geom::Aeqd::new(at), 50_000.0, &mut steps); + assert!(d == 0.0 && steps.bound == 0.0 && steps.left == REACH_STEPS_TOTAL, "{d} {} {}", steps.bound, steps.left); + // A center in the polygon at its first turn is found there, exactly, however many turns + // its box spans after. + let mut wide = Dataset::new("wide", GeometryType::Polygon, vec![]); + let ring = vec![c(-10.0, 86.0), c(1e9, 86.0), c(1e9, 87.0), c(-10.0, 87.0), c(-10.0, 86.0)]; + wide.push(Some(Geometry::Polygon(vec![vec![ring]])), vec![]); + assert_eq!(within_counted(&wide, &[0], c(0.0, 86.5), 1_000.0), (vec![(0, 0.0)], 0)); + } + + #[test] + fn a_distance_counted_at_a_bound_says_so() { + // From 89.9°N, a line along 89.95°N: its first, short edge is measured 5.75 km away; its + // next, 180° long, needs more steps than are left and counts at its box's distance, nearer: + // the line's distance is that bound. With steps enough, it's measured, and no bound. An + // edge 600° long, needing more steps than one edge takes, counts at its box, taking none + // (in longer steps, chords could pass farther than the edge does). + let c = Coord::new; + let (at, ground) = (c(0.0, 89.9), rustgis_geom::Aeqd::new(c(0.0, 89.9))); + let line = Geometry::Line(vec![vec![c(10.0, 89.95), c(10.1, 89.95), c(-170.0, 89.95)]]); + let mut steps = Steps::new(5); + let d = polar_reach(&line, &line.extent(), at, &ground, 10_000.0, &mut steps); + assert!(steps.bound == d && d < 5_700.0, "{d}"); + let mut steps = Steps::new(REACH_STEPS_TOTAL); + let d = polar_reach(&line, &line.extent(), at, &ground, 10_000.0, &mut steps); + assert!(steps.bound.is_infinite() && d > 5_500.0, "{d}"); + let long = Geometry::Line(vec![vec![c(0.0, 89.95), c(600.0, 89.95)]]); + let mut steps = Steps::new(REACH_STEPS_TOTAL); + let d = polar_reach(&long, &long.extent(), at, &ground, 10_000.0, &mut steps); + let bound = least_distance_to_box(at, (0.0, 600.0), (89.95, 89.95)); + assert!(steps.bound == d && d == bound && steps.left == REACH_STEPS_TOTAL, "{d} {}", steps.bound); + } + + #[test] + fn a_search_counts_the_rows_it_found_past_its_steps() { + // 300 lines each going 348 times around 89.95°N, 5.6 km from a center at 89.9°N: more + // steps than one search has, so the last rows count at their boxes' distance, and say so; + // a search at 85°N, measured in Web Mercator, takes no steps. + let c = Coord::new; + let mut ds = Dataset::new("lines", GeometryType::Polyline, vec![]); + for _ in 0..300 { + ds.push(Some(Geometry::Line(vec![(0..=348).map(|i| c(if i % 2 == 0 { -179.5 } else { 179.5 }, 89.95)).collect()])), vec![]); + } + let rows: Vec = (0..300).collect(); + let (found, bounded) = within_counted(&ds, &rows, c(0.0, 89.9), 10_000.0); + assert!(found.len() == 300 && bounded > 0 && bounded < 300, "{} {bounded}", found.len()); + let mut near = Dataset::new("near", GeometryType::Point, vec![]); + near.push(Some(Geometry::Point(c(0.0, 85.0))), vec![]); + let (found, bounded) = within_counted(&near, &[0], c(0.01, 85.0), 10_000.0); + assert!(found.len() == 1 && bounded == 0, "{found:?} {bounded}"); + } + #[test] fn routes_points_profiles_and_similarity() { // A square of roads with a slow diagonal. @@ -552,8 +1014,8 @@ mod tests { assert_eq!(radius_distance_m(at, at), Some(0.0)); let boundary = radius_boundary_point(at, 1000.0, 0.0).unwrap(); assert!((rustgis_geom::lonlat_to_merc(boundary).x - m.x - 2000.0).abs() < 1e-8); - for latitude in [0.0, 60.0] { - let center = Coord::new(3.0, latitude); + for (longitude, latitude) in [(3.0, 0.0), (3.0, 60.0), (-73.98, 40.75)] { + let center = Coord::new(longitude, latitude); for angle in [0.0, 0.5, 1.0, 2.0, 3.0, 4.5] { let point = radius_boundary_point(center, 1000.0, angle).unwrap(); assert!((radius_distance_m(center, point).unwrap() - 1000.0).abs() < 1e-7); diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index a94161f08..a2a48ba81 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -13,7 +13,7 @@ use std::borrow::Cow; use std::collections::{HashMap, HashSet, TryReserveError}; use std::ops::ControlFlow; -use crate::convert::{from_line, from_multipolygon, from_point, to_geo, to_multiline, to_multipolygon}; +use crate::convert::{from_line, from_multipolygon, from_point, reaches_poles, to_geo, to_multiline, to_multipolygon}; use crate::geodesic::{PARALLEL_FEATURES, VERTEX_BUDGET_ERROR, geodesic_buffer_with_limit, geodesic_buffers_with_limit}; use crate::{Inputs, Output, OutputLayer, ParamKind as K, ParamSpec as P, ToolResult, ToolSpec}; @@ -1232,7 +1232,10 @@ fn out_name(p: &Json, default: String) -> String { if n.is_empty() { default } else { n.to_string() } } -/// Scale factor of Web Mercator at a latitude (metres on the map per metre on the ground). +/// Scale factor of Web Mercator at a latitude (meters on the map per meter on the ground), clamped +/// at 85°. It scales tolerances for the analysis plane too: beyond 85.05° the plane stretches +/// north–south by sec(85.05°) and east–west by sec(latitude), both more than sec(85°), so a +/// tolerance converted with it reaches no farther on the ground than asked there (#1033). fn merc_scale(lat: f64) -> f64 { 1.0 / lat.clamp(-85.0, 85.0).to_radians().cos().max(0.01) } @@ -1497,9 +1500,10 @@ fn dissolved(g: &Geometry) -> std::borrow::Cow<'_, Geometry> { /// overlap with it, a line's stretches in it, the points in it or on its edge; outside, the rest. /// `None` when nothing is left, or only a shared edge or corner. /// -/// The mask's polygons are in Web Mercator, oriented as `union_all` leaves them. The feature's may -/// wind either way (the reference corpus has a clockwise parcel), but must not overlap each other: -/// pass a feature of several polygons through [`dissolved`] first, once for all its masks. +/// The mask's polygons are in the analysis plane, oriented as `union_all` leaves them. The +/// feature's may wind either way (the reference corpus has a clockwise parcel), but must not +/// overlap each other: pass a feature of several polygons through [`dissolved`] first, once for all +/// its masks. fn part(g: &Geometry, mask: &MultiPolygon, inside: bool) -> Option { match g { Geometry::Polygon(_) => { @@ -1512,12 +1516,19 @@ fn part(g: &Geometry, mask: &MultiPolygon, inside: bool) -> Option v, _ => &[], }; + // The poles the mask reaches, found at the first point at one. + let mut reached = None; let kept: Vec = pts .iter() .copied() .filter(|c| { - let m = rustgis_geom::lonlat_to_merc(*c); - mask.intersects(&geo::Point::new(m.x, m.y)) == inside + let m = rustgis_geom::lonlat_to_plane(*c); + // A point at a pole touches a mask that reaches the pole, whatever its + // longitude. + let at = reaches_poles(c.y, c.y); + let touches = mask.intersects(&geo::Point::new(m.x, m.y)) + || ((at.0 || at.1) && meet_at_a_pole(at, *reached.get_or_insert_with(|| plane_poles(mask)))); + touches == inside }) .collect(); match kept.len() { @@ -1551,6 +1562,8 @@ fn clip(inp: &dyn Inputs, p: &Json) -> Result { .filter(|(m, _)| !m.0.is_empty()) .collect(); let index = ExtentIndex::new(masks.iter().map(|(_, ext)| *ext)); + // The masks reaching each pole: a point at one is in them at any longitude (#1033). + let at_poles = AtPoles::new(masks.iter().map(|(m, _)| plane_poles(m))); // Single candidates reuse their normalized mask. Retain eight recent multi-mask unions // so alternating neighbourhoods can reuse work too; evict the least recently used. const MAX_CACHED_MASKS: usize = 8; @@ -1566,7 +1579,11 @@ fn clip(inp: &dyn Inputs, p: &Json) -> Result { let domain = index.total.union(&ext); let margin = domain.width().max(domain.height()) / f64::from(1_u32 << 27); let ext = ext.padded(0.0, margin); - let near: Vec<_> = index.candidates(&ext).into_iter().filter(|i| masks.get(*i).is_some_and(|(_, e)| e.intersects(&ext))).collect(); + let mut near: Vec<_> = + index.candidates(&ext).into_iter().filter(|i| masks.get(*i).is_some_and(|(_, e)| e.intersects(&ext))).collect(); + if matches!(g, Geometry::Point(_) | Geometry::MultiPoint(_)) { + at_poles.add_to(poles(&g.extent()), &mut near); + } if near.is_empty() { continue; } @@ -1607,11 +1624,12 @@ fn erase(inp: &dyn Inputs, p: &Json) -> Result { } let (er, zones) = overlay_zones(inp, p, "erase")?; let index = ExtentIndex::new(zones.iter().map(|z| z.ext)); + let at_poles = AtPoles::new(zones.iter().map(|z| poles(&z.ext))); let mut out = copy_schema(ds, out_name(p, format!("{}_Erase", ds.name)), ds.geometry_type); keep_source_crs(&mut out, ds, er); for r in rows(inp, p, "input", ds) { let Some(f) = ds.features.get(r) else { continue }; - if let Some(rest) = f.geometry.as_ref().and_then(|g| outside(g, &zones, &index)) { + if let Some(rest) = f.geometry.as_ref().and_then(|g| outside(g, &zones, (&index, &at_poles))) { out.push_next(Some(rest), f.attrs.clone()); } } @@ -1619,8 +1637,8 @@ fn erase(inp: &dyn Inputs, p: &Json) -> Result { result(out, format!("Erase: {n} features remain")) } -/// An overlay polygon: its row in the overlay layer, its extent and its polygons in Web Mercator, -/// oriented for the overlay, its parts unioned when it has several (as Clip's mask is). +/// An overlay polygon: its row in the overlay layer, its extent and its polygons in the analysis +/// plane, oriented for the overlay, its parts unioned when it has several (as Clip's mask is). struct Zone { row: usize, ext: Extent, @@ -1728,10 +1746,9 @@ impl Overlay { /// zone order; and, when `outside` is set, its part outside them all. Parts that only touch are /// left out. A feature no zone overlaps is its own part outside, coordinates untouched; one that /// is cut comes out rebuilt, its overlapping parts unioned (see [`dissolved`]). -fn cut(g: &Geometry, zones: &[Zone], index: &ExtentIndex, outside: bool) -> (Vec<(usize, Geometry)>, Option) { +fn cut(g: &Geometry, zones: &[Zone], index: (&ExtentIndex, &AtPoles), outside: bool) -> (Vec<(usize, Geometry)>, Option) { let ext = g.extent(); - let meets: Vec = - index.candidates(&ext).into_iter().filter(|zi| zones.get(*zi).is_some_and(|z| ext.intersects(&z.ext))).collect(); + let meets = meeting(g, &ext, zones, index); // A feature of overlapping parts is unioned once, for all the zones it may meet, and only // when it may meet one. let whole = if meets.is_empty() { std::borrow::Cow::Borrowed(g) } else { dissolved(g) }; @@ -1744,6 +1761,18 @@ fn cut(g: &Geometry, zones: &[Zone], index: &ExtentIndex, outside: bool) -> (Vec (inside, rest) } +/// The zones that may meet a feature (`ext` its extent), in zone order: those whose extent meets +/// its, and, for points, those reaching a pole a point is at, where they meet at any longitude +/// (#1033; a line or polygon meeting a zone only at a pole has no part in it). +fn meeting(g: &Geometry, ext: &Extent, zones: &[Zone], (index, at_poles): (&ExtentIndex, &AtPoles)) -> Vec { + let mut found: Vec = + index.candidates(ext).into_iter().filter(|zi| zones.get(*zi).is_some_and(|z| ext.intersects(&z.ext))).collect(); + if matches!(g, Geometry::Point(_) | Geometry::MultiPoint(_)) { + at_poles.add_to(poles(ext), &mut found); + } + found +} + /// A feature's part outside the zones that cut it, for [`cut`] and [`outside`], which find those /// zones differently: `cut` makes each part inside anyway, so a zone cuts when that part is /// something; `outside` makes none, so it tests each zone with a predicate. `whole` is the feature @@ -1763,9 +1792,9 @@ fn rest_outside(g: &Geometry, whole: &Geometry, hit: &[&Zone]) -> Option Option { +fn outside(g: &Geometry, zones: &[Zone], index: (&ExtentIndex, &AtPoles)) -> Option { let ext = g.extent(); - let near: Vec<&Zone> = index.candidates(&ext).into_iter().filter_map(|zi| zones.get(zi)).filter(|z| ext.intersects(&z.ext)).collect(); + let near: Vec<&Zone> = meeting(g, &ext, zones, index).into_iter().filter_map(|zi| zones.get(zi)).collect(); if near.is_empty() { return Some(g.clone()); } @@ -1773,13 +1802,15 @@ fn outside(g: &Geometry, zones: &[Zone], index: &ExtentIndex) -> Option { let m = z.poly.relate(&geo_g); m.is_intersects() && !m.is_touches() } - _ => z.poly.intersects(&geo_g), + Geometry::Point(_) | Geometry::MultiPoint(_) => meet_at_a_pole(poles(&ext), poles(&z.ext)) || z.poly.intersects(&geo_g), + Geometry::Line(_) => z.poly.intersects(&geo_g), }; let hit: Vec<&Zone> = near.into_iter().filter(cuts).collect(); rest_outside(g, &whole, &hit) @@ -1804,6 +1835,7 @@ fn overlay(inp: &dyn Inputs, p: &Json, op: Overlay) -> Result Result Result { }) .collect::>, String>>()?; let index = ExtentIndex::new(zones.iter().map(|z| z.ext)); + let at_poles = AtPoles::new(zones.iter().map(|z| poles(&z.ext))); let mut out = copy_schema(ds, out_name(p, format!("{}_Update", ds.name)), GeometryType::Polygon); keep_source_crs(&mut out, ds, up); for r in rows(inp, p, "input", ds) { @@ -1940,7 +1973,7 @@ fn update(inp: &dyn Inputs, p: &Json) -> Result { out.push_next(None, f.attrs.clone()); continue; } - Some(g) => outside(g, &zones, &index), + Some(g) => outside(g, &zones, (&index, &at_poles)), }; if let Some(rest) = rest { out.push_next(Some(rest), f.attrs.clone()); @@ -2029,6 +2062,11 @@ fn split(inp: &dyn Inputs, p: &Json) -> Result { let located: Vec<(usize, Extent)> = rows(inp, p, "input", ds).into_iter().filter_map(|r| Some((r, ds.features.get(r)?.geometry.as_ref()?.extent()))).collect(); let index = ExtentIndex::new(located.iter().map(|(_, e)| *e)); + // The input points at each pole: in a split polygon reaching it, at any longitude (#1033). + let point = |r: usize| { + ds.features.get(r).and_then(|f| f.geometry.as_ref()).is_some_and(|g| matches!(g, Geometry::Point(_) | Geometry::MultiPoint(_))) + }; + let points_at_poles = AtPoles::new(located.iter().map(|(r, e)| if point(*r) { poles(e) } else { (false, false) })); // Each input feature as the overlay takes it, made once, when a split polygon first meets it. let mut whole: Vec>> = vec![None; located.len()]; let mut layers = Vec::new(); @@ -2041,13 +2079,16 @@ fn split(inp: &dyn Inputs, p: &Json) -> Result { _ => std::borrow::Cow::Owned(union_all(polys.iter().flat_map(|z| z.poly.0.iter()))), }; let ext = polys.iter().fold(Extent::EMPTY, |e, z| e.union(&z.ext)); + let reached = polys.iter().map(|z| poles(&z.ext)).fold((false, false), |(n, s), (zn, zs)| (n || zn, s || zs)); let mut out = copy_schema(ds, format!("{prefix}_{name}"), ds.geometry_type); keep_source_crs(&mut out, ds, sp); index.candidates_into(&ext, &mut cand_buf); + points_at_poles.add_to(reached, &mut cand_buf); for &li in &cand_buf { let Some((r, e)) = located.get(li) else { continue }; let Some(f) = ds.features.get(*r) else { continue }; - let (Some(g), Some(slot)) = (f.geometry.as_ref().filter(|_| e.intersects(&ext)), whole.get_mut(li)) else { continue }; + let meets = e.intersects(&ext) || (point(*r) && meet_at_a_pole(poles(e), reached)); + let (Some(g), Some(slot)) = (f.geometry.as_ref().filter(|_| meets), whole.get_mut(li)) else { continue }; if let Some(c) = part(slot.get_or_insert_with(|| dissolved(g)), &mask, true) { out.push_next(Some(c), f.attrs.clone()); } @@ -2132,7 +2173,7 @@ fn convex_hull(inp: &dyn Inputs, p: &Json) -> Result { let Some(g) = &f.geometry else { continue }; let mut these = Vec::new(); g.for_each_coord(|c| { - let m = rustgis_geom::lonlat_to_merc(c); + let m = rustgis_geom::lonlat_to_plane(c); these.push(geo::Point::new(m.x, m.y)); }); if each { @@ -3033,18 +3074,36 @@ fn simplify(inp: &dyn Inputs, p: &Json) -> Result { let ds = input(inp, p, "input")?; let tol = num(p, "tolerance").ok_or("Tolerance must be a number")?.abs() * 1000.0; let mut out = copy_schema(ds, out_name(p, format!("{}_Simplified", ds.name)), ds.geometry_type); - let (mut before, mut after) = (0usize, 0usize); + let (mut before, mut after, mut polar) = (0usize, 0usize, 0usize); for r in rows(inp, p, "input", ds) { let Some(f) = ds.features.get(r) else { continue }; let Some(g) = &f.geometry else { continue }; before += g.coord_count(); + // A line or polygon reaching beyond 85.05° has its tolerance taken east–west at the scale + // of its label point's latitude, at most 85° (`merc_scale`), where the ground beyond is far + // narrower: less is taken out there than asked (#1033). Points aren't simplified. + let mut beyond = false; + if !matches!(g, Geometry::Point(_) | Geometry::MultiPoint(_)) { + g.for_each_coord(|c| beyond |= c.y.is_finite() && c.y.abs() > rustgis_geom::MAX_LAT); + } + polar += usize::from(beyond && tol > 0.0); let sg = simplify_geometry(g, tol); after += sg.coord_count(); out.push(Some(sg), f.attrs.clone()); } - result(out, format!("Simplify: {before} → {after} vertices")) + let mut messages = vec![format!("Simplify: {before} → {after} vertices")]; + if polar > 0 { + messages.push(format!( + "{polar} line or polygon feature(s) reach beyond 85.05°, where the tolerance is taken east–west at the scale of the feature's middle, at most 85°'s: less is simplified away there than asked" + )); + } + Ok(ToolResult { output: Output::Layer(out), messages }) } +/// A Select by Location selector: its shape and extent in the analysis plane, and the poles it +/// reaches (north, south). +type Selector = (geo::Geometry, Extent, (bool, bool)); + /// Conservative heap accounting for the geometry variants produced by `to_geo`. Coordinate /// capacities are exact; private polygon-interior arrays use an upper allocation allowance. fn selector_heap_bytes(g: &geo::Geometry) -> usize { @@ -3096,7 +3155,7 @@ fn select_by_location_with_work_limit(inp: &dyn Inputs, p: &Json, limit: usize, // Each selector has its own complexity limit. Retained selectors share a separate byte // budget; zero-distance selection keeps its existing unbuffered behavior. let storage = BufferStorage::new(if dist > 0.0 { bytes } else { usize::MAX }); - let mut selectors = Vec::new(); + let mut selectors: Vec = Vec::new(); crate::control::begin(picked.len()); for batch in picked.chunks(BUFFER_BATCH) { let made = crate::parallel::try_map( @@ -3105,12 +3164,13 @@ fn select_by_location_with_work_limit(inp: &dyn Inputs, p: &Json, limit: usize, || (), |_, _, g| { let b = if dist > 0.0 { geodesic_buffer_with_limit(g, dist, limit)?.geometry } else { None }; - let gg = to_geo(b.as_ref().unwrap_or(g)); + let g = b.as_ref().unwrap_or(g); + let gg = to_geo(g); if dist > 0.0 { storage.charge(selector_heap_bytes(&gg))?; } let ext = merc_extent(&gg); - Ok::<_, &'static str>((gg, ext)) + Ok::<_, &'static str>((gg, ext, poles(&g.extent()))) }, )?; storage.reserve(&mut selectors, made.len())?; @@ -3121,9 +3181,11 @@ fn select_by_location_with_work_limit(inp: &dyn Inputs, p: &Json, limit: usize, if dist > 0.0 { storage.charge(selectors.len().saturating_mul(std::mem::size_of::()))?; } - let index = ExtentIndex::new(selectors.iter().map(|(_, ext)| *ext)); + let index = ExtentIndex::new(selectors.iter().map(|(_, ext, _)| *ext)); + // The selectors reaching each pole, where they meet whatever else reaches it (#1033). + let at_poles = AtPoles::new(selectors.iter().map(|(_, _, at)| *at)); if dist > 0.0 { - storage.charge(index.storage_bytes())?; + storage.charge(index.storage_bytes().saturating_add(at_poles.storage_bytes()))?; } let invert = s(p, "invert") == "Yes"; let bounded_checks = dist > 0.0 && selectors.len() > ExtentIndex::MAX_LINEAR_FEATURES; @@ -3137,17 +3199,18 @@ fn select_by_location_with_work_limit(inp: &dyn Inputs, p: &Json, limit: usize, let gg = to_geo(g); let ext = merc_extent(&gg); let mut exhausted = false; - let hit = index.any_candidate(&ext, |i| { - if bounded_checks { - if remaining_checks == 0 { - exhausted = true; - // Stop the candidate traversal; this sentinel is never used as a match. - return true; + let hit = at_poles.any(poles(&g.extent())) + || index.any_candidate(&ext, |i| { + if bounded_checks { + if remaining_checks == 0 { + exhausted = true; + // Stop the candidate traversal; this sentinel is never used as a match. + return true; + } + remaining_checks -= 1; } - remaining_checks -= 1; - } - selectors.get(i).is_some_and(|(s, e)| e.intersects(&ext) && s.intersects(&gg)) - }); + selectors.get(i).is_some_and(|(s, e, _)| e.intersects(&ext) && s.intersects(&gg)) + }); if exhausted { return Err(format!("Distance selection exceeds its limit of {checks} selector checks; select fewer features or selectors")); } @@ -3470,6 +3533,21 @@ fn kernel_density_tool(inp: &dyn Inputs, p: &Json) -> Result let mut pts: Vec<(Coord, f64)> = Vec::new(); let mut skipped: usize = 0; + // A point the system can't hold is left out, and counted, as one past 85.05° in Web Mercator, + // which would land on the map's edge rather than where it is (#1033): its latitude doesn't + // come back. + let mut beyond: usize = 0; + let project = |c: Coord, beyond: &mut usize| { + // Web Mercator's limit is known exactly: the round trip alone would let a latitude a hair + // past it through, onto the edge. + let held = crs != Crs::WebMercator || c.y.abs() <= rustgis_geom::MAX_LAT; + let p = held + .then(|| transformer.from_lonlat(c)) + .flatten() + .filter(|p| p.is_finite() && transformer.to_lonlat(*p).is_some_and(|back| (back.y - c.y).abs() < 1e-4)); + *beyond += usize::from(p.is_none()); + p + }; for r in rows(inp, p, "input", ds) { let Some(g) = ds.features.get(r).and_then(|f| f.geometry.as_ref()) else { skipped += 1; @@ -3488,21 +3566,17 @@ fn kernel_density_tool(inp: &dyn Inputs, p: &Json) -> Result match g { Geometry::Point(c) => { - if c.is_finite() - && let Some(p_proj) = transformer.from_lonlat(*c) - && p_proj.is_finite() - { - pts.push((p_proj, w)); - } else { + if !c.is_finite() { skipped += 1; + } else if let Some(p_proj) = project(*c, &mut beyond) { + pts.push((p_proj, w)); } } Geometry::MultiPoint(cs) => { let mut any = false; for c in cs { if c.is_finite() - && let Some(p_proj) = transformer.from_lonlat(*c) - && p_proj.is_finite() + && let Some(p_proj) = project(*c, &mut beyond) { pts.push((p_proj, w)); any = true; @@ -3518,8 +3592,12 @@ fn kernel_density_tool(inp: &dyn Inputs, p: &Json) -> Result } } + // Only Web Mercator's limit is known to be a latitude. + let why = if crs == Crs::WebMercator { " (past 85.05° of latitude)" } else { "" }; + let beyond_note = (beyond > 0).then(|| format!("{beyond} point(s) beyond what {} can hold{why} were left out", crs.label())); if pts.is_empty() { - return Err("input layer has no valid points with positive weights to analyze".into()); + let none = "input layer has no valid points with positive weights to analyze"; + return Err(beyond_note.map_or_else(|| none.to_string(), |note| format!("{none}: {note}"))); } let pt_extent = crate::density::points_extent(&pts).ok_or("input layer has no valid points with positive weights to analyze")?; @@ -3582,7 +3660,7 @@ fn kernel_density_tool(inp: &dyn Inputs, p: &Json) -> Result }; Ok(ToolResult { output: Output::Layers(vec![OutputLayer::Raster { file_name: format!("{name}.tif"), name, bytes }]), - messages: vec![msg], + messages: std::iter::once(msg).chain(beyond_note).collect(), }) } @@ -3590,11 +3668,11 @@ fn kernel_density_tool(inp: &dyn Inputs, p: &Json) -> Result /// same transform and non-finite-coordinate filtering as [`to_geo`] followed by [`merc_extent`]. fn projected_extent(g: &Geometry) -> Extent { let mut ext = Extent::EMPTY; - g.for_each_coord(|c| ext.expand(rustgis_geom::lonlat_to_merc(c))); + g.for_each_coord(|c| ext.expand(rustgis_geom::lonlat_to_plane(c))); ext } -/// The extent of a geometry in Web Mercator, as [`to_geo`] gives it, leaving out non-finite +/// The extent of a geometry in the analysis plane, as [`to_geo`] gives it, leaving out non-finite /// coordinates as [`Extent::expand`] does: one NaN doesn't hide the rest. fn merc_extent(g: &impl geo::CoordsIter) -> Extent { g.coords_iter().fold(Extent::EMPTY, |mut e, c| { @@ -3605,10 +3683,11 @@ fn merc_extent(g: &impl geo::CoordsIter) -> Extent { /// Summarize Within: each polygon's count of the (selected) point features that meet it, inside it /// or on its edge; a feature of several points counts once in each polygon it meets. Each point is -/// tested only against the polygons whose extents hold it, through an [`ExtentIndex`]. The extents -/// are in Web Mercator, where the test runs, so a point beyond the map's latitude limit still -/// meets a polygon reaching that limit. The statistics asked for, of the point layer's fields, -/// are over the same points, as `_` after `point_count`. +/// tested only against the polygons whose extents hold it, through an [`ExtentIndex`], and, at a +/// pole, those reaching it, which it touches at any longitude. The extents are in the analysis +/// plane, where the test runs, which goes on past 85.05° to the poles (#1033). The statistics +/// asked for, of the point layer's fields, are over the same points, as `_` after +/// `point_count`. fn count_points(inp: &dyn Inputs, p: &Json) -> Result { let polys = input(inp, p, "polygons")?; let pts = input(inp, p, "points")?; @@ -3618,18 +3697,21 @@ fn count_points(inp: &dyn Inputs, p: &Json) -> Result { out.source = None; let fi = out.ensure_field("point_count", FieldType::Integer); let added = statistic_fields(pts, &requests, false, &out.fields)?; - // Each polygon's row, shape and extent; a feature without polygons meets no point. - let zones: Vec<(usize, MultiPolygon, Extent)> = out + // Each polygon's row, shape, extent and the poles it reaches; a feature without polygons meets + // no point. + let zones: Vec<(usize, MultiPolygon, Extent, (bool, bool))> = out .features .iter() .enumerate() .filter_map(|(row, f)| { - let mp = f.geometry.as_ref().and_then(to_multipolygon)?; + let g = f.geometry.as_ref()?; + let mp = to_multipolygon(g)?; let ext = merc_extent(&mp); - Some((row, mp, ext)) + Some((row, mp, ext, poles(&g.extent()))) }) .collect(); - let index = ExtentIndex::new(zones.iter().map(|(_, _, ext)| *ext)); + let index = ExtentIndex::new(zones.iter().map(|(_, _, ext, _)| *ext)); + let at_poles = AtPoles::new(zones.iter().map(|(.., at)| *at)); let mut counts = vec![0i64; out.features.len()]; // Each polygon's points, kept only for statistics. let mut members: Vec> = if requests.is_empty() { Vec::new() } else { vec![Vec::new(); out.features.len()] }; @@ -3642,10 +3724,13 @@ fn count_points(inp: &dyn Inputs, p: &Json) -> Result { let ext = merc_extent(&pg); let mut met = false; index.candidates_into(&ext, &mut cand_buf); + // A point at a pole touches a polygon reaching it, at whatever longitude (a line or + // polygon passed as points, meeting one only there, has no part in it). + let at = if matches!(g, Geometry::Point(_) | Geometry::MultiPoint(_)) { poles(&g.extent()) } else { (false, false) }; + at_poles.add_to(at, &mut cand_buf); for &zi in &cand_buf { - let Some((row, mp, z_ext)) = zones.get(zi) else { continue }; - if z_ext.intersects(&ext) - && mp.intersects(&pg) + let Some((row, mp, z_ext, z_at)) = zones.get(zi) else { continue }; + if (meet_at_a_pole(at, *z_at) || (z_ext.intersects(&ext) && mp.intersects(&pg))) && let Some(c) = counts.get_mut(*row) { *c += 1; @@ -3745,54 +3830,87 @@ fn lonlat_ring_area(ring: &[Coord]) -> f64 { (ring.iter().zip(ring.iter().cycle().skip(1)).map(|(a, b)| d(a).0 * d(b).1 - d(b).0 * d(a).1).sum::() / 2.0).abs() } -/// The area (square metres, holes taken out) of a lon/lat polygon geometry on a sphere of the -/// Earth's radius, as the map draws it: each edge a straight line in Web Mercator, which is a rhumb -/// line, with longitudes as they are (not folded across the antimeridian) and latitudes clamped to -/// the map's limit, as the overlay takes them. A polygon cut along any straight line on the map -/// splits this area exactly, so the parts' areas add up to the whole's, up to the overlay's -/// rounding of the cut ([`OverlayGrid`]). Zero for any other geometry. +/// The area (square meters, holes taken out) of a lon/lat polygon geometry on a sphere of the +/// Earth's radius, as the overlay takes it: each edge straight in the analysis plane +/// ([`rustgis_geom::lonlat_to_plane`]), a rhumb line within ±85.05° and straight in longitude and +/// latitude beyond (#1033), with longitudes as they are (not folded across the antimeridian). A +/// polygon cut along any straight line in that plane splits this area exactly, so the parts' areas +/// add up to the whole's, up to the overlay's rounding of the cut ([`OverlayGrid`]). Zero for any +/// other geometry. fn map_area_m2(g: &Geometry) -> f64 { let Geometry::Polygon(polys) = g else { return 0.0 }; let r2 = rustgis_geom::EARTH_RADIUS_M * rustgis_geom::EARTH_RADIUS_M; polys.iter().flat_map(|rings| rings.iter().enumerate()).map(|(i, r)| if i == 0 { 1.0 } else { -1.0 } * rhumb_ring_area(r) * r2).sum() } -/// The area of a ring of rhumb lines on the unit sphere, |∮ sin φ dλ|, for [`map_area_m2`]. The -/// integral is the area element cos φ dφ dλ taken round the ring by Green's theorem, the one -/// Chamberlain & Duquette (JPL, 2007) evaluate with edges straight in longitude and latitude -/// ([`rustgis_geom::geodesic_ring_area_m2`]). Along a rhumb line the isometric latitude -/// ψ = ln tan(π/4 + φ/2), the Mercator y over the radius (Snyder, *Map Projections: A Working -/// Manual*, USGS Professional Paper 1395, 1987, on the Mercator projection), changes linearly with -/// the longitude λ, and sin φ = tanh ψ, so an edge's integral is Δλ/Δψ × (ln cosh ψ₂ − ln cosh ψ₁); -/// for a short Δψ, where that difference cancels, its Taylor series about the middle. Longitudes -/// are not folded, as on the map, so the ring's Δλ add up to 0 and no ring winds around a pole: a -/// ring stored around one with a jump across the antimeridian is the flat shape the map draws. +/// The area of a ring on the unit sphere, |∮ sin φ dλ|, each edge straight in the analysis plane, +/// for [`map_area_m2`]. The integral is the area element cos φ dφ dλ taken round the ring by +/// Green's theorem, the one Chamberlain & Duquette (JPL, 2007) evaluate with edges straight in +/// longitude and latitude ([`rustgis_geom::geodesic_ring_area_m2`]). An edge crossing ±85.05° is +/// cut where it crosses, found along it in the plane, where longitude and the plane's y change +/// together. Within ±85.05° a piece is a rhumb line: the isometric latitude ψ = ln tan(π/4 + φ/2), +/// the Mercator y over the radius (Snyder, *Map Projections: A Working Manual*, USGS Professional +/// Paper 1395, 1987, on the Mercator projection), changes linearly with the longitude λ, and +/// sin φ = tanh ψ, so its integral is Δλ/Δψ × (ln cosh ψ₂ − ln cosh ψ₁). Beyond, φ changes +/// linearly with λ, so it is Δλ × (cos φ₁ − cos φ₂)/Δφ. For a short Δψ or Δφ, where those +/// differences cancel, each takes its Taylor series about the middle. Longitudes are not folded, +/// as on the map, so the ring's Δλ add up to 0 and no ring winds around a pole: a ring stored +/// around one with a jump across the antimeridian is the flat shape the map draws. fn rhumb_ring_area(ring: &[Coord]) -> f64 { - let iso = |c: &Coord| c.y.clamp(-rustgis_geom::MAX_LAT, rustgis_geom::MAX_LAT).to_radians().tan().asinh(); + let max_lat = rustgis_geom::MAX_LAT; + let iso = |lat: f64| lat.clamp(-max_lat, max_lat).to_radians().tan().asinh(); let ln_cosh = |x: f64| x.abs() + (-2.0 * x.abs()).exp().ln_1p() - std::f64::consts::LN_2; + // The mean of sin φ along a piece from `p` to `q`, straight in the plane on one side of ±85.05°. + let mean_sin = |p: Coord, q: Coord| { + if p.y.abs().max(q.y.abs()) <= max_lat { + let (p1, p2) = (iso(p.y), iso(q.y)); + let dpsi = p2 - p1; + if dpsi.abs() < 1e-4 { + // tanh over the piece, averaged: tanh m + Δψ²/24 × tanh″ m, with tanh″ = −2 tanh sech². + let t = ((p1 + p2) / 2.0).tanh(); + t - dpsi * dpsi / 12.0 * t * (1.0 - t * t) + } else { + (ln_cosh(p2) - ln_cosh(p1)) / dpsi + } + } else { + let (f1, f2) = (p.y.to_radians(), q.y.to_radians()); + let dphi = f2 - f1; + if dphi.abs() < 1e-4 { + // sin over the piece, averaged: sin m × (1 − Δφ²/24). + ((f1 + f2) / 2.0).sin() * (1.0 - dphi * dphi / 24.0) + } else { + (f1.cos() - f2.cos()) / dphi + } + } + }; + let edge_y = |lat: f64| rustgis_geom::lonlat_to_plane(Coord::new(0.0, lat)).y; + let (north, south) = (edge_y(max_lat), edge_y(-max_lat)); let mut total = 0.0; for (a, b) in ring.iter().zip(ring.iter().cycle().skip(1)) { - let dlon = b.x - a.x; - let (p1, p2) = (iso(a), iso(b)); - let dpsi = p2 - p1; - let mean_sin = if dpsi.abs() < 1e-4 { - // tanh over the edge, averaged: tanh m + Δψ²/24 × tanh″ m, with tanh″ = −2 tanh sech². - let t = ((p1 + p2) / 2.0).tanh(); - t - dpsi * dpsi / 12.0 * t * (1.0 - t * t) - } else { - (ln_cosh(p2) - ln_cosh(p1)) / dpsi - }; - total += dlon.to_radians() * mean_sin; + let (a, b) = (Coord::new(a.x, a.y.clamp(-90.0, 90.0)), Coord::new(b.x, b.y.clamp(-90.0, 90.0))); + let (ya, yb) = (rustgis_geom::lonlat_to_plane(a).y, rustgis_geom::lonlat_to_plane(b).y); + // Where the edge crosses 85.05° north or south, in order along it. + let mut cuts: Vec<(f64, f64)> = [(north, max_lat), (south, -max_lat)] + .into_iter() + .map(|(y, lat)| ((y - ya) / (yb - ya), lat)) + .filter(|(t, _)| *t > 0.0 && *t < 1.0) + .collect(); + cuts.sort_by(|p, q| p.0.total_cmp(&q.0)); + let mut from = a; + for to in cuts.into_iter().map(|(t, lat)| Coord::new(a.x + t * (b.x - a.x), lat)).chain(std::iter::once(b)) { + total += (to.x - from.x).to_radians() * mean_sin(from, to); + from = to; + } } total.abs() } -/// The overlay's grid for a feature cut by masks within a Web Mercator extent. The overlay keeps 29 +/// The overlay's grid for a feature cut by masks within an analysis-plane extent. The overlay keeps 29 /// integer bits about the centre of its operands' extent (see [`clip`]), a grid step of the extent /// over 2^30, so it moves each vertex by less than a step: the area by less than the perimeter times /// a step, and a sliver it leaves along an edge is no wider than a step. struct OverlayGrid { - /// The grid step, in Web Mercator metres, taken over the feature and all its masks. + /// The grid step, in the analysis plane's meters, taken over the feature and all its masks. step: f64, /// How far rounding can move the feature's area, relative to it: twice the perimeter times the /// step over the area, from 1e-9 to 1. A part closer than this to the whole may be the whole. @@ -3826,7 +3944,7 @@ impl OverlayGrid { } } -/// The planar area and perimeter (holes included) of polygons in Web Mercator metres. +/// The planar area and perimeter (holes included) of polygons in the analysis plane's meters. fn merc_area_perimeter(mp: &MultiPolygon) -> (f64, f64) { use geo::{Area, Euclidean, Length}; let perimeter = mp.0.iter().flat_map(|p| std::iter::once(p.exterior()).chain(p.interiors())).map(|r| Euclidean.length(r)).sum(); @@ -4316,11 +4434,12 @@ struct TabulateZone { key: Vec, /// The lon/lat extent of its features, for the index. ext: Extent, - /// Its features' polygons in Web Mercator, unioned when there are several, as [`mask_of`] does. + /// Its features' polygons in the analysis plane, unioned when there are several, as + /// [`mask_of`] does. mask: MultiPolygon, - /// The Web Mercator extent of its mask, for [`OverlayGrid`]. + /// The analysis-plane extent of its mask, for [`OverlayGrid`]. merc_ext: Extent, - /// Its area in square metres as the map draws it ([`map_area_m2`]), overlaps between its + /// Its area in square meters as the overlay takes it ([`map_area_m2`]), overlaps between its /// features counted once. map_area_m2: f64, } @@ -4381,8 +4500,7 @@ fn tabulate_intersection(inp: &dyn Inputs, p: &Json) -> Result 0.0))); let shapeless = |ds: &Dataset, rs: &[usize]| { rs.iter().filter(|&&r| ds.features.get(r).and_then(|f| f.geometry.as_ref()).is_none_or(|g| !has_shape(g))).count() @@ -4390,15 +4508,6 @@ fn tabulate_intersection(inp: &dyn Inputs, p: &Json) -> Result rustgis_geom::MAX_LAT || e.min_y < -rustgis_geom::MAX_LAT) - }; - rs.iter().filter(|&&r| ds.features.get(r).and_then(|f| f.geometry.as_ref()).is_some_and(beyond)).count() - }; - let (zones_beyond, classes_beyond) = (beyond_limit(zones, &zone_rows), beyond_limit(classes, &class_rows)); // Features with a ring that jumps across the antimeridian unsplit: the overlay, like the map, // takes such a ring the long way round, while `area_km2` takes its ground shape. let unsplit = |ds: &Dataset, rs: &[usize]| { @@ -4411,6 +4520,20 @@ fn tabulate_intersection(inp: &dyn Inputs, p: &Json) -> Result max || b.y.abs() > max) && a.x != b.x) + })) + }; + rs.iter().filter(|&&r| ds.features.get(r).and_then(|f| f.geometry.as_ref()).is_some_and(across)).count() + }; + let (zones_polar, classes_polar) = (polar(zones, &zone_rows), polar(classes, &class_rows)); let zone_list: Vec = crate::stats::group_rows(zones, &zone_rows, &zone_keys) .into_iter() .map(|group| { @@ -4648,29 +4771,28 @@ fn tabulate_intersection(inp: &dyn Inputs, p: &Json) -> Result = [(zones_beyond, "zone"), (classes_beyond, "class")] + let unsplit: Vec = [(zones_unsplit, "zone"), (classes_unsplit, "class")] .into_iter() .filter(|(n, _)| *n > 0) .map(|(n, layer)| features(n, layer)) .collect(); - if !beyond.is_empty() { - let verb = if zones_beyond + classes_beyond == 1 { "reaches" } else { "reach" }; + if !unsplit.is_empty() { + let verb = if zones_unsplit + classes_unsplit == 1 { "crosses" } else { "cross" }; messages.push(format!( - "{} {verb} beyond {:.2}° latitude, where the map ends: parts there are taken as lying along it, with no area", - beyond.join(" and "), - rustgis_geom::MAX_LAT + "{} {verb} the antimeridian without being split there: the map draws them the long way round, and their percentages and shares follow the map, while area_km2 follows the ground", + unsplit.join(" and ") )); } - let unsplit: Vec = [(zones_unsplit, "zone"), (classes_unsplit, "class")] + let polar: Vec = [(zones_polar, "zone"), (classes_polar, "class")] .into_iter() .filter(|(n, _)| *n > 0) .map(|(n, layer)| features(n, layer)) .collect(); - if !unsplit.is_empty() { - let verb = if zones_unsplit + classes_unsplit == 1 { "crosses" } else { "cross" }; + if !polar.is_empty() { + let verb = if zones_polar + classes_polar == 1 { "has" } else { "have" }; messages.push(format!( - "{} {verb} the antimeridian without being split there: the map draws them the long way round, and their percentages and shares follow the map, while area_km2 follows the ground", - unsplit.join(" and ") + "{} {verb} an edge beyond 85.05° running east or west: the analysis plane takes it straight in longitude and latitude beyond 85.05°, and their percentages and shares follow the plane, while area_km2 takes the shortest way on the ground, which near a pole can differ widely", + polar.join(" and ") )); } match unresolved { @@ -4777,14 +4899,85 @@ impl JoinMatch { /// The least distance Closest's first search reaches, in metres. const MIN_REACH_M: f64 = 10.0; +/// The poles a feature with this extent reaches: north, then south. +fn poles(e: &Extent) -> (bool, bool) { + reaches_poles(e.min_y, e.max_y) +} + +/// The features of a list reaching the North Pole and the South Pole, by their place in it: a +/// feature reaching a pole meets them there, however far apart their extents' longitudes are, so +/// they're candidates for it beside those an index finds by extent (#1033). +#[derive(Default)] +struct AtPoles { + north: Vec, + south: Vec, +} + +impl AtPoles { + /// From the poles each feature reaches (north, south), in list order. + fn new(reached: impl IntoIterator) -> Self { + let mut at = Self::default(); + for (i, (north, south)) in reached.into_iter().enumerate() { + if north { + at.north.push(i); + } + if south { + at.south.push(i); + } + } + at + } + + /// Adds to `found` those reaching the poles `at` (north, south) reaches, then sorts it, each + /// once; leaves it as it is when there are none to add. + fn add_to(&self, (north, south): (bool, bool), found: &mut Vec) { + let before = found.len(); + if north { + found.extend_from_slice(&self.north); + } + if south { + found.extend_from_slice(&self.south); + } + if found.len() != before { + found.sort_unstable(); + found.dedup(); + } + } + + /// The bytes its lists hold. + fn storage_bytes(&self) -> usize { + self.north.capacity().saturating_add(self.south.capacity()).saturating_mul(std::mem::size_of::()) + } + + /// Whether any feature reaches a pole `at` (north, south) reaches. + fn any(&self, (north, south): (bool, bool)) -> bool { + (north && !self.north.is_empty()) || (south && !self.south.is_empty()) + } +} + +/// The poles polygons in the analysis plane reach, from their box there: north, then south. +fn plane_poles(mp: &MultiPolygon) -> (bool, bool) { + use geo::BoundingRect; + let lat = |y: f64| rustgis_geom::plane_to_lonlat(Coord::new(0.0, y)).y; + mp.bounding_rect().map_or((false, false), |r| reaches_poles(lat(r.min().y), lat(r.max().y))) +} + +/// Whether two features meet at a pole, given the poles each reaches. The analysis plane, like +/// longitude and latitude, draws a pole as a line, so (0°, −90°) and (90°, −90°) are apart there +/// though they are one place on the ground (#1033). (Distances already meet there: they're +/// measured on the ground.) +fn meet_at_a_pole((a_north, a_south): (bool, bool), (b_north, b_south): (bool, bool)) -> bool { + (a_north && b_north) || (a_south && b_south) +} + /// The ground distance between two features in metres: zero when they intersect. fn ground_distance(t: &Geometry, t_geo: &geo::Geometry, j: &Geometry, j_geo: &geo::Geometry) -> Option { if t_geo.intersects(j_geo) { Some(0.0) } else { crate::join::distance_m(t, j) } } /// How much of a target lies in a join polygon: the ground area of a polygon's overlap, the ground -/// length of a line's. The geometries are in Web Mercator, polygons oriented for the overlay; the -/// overlap found is measured on the ground, in longitude/latitude. +/// length of a line's. The geometries are in the analysis plane, polygons oriented for the overlay; +/// the overlap found is measured on the ground, in longitude/latitude. fn overlap(t_geo: &geo::Geometry, j_geo: &geo::Geometry) -> f64 { match (t_geo, j_geo) { (geo::Geometry::MultiPolygon(t), geo::Geometry::MultiPolygon(j)) => { @@ -4797,7 +4990,7 @@ fn overlap(t_geo: &geo::Geometry, j_geo: &geo::Geometry) -> f64 { } } -/// A geometry in Web Mercator, polygons oriented as the overlay needs. +/// A geometry in the analysis plane, polygons oriented as the overlay needs. fn oriented(g: &Geometry) -> geo::Geometry { match to_geo(g) { geo::Geometry::MultiPolygon(m) => geo::Geometry::MultiPolygon(m.orient(Direction::Default)), @@ -4954,8 +5147,8 @@ impl ExtentIndex { } } -/// A feature to join or measure to: its row in its layer, its extent, and its geometry in Web -/// Mercator (`geo`) and in lon/lat (`geom`). +/// A feature to join or measure to: its row in its layer, its extent, and its geometry in the +/// analysis plane (`geo`) and in lon/lat (`geom`). struct Located<'a> { row: usize, ext: Extent, @@ -4971,11 +5164,11 @@ struct Nearby<'a, 'b> { } impl Nearby<'_, '_> { - /// The features within `reach` metres of a target (`t`, with its extent, and its geometry in Web - /// Mercator), in order, with their distances. With `nearest`, only the nearest is wanted: a - /// feature that can't be as near as the nearest so far isn't measured, so some farther ones - /// may be missing. The feature in row `exclude` is left out (the target itself, when both are - /// one layer). `measured` keeps each feature's distance, measured once per target. + /// The features within `reach` meters of a target (`t`, with its extent, and its geometry in + /// the analysis plane), in order, with their distances. With `nearest`, only the nearest is + /// wanted: a feature that can't be as near as the nearest so far isn't measured, so some + /// farther ones may be missing. The feature in row `exclude` is left out (the target itself, + /// when both are one layer). `measured` keeps each feature's distance, once per target. fn within( &self, t: &Geometry, @@ -5056,11 +5249,11 @@ impl Nearby<'_, '_> { /// A point (lon/lat) two intersecting features share: one of the points of either on the other; /// else a point of their overlap, when one is polygons; else a point where their edges meet. The -/// geometries come with their Web Mercator forms, where the test and the overlay work. +/// geometries come with their forms in the analysis plane, where the test and the overlay work. fn shared_point(a: &Geometry, a_geo: &geo::Geometry, b: &Geometry, b_geo: &geo::Geometry) -> Option { use geo::LinesIter; - let back = |c: geo::Coord| rustgis_geom::merc_to_lonlat(Coord::new(c.x, c.y)); - // Points answer with their own coordinates, not ones back from Web Mercator. + let back = |c: geo::Coord| rustgis_geom::plane_to_lonlat(Coord::new(c.x, c.y)); + // Points answer with their own coordinates, not ones back from the plane. let point_on = |pts: &Geometry, other: &geo::Geometry| -> Option> { let pts: &[Coord] = match pts { Geometry::Point(c) => std::slice::from_ref(c), @@ -5068,7 +5261,7 @@ fn shared_point(a: &Geometry, a_geo: &geo::Geometry, b: &Geometry, b_geo: & _ => return None, }; Some(pts.iter().copied().find(|c| { - let m = rustgis_geom::lonlat_to_merc(*c); + let m = rustgis_geom::lonlat_to_plane(*c); other.intersects(&geo::Point::new(m.x, m.y)) })) }; @@ -5305,6 +5498,12 @@ fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result { let index = ExtentIndex::new(join_geoms.iter().map(|j| j.ext)); let nearby = Nearby { feats: &join_geoms, index: &index }; + // For Intersects and Center In, the join features reaching each pole. + let joins_at_poles = if matches!(mode, JoinMatch::Intersects | JoinMatch::CentreIn) { + AtPoles::new(join_geoms.iter().map(|j| poles(&j.ext))) + } else { + AtPoles::default() + }; let mut total_matches = 0usize; let mut written_features = 0usize; @@ -5323,11 +5522,13 @@ fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result { match mode { JoinMatch::Intersects | JoinMatch::Contains | JoinMatch::Within => { index.candidates_into(&t_ext, &mut cand_buf); + let t_poles = poles(&t_ext); + joins_at_poles.add_to(t_poles, &mut cand_buf); for j in cand_buf.iter().copied().filter_map(at) { let hit = match mode { JoinMatch::Contains => extent_contains(&t_ext, &j.ext) && Contains::contains(&t_geo, &j.geo), JoinMatch::Within => extent_contains(&j.ext, &t_ext) && Contains::contains(&j.geo, &t_geo), - _ => t_ext.intersects(&j.ext) && t_geo.intersects(&j.geo), + _ => meet_at_a_pole(t_poles, poles(&j.ext)) || (t_ext.intersects(&j.ext) && t_geo.intersects(&j.geo)), }; if hit { matches.push((j.row, None)); @@ -5338,8 +5539,10 @@ fn spatial_join(inp: &dyn Inputs, p: &Json) -> Result { if let Some(pt) = interior_point(tg) { let (p_ext, p_geo) = (pt.extent(), to_geo(&pt)); index.candidates_into(&p_ext, &mut cand_buf); + let p_poles = poles(&p_ext); + joins_at_poles.add_to(p_poles, &mut cand_buf); for j in cand_buf.iter().copied().filter_map(at) { - if j.ext.intersects(&p_ext) && j.geo.intersects(&p_geo) { + if meet_at_a_pole(p_poles, poles(&j.ext)) || (j.ext.intersects(&p_ext) && j.geo.intersects(&p_geo)) { matches.push((j.row, None)); } } @@ -6023,6 +6226,28 @@ mod tests { assert_eq!(previous, 12, "three closed triangle rings at a very large tolerance"); } + #[test] + fn simplify_says_how_many_features_reach_beyond_85_degrees() { + // A box from 86°N to 87°N reaches beyond 85.05°, where less is taken out east–west than + // asked: the tool says so. One from 80°N to 84°N, or a tolerance of 0, says nothing. + let said = |y0: f64, y1: f64, tolerance: f64| { + let ring = [(0.0, y0), (1.0, y0), (1.0, y1), (0.0, y1), (0.0, y0)].map(|(x, y)| Coord::new(x, y)).to_vec(); + let mut ds = Dataset::new("d", GeometryType::Polygon, vec![]); + ds.push(Some(Geometry::Polygon(vec![vec![ring]])), vec![]); + let m = Map(HashMap::from([("d".to_string(), ds)])); + let res = crate::run("simplify", &m, &serde_json::json!({"input": "d", "tolerance": tolerance})).unwrap(); + res.messages.iter().any(|m| m.contains("beyond 85.05°")) + }; + assert!(said(86.0, 87.0, 1.0)); + assert!(!said(80.0, 84.0, 1.0) && !said(86.0, 87.0, 0.0)); + // A point beyond 85.05° isn't simplified, so isn't counted. + let mut ds = Dataset::new("d", GeometryType::Point, vec![]); + ds.push(Some(Geometry::Point(Coord::new(0.0, 88.0))), vec![]); + let m = Map(HashMap::from([("d".to_string(), ds)])); + let res = crate::run("simplify", &m, &serde_json::json!({"input": "d", "tolerance": 1.0})).unwrap(); + assert!(res.messages.iter().all(|m| !m.contains("beyond 85.05°")), "{:?}", res.messages); + } + #[test] fn simplify_keeps_points_lines_and_degenerate_rings() { use geo::Area; @@ -7608,14 +7833,14 @@ mod tests { } /// Points along each edge of each ring at the fractions `ts`, as the tools draw edges: straight - /// in Web Mercator. + /// in the analysis plane (Web Mercator within ±85.05°, straight in lon/lat beyond). fn along_edges(polys: &[Vec>], ts: &[f64]) -> Vec { - use rustgis_geom::{lonlat_to_merc, merc_to_lonlat}; + use rustgis_geom::{lonlat_to_plane, plane_to_lonlat}; polys .iter() .flatten() - .flat_map(|ring| ring.windows(2).map(|w| (lonlat_to_merc(w[0]), lonlat_to_merc(w[1]))).collect::>()) - .flat_map(|(a, b)| ts.iter().map(move |t| merc_to_lonlat(Coord::new(a.x + t * (b.x - a.x), a.y + t * (b.y - a.y))))) + .flat_map(|ring| ring.windows(2).map(|w| (lonlat_to_plane(w[0]), lonlat_to_plane(w[1]))).collect::>()) + .flat_map(|(a, b)| ts.iter().map(move |t| plane_to_lonlat(Coord::new(a.x + t * (b.x - a.x), a.y + t * (b.y - a.y))))) .collect() } @@ -7664,7 +7889,7 @@ mod tests { /// The outline keeps its distance between its vertices too, not only at them (#1159): edges /// drawn straight in the projection came back to lon/lat as straight lines there, so the long /// sides of a 2 km buffer of a 96 km east–west line sat up to 190 m off. Every vertex, and every - /// edge's quarter points and midpoint as the tools draw edges (straight in Web Mercator), is now + /// edge's quarter points and midpoint as the tools draw edges (straight in the analysis plane), is now /// within 0.012% of the distance (or 1 cm) along a line's sides, and no further inside an arc /// than 0.031% of the distance in the projection plane (#1221). The old 32-edge circles /// and 15-edge caps cut about 0.5% inside and fail these bounds. @@ -7716,11 +7941,13 @@ mod tests { } } - // Near a pole: a circle reaching past Web Mercator's 85.05° keeps its distance up to 85° - // (an edge across 85.05° was left whole, its visible part 80 km out). Within the last - // twentieth of a degree, where Mercator stretches its edges most, it can be 0.03% out. + // Near a pole: a circle reaching past 85.05°, and over the pole, keeps its distance all the + // way round (an edge across 85.05° was once left whole, its visible part 80 km out; edges + // beyond it were left whole until #1033). The edges closing the cap, up ±180° to the pole + // and along it, aren't on the circle. let north = Coord::new(10.0, 83.0); - for c in points(&buffered(Geometry::Point(north), 1000.0)).into_iter().filter(|c| c.y.abs() < 85.0) { + let closing = |c: &Coord| c.y.abs() >= 90.0 || c.x.abs() >= 180.0 - 1e-9; + for c in points(&buffered(Geometry::Point(north), 1000.0)).into_iter().filter(|c| !closing(c)) { let d = on_wgs84_m(north, c); assert!(d < 1_000_000.0 * (1.0 + 1.2e-4) && d > 1_000_000.0 - arc_error(1_000_000.0) - 100.0, "{c:?} is {d} m from {north:?}"); } @@ -8050,6 +8277,48 @@ mod tests { assert_eq!(with_budget("Distance must be a number".into(), 512 << 20), "Distance must be a number"); } + #[test] + fn kernel_density_skips_points_its_system_cant_hold() { + // In Web Mercator, points past 85.05° would land on the map's edge: they're left out and + // counted, not piled up on the raster's top row (#1033), a multipoint's members each, and + // one a hair past the limit (which the projection's round trip alone would let through). + let mut pts = Dataset::new("pts", GeometryType::Point, vec![]); + let hair = Coord::new(0.0, rustgis_geom::MAX_LAT + 5e-5); + for c in [Coord::new(0.0, 87.0), Coord::new(1.0, 89.0), Coord::new(0.0, 60.0), Coord::new(0.01, 60.0), hair] { + pts.push(Some(Geometry::Point(c)), vec![]); + } + pts.push(Some(Geometry::MultiPoint(vec![Coord::new(0.0, 60.02), Coord::new(0.0, 87.0)])), vec![]); + let m = Map(HashMap::from([("pts".to_string(), pts)])); + let params = serde_json::json!({"input": "pts", "crs": "EPSG:3857", "radius": 5000.0, "pixelSize": 1000.0}); + let res = crate::run("kernelDensity", &m, ¶ms).unwrap(); + assert!(res.messages.iter().any(|m| m.contains("from 3 points (bandwidth")), "{:?}", res.messages); + assert!(res.messages.iter().any(|m| m.starts_with("4 point(s) beyond what") && m.contains("85.05°")), "{:?}", res.messages); + // With none it can hold, the error says why. + let mut far = Dataset::new("far", GeometryType::Point, vec![]); + far.push(Some(Geometry::Point(Coord::new(0.0, 87.0))), vec![]); + let m = Map(HashMap::from([("far".to_string(), far)])); + let params = serde_json::json!({"input": "far", "crs": "EPSG:3857", "radius": 5000.0, "pixelSize": 1000.0}); + let err = crate::run("kernelDensity", &m, ¶ms).err().unwrap_or_default(); + assert!(err.contains("1 point(s) beyond what"), "{err}"); + // Points it can hold say nothing of the kind. + let mut near = Dataset::new("near", GeometryType::Point, vec![]); + near.push(Some(Geometry::Point(Coord::new(0.0, 60.0))), vec![]); + let m = Map(HashMap::from([("near".to_string(), near)])); + let params = serde_json::json!({"input": "near", "crs": "EPSG:3857", "radius": 5000.0, "pixelSize": 1000.0}); + assert!(crate::run("kernelDensity", &m, ¶ms).unwrap().messages.iter().all(|m| !m.contains("beyond"))); + // In another system the note names it, not Web Mercator's limit: UTM zone 18N can't hold a + // point a quarter of the way round the equator from its central meridian (75°W). + let mut off = Dataset::new("off", GeometryType::Point, vec![]); + for c in [Coord::new(-75.0, 40.0), Coord::new(-74.99, 40.0), Coord::new(15.0, 0.0)] { + off.push(Some(Geometry::Point(c)), vec![]); + } + let m = Map(HashMap::from([("off".to_string(), off)])); + let params = serde_json::json!({"input": "off", "crs": "EPSG:32618", "radius": 5000.0, "pixelSize": 1000.0}); + let res = crate::run("kernelDensity", &m, ¶ms).unwrap(); + let note = res.messages.iter().find(|m| m.contains("beyond what")).cloned().unwrap_or_default(); + assert!(note.starts_with("1 point(s) beyond what") && !note.contains("85.05°"), "{:?}", res.messages); + } + #[test] fn kernel_density_analysis_tool_tests() { // Create point dataset with UTM projection @@ -8585,6 +8854,30 @@ mod tests { assert!(around_pole.abs() < 1.0, "{around_pole} m²"); } + /// An edge crossing 85.05° and running east or west is cut where it crosses: a rhumb line on one + /// side, straight in longitude and latitude on the other. The map measure of a triangle with two + /// such edges agrees with ∮ sin φ dλ summed in a million steps along each edge as the plane draws + /// it, independently of the cut. + #[test] + fn map_area_cuts_edges_across_85_degrees() { + let ring: Vec = [(0.0, 80.0), (20.0, 88.0), (-15.0, 87.0), (0.0, 80.0)].iter().map(|&(x, y)| Coord::new(x, y)).collect(); + let measured = map_area_m2(&Geometry::Polygon(vec![vec![ring.clone()]])); + let steps = 1_000_000; + let mut sum = 0.0; + for w in ring.windows(2) { + let (a, b) = (rustgis_geom::lonlat_to_plane(w[0]), rustgis_geom::lonlat_to_plane(w[1])); + let at = |t: f64| rustgis_geom::plane_to_lonlat(Coord::new(a.x + (b.x - a.x) * t, a.y + (b.y - a.y) * t)); + let mut last = at(0.0); + for k in 1..=steps { + let (next, mid) = (at(f64::from(k) / f64::from(steps)), at((f64::from(k) - 0.5) / f64::from(steps))); + sum += mid.y.to_radians().sin() * (next.x - last.x).to_radians(); + last = next; + } + } + let expected = sum.abs() * rustgis_geom::EARTH_RADIUS_M * rustgis_geom::EARTH_RADIUS_M; + assert!((measured / expected - 1.0).abs() < 1e-9, "{measured} m², not {expected}"); + } + /// Parcels and zones for Intersect: `p1` is a 2° square that zone `z1` covers the right half of; /// `z2` overlaps `z1` but only touches `p1` along its right edge, and `p2` along its left. /// Both layers have an `id` field. @@ -9038,8 +9331,7 @@ mod tests { let gg = to_geo(&buffered); let ext = merc_extent(&gg); let n = 32; - let before_index = - n * (selector_heap_bytes(&gg) + std::mem::size_of::<(geo::Geometry, Extent)>() + std::mem::size_of::()); + let before_index = n * (selector_heap_bytes(&gg) + std::mem::size_of::() + std::mem::size_of::()); let index_bytes = ExtentIndex::new((0..n).map(|_| ext)).storage_bytes(); assert!(index_bytes > 0); let mut selectors = Dataset::new("selectors", GeometryType::Point, vec![]); diff --git a/crates/analysis/src/topology.rs b/crates/analysis/src/topology.rs index 028b0c319..761318a1a 100644 --- a/crates/analysis/src/topology.rs +++ b/crates/analysis/src/topology.rs @@ -1,9 +1,9 @@ //! Editing geometry: merging features, and topology checks (overlaps, gaps, self-intersections, -//! dangles, coverage). Computed in Web Mercator metres. +//! dangles, coverage). Computed in the analysis plane ([`rustgis_geom::lonlat_to_plane`]). use geo::{Area, BooleanOps, BoundingRect, Contains, Intersects, Orient, orient::Direction}; use rustgis_data::{Dataset, Geometry}; -use rustgis_geom::{Coord, lonlat_to_merc, merc_to_lonlat}; +use rustgis_geom::{Coord, MAX_LAT, WGS84_POLAR_RADIUS_M, lonlat_to_plane, plane_to_lonlat}; use crate::convert::{from_multipolygon, to_multipolygon}; @@ -14,19 +14,22 @@ pub struct TopologyError { pub at: Coord, } -/// A check's errors, and whether it finished (a check stops at an error or work limit). +/// A check's errors, whether it finished (a check stops at an error or work limit) with every part +/// measured as closely as it's meant to be, and how many polygons it had to leave out (see the +/// check for why). #[derive(Clone, Debug, Default, PartialEq)] pub struct Checked { pub errors: Vec, pub complete: bool, + pub skipped: usize, } impl Checked { fn done(errors: Vec) -> Checked { - Checked { errors, complete: true } + Checked { errors, complete: true, skipped: 0 } } fn partial(errors: Vec) -> Checked { - Checked { errors, complete: false } + Checked { errors, complete: false, skipped: 0 } } } @@ -85,7 +88,7 @@ fn polygons(ds: &Dataset, rows: &[usize]) -> Vec<(usize, geo::MultiPolygon, fn centre(mp: &geo::MultiPolygon) -> Coord { use geo::InteriorPoint; - mp.interior_point().map(|p| merc_to_lonlat(Coord::new(p.x(), p.y()))).unwrap_or_default() + mp.interior_point().map(|p| plane_to_lonlat(Coord::new(p.x(), p.y()))).unwrap_or_default() } fn rects_touch(a: &geo::Rect, b: &geo::Rect) -> bool { @@ -215,9 +218,9 @@ pub fn self_intersections(ds: &Dataset, rows: &[usize]) -> Checked { _ => Vec::new(), }; for (p, closed) in paths { - let m: Vec = p.iter().map(|c| lonlat_to_merc(*c)).collect(); + let m: Vec = p.iter().map(|c| lonlat_to_plane(*c)).collect(); if let Some(c) = self_crossing(&m, closed) { - out.push(TopologyError { rows: vec![*r], at: merc_to_lonlat(c) }); + out.push(TopologyError { rows: vec![*r], at: plane_to_lonlat(c) }); break; } } @@ -229,25 +232,27 @@ pub fn self_intersections(ds: &Dataset, rows: &[usize]) -> Checked { } /// Line ends that touch no other line (within `tol_m` metres). Line parts are indexed on a grid, -/// so each end is measured only against parts near it. +/// so each end is measured only against parts near it; an end beyond 85.05° in a frame about its +/// pole (`PolarEnds`). pub fn dangles(ds: &Dataset, rows: &[usize], tol_m: f64) -> Checked { let lines: Vec<(usize, Vec>)> = rows .iter() .filter_map(|r| match ds.features.get(*r)?.geometry.as_ref()? { - Geometry::Line(parts) => Some((*r, parts.iter().map(|p| p.iter().map(|c| lonlat_to_merc(*c)).collect()).collect())), + Geometry::Line(parts) => Some((*r, parts.iter().map(|p| p.iter().map(|c| lonlat_to_plane(*c)).collect()).collect())), _ => None, }) .collect(); - // Tolerance in Mercator metres grows with latitude; index with the largest it can be here. - let max_tol = tol_m / 85f64.to_radians().cos(); + // Within ±85.05°, tolerance in the plane's meters grows with latitude; index with the largest + // it can be there. + let max_tol = tol_m / MAX_LAT.to_radians().cos(); let mut ext = (f64::MAX, f64::MAX, f64::MIN, f64::MIN); for c in lines.iter().flat_map(|(_, ps)| ps.iter().flatten()) { ext = (ext.0.min(c.x), ext.1.min(c.y), ext.2.max(c.x), ext.3.max(c.y)); } let span = (ext.2 - ext.0).max(ext.3 - ext.1).max(1.0); - let cell = (span / 512.0).max(max_tol.min(span)).max(1e-6); - let key = |x: f64, y: f64| (((x - ext.0) / cell).floor() as i64, ((y - ext.1) / cell).floor() as i64); - // Each part's segments go into the cells their (widened) boxes cover. + let cells = Cells::new((ext.0, ext.1), (span / 512.0).max(max_tol.min(span)).max(1e-6), (ext.2, ext.3)); + // Each part's segments go into the cells their (widened) boxes cover, within the lines' extent: + // only those are ever looked in, as every end lies in it. let mut grid: std::collections::HashMap<(i64, i64), Vec<(usize, usize)>> = std::collections::HashMap::new(); let mut work = 0usize; crate::control::begin(lines.len()); @@ -259,10 +264,10 @@ pub fn dangles(ds: &Dataset, rows: &[usize], tol_m: f64) -> Checked { for (pi, part) in parts.iter().enumerate() { for w in part.windows(2) { let (Some(a), Some(b)) = (w.first(), w.get(1)) else { continue }; - let (x0, y0) = key(a.x.min(b.x) - max_tol, a.y.min(b.y) - max_tol); - let (x1, y1) = key(a.x.max(b.x) + max_tol, a.y.max(b.y) + max_tol); - for gx in x0..=x1 { - for gy in y0..=y1 { + let (x0, y0) = cells.key(a.x.min(b.x) - max_tol, a.y.min(b.y) - max_tol); + let (x1, y1) = cells.key(a.x.max(b.x) + max_tol, a.y.max(b.y) + max_tol); + for gx in x0.max(0)..=x1.min(cells.last.0) { + for gy in y0.max(0)..=y1.min(cells.last.1) { work += 1; if work > MAX_WORK { return Checked::partial(Vec::new()); @@ -278,6 +283,8 @@ pub fn dangles(ds: &Dataset, rows: &[usize], tol_m: f64) -> Checked { } let mut out = Vec::new(); crate::control::begin(lines.len()); + // Made at the first end beyond 85.05°. + let mut polar: Option = None; for (li, (r, parts)) in lines.iter().enumerate() { if crate::control::cancelled() { return Checked::partial(out); @@ -288,27 +295,37 @@ pub fn dangles(ds: &Dataset, rows: &[usize], tol_m: f64) -> Checked { if part.len() < 2 || part.first() == part.last() { continue; } - for end in [part.first(), part.last()].into_iter().flatten() { - let tol = tol_m / merc_to_lonlat(*end).y.to_radians().cos().max(0.01); - let near = grid.get(&key(end.x, end.y)).map(Vec::as_slice).unwrap_or(&[]); - // The end's own part counts too, past the stretch next to the end (a line curling - // back onto itself, like a 6, doesn't dangle there). - let at_start = part.first() == Some(end); - let own_rest = rest_of_part(part, at_start, 2.0 * tol); - let touches = !own_rest.is_empty() && rustgis_data::geometry::seg_dist2(&own_rest, *end) <= tol * tol - || near.iter().any(|&(oli, opi)| { - work += 1; - (oli != li || opi != pi) - && lines - .get(oli) - .and_then(|(_, ps)| ps.get(opi)) - .is_some_and(|op| op.len() >= 2 && rustgis_data::geometry::seg_dist2(op, *end) <= tol * tol) - }); + for (end, at_start) in [(part.first(), true), (part.last(), false)] { + let Some(end) = end else { continue }; + let at = plane_to_lonlat(*end); + let lat = at.y; + let touches = if lat.abs() > MAX_LAT { + let polar = polar.get_or_insert_with(|| PolarEnds::new(&lines, &grid, cells)); + match polar.meets((li, pi, at_start), (*end, at), tol_m, &mut work) { + Ok(touches) => touches, + Err(WorkLimit) => return Checked::partial(out), + } + } else { + let tol = tol_m / lat.to_radians().cos().max(0.01); + let near = grid.get(&cells.key(end.x, end.y)).map(Vec::as_slice).unwrap_or(&[]); + // The end's own part counts too, past the stretch next to the end (a line + // curling back onto itself, like a 6, doesn't dangle there). + let own_rest = rest_of_part(part, at_start, 2.0 * tol); + !own_rest.is_empty() && rustgis_data::geometry::seg_dist2(&own_rest, *end) <= tol * tol + || near.iter().any(|&(oli, opi)| { + work += 1; + (oli != li || opi != pi) + && lines + .get(oli) + .and_then(|(_, ps)| ps.get(opi)) + .is_some_and(|op| op.len() >= 2 && rustgis_data::geometry::seg_dist2(op, *end) <= tol * tol) + }) + }; if work > MAX_WORK { return Checked::partial(out); } if !touches { - out.push(TopologyError { rows: vec![*r], at: merc_to_lonlat(*end) }); + out.push(TopologyError { rows: vec![*r], at }); } } } @@ -316,7 +333,457 @@ pub fn dangles(ds: &Dataset, rows: &[usize], tol_m: f64) -> Checked { return Checked::partial(out); } } - Checked::done(out) + // Every end was checked, but one measured coarsely may have been found wrongly. + if polar.is_some_and(|p| p.coarse) { Checked::partial(out) } else { Checked::done(out) } +} + +/// The most a piece of an edge may stray from the curve it stands for in [`edge_from`], as a share +/// of the tolerance; at least [`POLAR_STRAY_MIN_M`]. +const POLAR_STRAY_SHARE: f64 = 1e-3; +/// The least that allowance, in meters. +const POLAR_STRAY_MIN_M: f64 = 1e-3; +/// The most pieces [`edge_from`] cuts one stretch into, and [`walk_edge`] an edge. Only one winding +/// round a pole many times needs more (a hostile edge thousands of degrees long): it's measured in +/// this many, which stray farther, and the check says it's incomplete. +const POLAR_PIECES_MAX: f64 = 65_536.0; +/// The longest piece, in radians of longitude, that [`walk_edge`] measures an edge in: its chords +/// fall short of the edge by under 0.006%. +const POLAR_LENGTH_STEP: f64 = 0.0175; +/// The most times [`walk_edge`] halves the piece it passes its distance in, to place it there (a +/// 2⁻⁶⁰th of the piece, below an f64's precision); it stops sooner, once the place is within +/// [`POLAR_STRAY_MIN_M`] past it, about 25 halvings for a piece of 1° at 2° from the pole. +const PLACE_HALVINGS: usize = 60; +/// How many cheap reads count as one unit of work: an edge's box found out of an end's reach, which +/// at worst costs about a quarter of a piece measured in [`edge_from`]. Reads add up across ends. +const READS_A_WORK: usize = 4; + +/// The dangle check passed [`MAX_WORK`]: nothing more is compared. +struct WorkLimit; + +/// Adds `n` to `work`, or says it passed [`MAX_WORK`]. +fn charge(work: &mut usize, n: usize) -> Result<(), WorkLimit> { + *work = work.saturating_add(n); + if *work > MAX_WORK { Err(WorkLimit) } else { Ok(()) } +} + +/// Counts a cheap read ([`READS_A_WORK`] to a unit of work). +fn read(reads: &mut usize, work: &mut usize) -> Result<(), WorkLimit> { + *reads = reads.wrapping_add(1); + if reads.is_multiple_of(READS_A_WORK) { charge(work, 1) } else { Ok(()) } +} + +/// The dangle check's grid: which cell a point of the analysis plane falls in. +#[derive(Clone, Copy)] +struct Cells { + /// The lowest corner of the lines' extent, a cell's width and height, and the last cell. + origin: (f64, f64), + size: f64, + last: (i64, i64), +} + +impl Cells { + /// Cells of `size` from `origin` (the lowest corner) to `far` (the highest). + fn new(origin: (f64, f64), size: f64, far: (f64, f64)) -> Self { + let mut cells = Self { origin, size, last: (0, 0) }; + cells.last = cells.key(far.0, far.1); + cells + } + + fn key(&self, x: f64, y: f64) -> (i64, i64) { + (((x - self.origin.0) / self.size).floor() as i64, ((y - self.origin.1) / self.size).floor() as i64) + } +} + +/// Where `c` (longitude/latitude) lies in the azimuthal equidistant projection centered on a pole +/// (the North Pole when `north`): its distance from the pole, at WGS 84's curvature there, in the +/// direction of its longitude. Near the pole that's the ground, within 0.01% along meridians and +/// 0.13% across them out to 85°, right up to the pole and across ±180°. (Not geom's exact `Aeqd`, +/// which solves a geodesic for each point: an edge can be measured in thousands of pieces.) +fn about_pole(c: Coord, north: bool) -> Coord { + let rho = WGS84_POLAR_RADIUS_M * (90.0 - if north { c.y } else { -c.y }).to_radians(); + let (sin, cos) = c.x.to_radians().sin_cos(); + Coord::new(rho * sin, rho * cos) +} + +/// What an end beyond 85.05° can reach in the analysis plane: every point within the tolerance of +/// it in the frame about its pole lies within `wy` of it north–south and, unless the tolerance +/// reaches the pole, within `wx` east–west, also across ±180° (one turn each way: longitudes more +/// than a turn apart aren't compared). +struct Reach { + end: Coord, + wx: f64, + wy: f64, +} + +impl Reach { + /// The reach of `end`, at latitude `lat`, for a tolerance of `tol_m`: within `tol_m` of the end + /// in the frame, a point's distance from the pole is within `tol_m` of the end's, which the + /// plane stretches by at most sec(85.05°) north–south, and its longitude within asin(tol/ρ). + fn new(end: Coord, lat: f64, tol_m: f64) -> Self { + let rho = WGS84_POLAR_RADIUS_M * (90.0 - lat.abs()).to_radians(); + let wy = tol_m / MAX_LAT.to_radians().cos() * 1.001; + let wx = if tol_m < rho { rustgis_geom::MERCATOR_R * (tol_m / rho).asin() * 1.001 } else { f64::INFINITY }; + Self { end, wx, wy } + } + + /// The window's offsets east–west: itself, and a turn either way while it's narrower than the + /// world. + fn shifts(&self) -> impl Iterator { + let round = std::f64::consts::TAU * rustgis_geom::MERCATOR_R; + [0.0, round, -round].into_iter().take(if self.wx.is_finite() { 3 } else { 1 }) + } + + /// Whether the box of the edge `a`–`b` meets the window or one of its copies: a cheap test + /// before [`Reach::spans`]. + fn boxed(&self, a: &Coord, b: &Coord) -> bool { + a.y.min(b.y) <= self.end.y + self.wy + && a.y.max(b.y) >= self.end.y - self.wy + && (self.wx.is_infinite() + || self.shifts().any(|s| a.x.min(b.x) <= self.end.x + s + self.wx && a.x.max(b.x) >= self.end.x + s - self.wx)) + } + + /// The stretches of the edge `a`–`b` (straight in the plane) within reach, one for the window + /// and each of its copies, as ranges of the edge's parameter (0 at `a`, 1 at `b`); none when a + /// coordinate is missing or infinite. + fn spans(&self, a: &Coord, b: &Coord) -> [Option<(f64, f64)>; 3] { + // Where p + (q − p)·t lies within [lo, hi], for t from 0 to 1. + let clip = |lo: f64, hi: f64, p: f64, q: f64| { + if !(p.is_finite() && q.is_finite()) { + return None; + } + let d = q - p; + if d.abs() > 0.0 { + let (t0, t1) = ((lo - p) / d, (hi - p) / d); + let (t0, t1) = (t0.min(t1).max(0.0), t0.max(t1).min(1.0)); + (t0 <= t1).then_some((t0, t1)) + } else { + (lo <= p && p <= hi).then_some((0.0, 1.0)) + } + }; + let mut spans = [None; 3]; + let Some((y0, y1)) = clip(self.end.y - self.wy, self.end.y + self.wy, a.y, b.y) else { return spans }; + if self.wx.is_infinite() { + spans[0] = (a.x.is_finite() && b.x.is_finite()).then_some((y0, y1)); + return spans; + } + for (span, s) in spans.iter_mut().zip(self.shifts()) { + *span = clip(self.end.x + s - self.wx, self.end.x + s + self.wx, a.x, b.x).and_then(|(x0, x1)| { + let (t0, t1) = (x0.max(y0), x1.min(y1)); + (t0 <= t1).then_some((t0, t1)) + }); + } + spans + } +} + +/// How far `e` (placed by [`about_pole`], about the North Pole when `north`) is, squared, from the +/// stretch `t` of the edge `a`–`b` (in the analysis plane). Beyond 85.05° the edge, straight in +/// the plane, is a spiral in the frame: it's measured in pieces short enough that none strays from +/// it by more than a thousandth of `tol_m`, at least a millimeter. Where the stretch reaches into +/// the band it isn't a spiral there, so that bound doesn't hold as such: the pieces are cut to +/// half the allowance, far finer than the edge bends there. The pieces count against `work`; a +/// stretch needing more than [`POLAR_PIECES_MAX`] is measured in that many and sets `coarse`. At +/// t = 0 and 1 the stretch ends at the edge's vertices `la` and `lb` (longitude/latitude), +/// converted as an end is (not at a + (b − a)·1), so an end on a vertex is placed on it. +#[allow(clippy::too_many_arguments)] +fn edge_from( + e: Coord, + (a, b): (Coord, Coord), + (la, lb): (Coord, Coord), + (t0, t1): (f64, f64), + north: bool, + tol_m: f64, + work: &mut usize, + coarse: &mut bool, +) -> Result { + let at = |t: f64| { + if t <= 0.0 { + la + } else if t >= 1.0 { + lb + } else { + plane_to_lonlat(Coord::new(a.x + (b.x - a.x) * t, a.y + (b.y - a.y) * t)) + } + }; + let (p, q) = (at(t0), at(t1)); + // Pieces short enough that none strays from the spiral by more than `stray`. + let colat = |c: Coord| (90.0 - if north { c.y } else { -c.y }).to_radians(); + let rho = WGS84_POLAR_RADIUS_M * colat(p).max(colat(q)); + let drho = WGS84_POLAR_RADIUS_M * (colat(p) - colat(q)).abs(); + let dl = ((b.x - a.x) * (t1 - t0) / rustgis_geom::MERCATOR_R).abs(); + let reaches_band = p.y.abs() <= MAX_LAT || q.y.abs() <= MAX_LAT; + let stray = (tol_m * POLAR_STRAY_SHARE).max(POLAR_STRAY_MIN_M) / if reaches_band { 2.0 } else { 1.0 }; + let n = (crate::geodesic::spiral_stray(rho, drho, dl) / stray).sqrt(); + // Past the cap, the cap; a count that overflows or isn't a number (a coordinate far out of + // range, or missing), one piece. Either is coarse. Along a meridian, one piece, exactly. + let pieces = if !n.is_finite() { + *coarse = true; + 1.0 + } else if n > POLAR_PIECES_MAX { + *coarse = true; + POLAR_PIECES_MAX + } else if n > 1.0 { + n.ceil() + } else { + 1.0 + }; + charge(work, pieces as usize)?; + let (mut last, mut d2, mut k) = (about_pole(p, north), f64::INFINITY, 0.0); + while k < pieces { + k += 1.0; + let t = if k >= pieces { t1 } else { t0 + (t1 - t0) * k / pieces }; + let next = about_pole(at(t), north); + d2 = d2.min(rustgis_data::geometry::seg_dist2(&[last, next], e)); + last = next; + } + Ok(d2) +} + +/// How far a walk along an edge from its start got: [`walk_edge`]. +#[derive(Debug, PartialEq)] +enum Walked { + /// It passed the distance asked for here, the edge's parameter (0 at its start, 1 at its end). + Past(f64), + /// It didn't: the edge is this long (meters). + Short(f64), +} + +/// Where the walk from `a` along the edge `a`–`b` (in the analysis plane), in the frame about a +/// pole (the North Pole when `north`), passes `skip` meters: in pieces of at most +/// [`POLAR_LENGTH_STEP`] of longitude, each counting against `work`, up to the one that passes it, +/// within which the place is found by halving the piece, to a millimeter, each halving counting +/// against `work` too. An edge whose length overflows or isn't a number (a coordinate far out of +/// range, or missing) is one piece, and sets `coarse`, as does one needing more than +/// [`POLAR_PIECES_MAX`], walked in that many, longer (their chords fall shorter, so the walk passes +/// later, never early). +fn walk_edge((a, b): (Coord, Coord), north: bool, skip: f64, work: &mut usize, coarse: &mut bool) -> Result { + let n = ((b.x - a.x) / rustgis_geom::MERCATOR_R).abs() / POLAR_LENGTH_STEP; + let pieces = if !n.is_finite() { + *coarse = true; + 1.0 + } else if n > POLAR_PIECES_MAX { + *coarse = true; + POLAR_PIECES_MAX + } else { + n.ceil().max(1.0) + }; + // At its ends, the edge's own vertices (a + (b − a)·t doesn't come back to them when b − a + // overflows). + let at = |t: f64| { + let c = if t <= 0.0 { + a + } else if t >= 1.0 { + b + } else { + Coord::new(a.x + (b.x - a.x) * t, a.y + (b.y - a.y) * t) + }; + about_pole(plane_to_lonlat(c), north) + }; + let (mut last, mut length, mut k) = (at(0.0), 0.0, 0.0); + while k < pieces { + charge(work, 1)?; + k += 1.0; + let next = at(k / pieces); + let piece = (next.x - last.x).hypot(next.y - last.y); + if length + piece > skip { + // Where in the piece it's `skip - length` from the piece's start, found by halving: its + // parameter isn't in proportion to its length (a piece reaching from beyond 85.05° far + // into the band stretches in the plane about 11.6 times more per meter beyond it). At + // the piece's end when nothing is to be skipped, so a zero tolerance never measures + // the end itself; at or just past the place, never before it. + let (mut lo, mut hi) = ((k - 1.0) / pieces, k / pieces); + if skip > 0.0 { + let rest = skip - length; + // How far past the place `hi` is; halved until it's within a millimeter, each + // halving a unit of work. + let (mut over, mut halvings) = (piece - rest, 0); + while over > POLAR_STRAY_MIN_M && halvings < PLACE_HALVINGS { + let mid = 0.5 * (lo + hi); + if mid <= lo || mid >= hi { + break; + } + halvings += 1; + charge(work, 1)?; + let p = at(mid); + let d = (p.x - last.x).hypot(p.y - last.y); + if d < rest { + lo = mid; + } else { + hi = mid; + over = d - rest; + } + } + } + return Ok(Walked::Past(hi)); + } + length += piece; + last = next; + } + Ok(Walked::Short(length)) +} + +/// The dangle check's ends beyond 85.05° (#1033). There the analysis plane stretches east–west by +/// sec(latitude) but north–south by a constant sec(85.05°), and a pole is a line, so an end is +/// measured in a frame about its pole ([`about_pole`]) instead, right up to the pole, where every +/// longitude meets, and across ±180°. The parts come from the grid's cells within the end's reach, +/// only those holding any (an end near a pole reaches every longitude, a whole row of cells). +struct PolarEnds<'a> { + lines: &'a [(usize, Vec>)], + cells: Cells, + /// Each row's cells that hold parts, by column, in order, with their parts. + occupied: Vec>, + /// Each line's first part's place among all the parts, and the end each part was last opened + /// for (so a part in many cells is measured once an end). + first: Vec, + seen: Vec, + end: usize, + /// Cheap reads so far ([`READS_A_WORK`]). + reads: usize, + /// Whether a stretch was measured in [`POLAR_PIECES_MAX`] pieces, coarser than it needed: the + /// check then says it's incomplete. + coarse: bool, +} + +impl<'a> PolarEnds<'a> { + fn new( + lines: &'a [(usize, Vec>)], + grid: &'a std::collections::HashMap<(i64, i64), Vec<(usize, usize)>>, + cells: Cells, + ) -> Self { + let mut first = Vec::with_capacity(lines.len()); + let mut parts = 0usize; + for (_, ps) in lines { + first.push(parts); + parts = parts.saturating_add(ps.len()); + } + // The grid spans at most 513 rows and columns; a cell beyond them is never looked in. (The + // cap on rows, well past 513, only guards against a grid built otherwise.) + let rows = usize::try_from(cells.last.1.saturating_add(1)).unwrap_or(0).min(1024); + let mut occupied: Vec> = vec![Vec::new(); rows]; + for (&(gx, gy), entries) in grid { + if (0..=cells.last.0).contains(&gx) + && let Some(row) = usize::try_from(gy).ok().and_then(|gy| occupied.get_mut(gy)) + { + row.push((gx, entries.as_slice())); + } + } + occupied.iter_mut().for_each(|row| row.sort_unstable_by_key(|(gx, _)| *gx)); + Self { lines, cells, occupied, first, seen: vec![0; parts], end: 0, reads: 0, coarse: false } + } + + /// Whether `end` (in the analysis plane; `at` in longitude/latitude), an end of part `pi` of + /// line `li` (its first point when `at_start`), comes within `tol_m` of a line: another part, + /// or its own past the stretch within twice the tolerance of the end, as within the band. Only + /// the stretches of edges within the end's reach ([`Reach`]) are measured. Each cell holding + /// parts, entry, part opened, edge whose box meets the reach and piece measured counts against + /// `work`, and edges out of reach count [`READS_A_WORK`] to a unit. The end's own part is read + /// once an end, as within the band: its edges are measured for length (in pieces) until the + /// walk passes twice the tolerance, and the rest, from that point on, as others' are. + fn meets( + &mut self, + (li, pi, at_start): (usize, usize, bool), + (end, at): (Coord, Coord), + tol_m: f64, + work: &mut usize, + ) -> Result { + let north = at.y > 0.0; + let e = about_pole(at, north); + let lines = self.lines; + let reach = Reach::new(end, at.y, tol_m); + // Whether the edge `a`–`b`, from its parameter `from` on, comes within the tolerance. + let mut near = |a: &Coord, b: &Coord, from: f64, work: &mut usize, coarse: &mut bool| -> Result { + if !reach.boxed(a, b) { + return read(&mut self.reads, work).map(|()| false); + } + charge(work, 1)?; + let mut lonlat = None; + for (t0, t1) in reach.spans(a, b).into_iter().flatten() { + let t0 = t0.max(from); + if t0 > t1 { + continue; + } + let ends = *lonlat.get_or_insert_with(|| (plane_to_lonlat(*a), plane_to_lonlat(*b))); + if edge_from(e, (*a, *b), ends, (t0, t1), north, tol_m, work, coarse)? <= tol_m * tol_m { + return Ok(true); + } + } + Ok(false) + }; + // Its own part, from the end on: the rest, past twice the tolerance along it, counts too (a + // line curling back onto itself, like a 6, doesn't dangle there). The rest starts where the + // walk passes that, inside an edge too: beyond 85.05° an edge is a spiral in the frame, and + // one winding round the pole comes back past the end. + let part = lines.get(li).and_then(|(_, ps)| ps.get(pi)).map(Vec::as_slice).unwrap_or(&[]); + let n = part.len(); + let (mut walked, mut past) = (0.0, false); + for i in 0..n.saturating_sub(1) { + let ends = if at_start { (part.get(i), part.get(i + 1)) } else { (part.get(n - 1 - i), part.get(n - 2 - i)) }; + let (Some(a), Some(b)) = ends else { break }; + let skip = 2.0 * tol_m - walked; + let from = if past { + 0.0 + } else if skip <= 0.0 && tol_m > 0.0 { + // The edges before came to exactly twice the tolerance: the rest is all of this one. + past = true; + 0.0 + } else { + match walk_edge((*a, *b), north, skip, work, &mut self.coarse)? { + Walked::Short(length) => { + walked += length; + continue; + } + Walked::Past(from) => { + past = true; + from + } + } + }; + if near(a, b, from, work, &mut self.coarse)? { + return Ok(true); + } + } + let cells = self.cells; + let rows = cells.key(end.x, end.y - reach.wy).1.max(0)..=cells.key(end.x, end.y + reach.wy).1.min(cells.last.1); + // The window's columns and its copies', in order, each column once where they overlap. + let mut columns = [(0, -1); 3]; + for (span, shift) in columns.iter_mut().zip(reach.shifts()) { + let left = cells.key(end.x + shift - reach.wx, end.y).0.max(0); + *span = (left, cells.key(end.x + shift + reach.wx, end.y).0.min(cells.last.0)); + } + columns.sort_unstable(); + let mut done = -1i64; + for span in &mut columns { + span.0 = span.0.max(done.saturating_add(1)); + done = done.max(span.1); + } + self.end = self.end.wrapping_add(1); + for gy in rows { + let Some(row) = usize::try_from(gy).ok().and_then(|gy| self.occupied.get(gy)) else { continue }; + for &(from, to) in &columns { + let start = row.partition_point(|(gx, _)| *gx < from); + for &(_, entries) in row.iter().skip(start).take_while(|(gx, _)| *gx <= to) { + charge(work, 1)?; + for &(oli, opi) in entries { + charge(work, 1)?; + let Some(seen) = self.first.get(oli).and_then(|f| self.seen.get_mut(f + opi)) else { continue }; + if *seen == self.end || (oli, opi) == (li, pi) { + continue; + } + *seen = self.end; + charge(work, 1)?; + let other = lines.get(oli).and_then(|(_, ps)| ps.get(opi)).map(Vec::as_slice).unwrap_or(&[]); + for w in other.windows(2) { + let [a, b] = w else { continue }; + if near(a, b, 0.0, work, &mut self.coarse)? { + return Ok(true); + } + } + } + } + } + } + Ok(false) + } } /// A part without the stretch within `skip` metres (along the line) of one end: from the first @@ -334,15 +801,224 @@ fn rest_of_part(part: &[Coord], from_start: bool, skip: f64) -> Vec { Vec::new() } -/// Points of `ds` that fall in no polygon of `by`. +/// Whether polygons (in the analysis plane, of one feature or several) hold the North Pole (the +/// South Pole when not `north`) inside them together, as the check counts the inside: on the +/// ground they surround it, so in the plane they run along that edge of the world a whole turn, +/// however their longitudes are written (from −180° to 180°, from 0° to 360°, or in parts +/// meeting or overlapping along meridians, written either way). They do if lines 0.13 m, 1 m, +/// 10 m, … 100 km off it in the plane (1.1 cm to 8.7 km on the ground; the nearest just past +/// [`rustgis_geom::POLE_TOLERANCE`], within which the overlays take a polygon to reach the pole) +/// lie inside them all the way around but for gaps of under a meter (9e-6° of longitude), as the +/// overlays' rounding leaves seams centimeters wide in the plane (nothing on the ground beside the +/// pole), from the nearest out to the first wholly outside them, where they end (#1033). A gap +/// reaching the pole (between polygons, in one's outline, or a hole no other polygon fills), +/// however narrow beside it, widens away from it, so it counts once it's wider than a meter on one +/// of the lines; one narrower all the way out to 100 km is taken for a seam. A gap near the pole +/// that doesn't reach it, crossing a line, counts too: the check leans to reporting. A polygon +/// whose crossings with a line can't be worked out (a coordinate missing or infinite there, or so +/// far around, a million turns, that a meter is lost in rounding) doesn't count. When the rest +/// don't surround the pole, or might not have with them (the lines end at one wholly outside the +/// rest that a polygon left out reaches), the rows (`polygons` come with their features' rows) of +/// those left out are returned with the answer, each once, for the check to say; otherwise those +/// left out changed nothing. +fn surrounds_pole<'a>(polygons: impl IntoIterator)>, north: bool) -> (bool, Vec) { + /// The widest gap, in the plane's meters, taken for a seam. + const SEAM_M: f64 = 1.0; + let (rows, polygons): (Vec, Vec<&geo::Polygon>) = polygons.into_iter().unzip(); + let mut skipped = vec![false; polygons.len()]; + let round = std::f64::consts::TAU * rustgis_geom::MERCATOR_R; + let inside = |gaps: &[(f64, f64)]| gaps.iter().all(|(a, b)| b - a <= SEAM_M); + let ended = |gaps: &[(f64, f64)]| gaps.iter().any(|(a, b)| b - a >= round); + let mut lines = [0.13, 1.0, 10.0, 100.0, 1e3, 1e4, 1e5].into_iter(); + let (mut surrounds, mut uncertain) = (false, false); + if let Some(near) = lines.next() { + surrounds = inside(&outside_line(&polygons, north, near, &mut skipped).0); + for off in lines { + if !surrounds { + break; + } + let (gaps, left) = outside_line(&polygons, north, off, &mut skipped); + if ended(&gaps) { + // Where the polygons end, unless one left out reaches the line. + uncertain = left; + break; + } + surrounds = inside(&gaps); + } + } + if surrounds && !uncertain { + return (true, Vec::new()); + } + let mut left_out: Vec = rows.iter().zip(&skipped).filter(|(_, s)| **s).map(|(r, _)| *r).collect(); + left_out.sort_unstable(); + left_out.dedup(); + (surrounds, left_out) +} + +/// The stretches of the line `off` meters off the North Pole (the South Pole when not `north`) in +/// the analysis plane that lie outside `polygons`, within one turn from 0, a stretch across its end +/// taken whole (so it may run past a turn), for [`surrounds_pole`], and whether it left out a +/// polygon. A polygon that may reach the line but that it can't use is marked in `skipped`; one +/// whose latitudes miss the line plays no part, however its longitudes are written. +fn outside_line(polygons: &[&geo::Polygon], north: bool, off: f64, skipped: &mut [bool]) -> (Vec<(f64, f64)>, bool) { + /// The farthest a polygon may lie east or west in the plane (meters, about a million turns), + /// where it still rounds to about a centimeter, well under a seam. + const FAR_M: f64 = 4.0e13; + let round = std::f64::consts::TAU * rustgis_geom::MERCATOR_R; + let pole = lonlat_to_plane(Coord::new(0.0, if north { 90.0 } else { -90.0 })).y; + let y = if north { pole - off } else { pole + off }; + // Each polygon's stretches of the line inside it: inside its outline, less those inside its + // holes (so a hole outside the outline takes nothing), between the places each ring crosses + // the line (worked out exactly, where an overlay would round them to a grid); each brought + // within one turn (split where it wraps). + let mut inside: Vec<(f64, f64)> = Vec::new(); + let (mut outline, mut holes) = (Vec::new(), Vec::new()); + let mut left = false; + for (i, polygon) in polygons.iter().enumerate() { + let Some(r) = polygon.bounding_rect() else { continue }; + let mut skip = || { + left = true; + if let Some(s) = skipped.get_mut(i) { + *s = true; + } + }; + let (y0, y1) = (r.min().y, r.max().y); + if y0.is_finite() && y1.is_finite() && !(y0 <= y && y <= y1) { + continue; + } + let usable = [r.min().x, y0, r.max().x, y1].iter().all(|v| v.is_finite()) && r.min().x >= -FAR_M && r.max().x <= FAR_M; + if !usable { + skip(); + continue; + } + outline.clear(); + holes.clear(); + // A coordinate missing would pair the crossings wrongly. + if !ring_stretches(polygon.exterior(), y, &mut outline) + || !polygon.interiors().iter().all(|hole| ring_stretches(hole, y, &mut holes)) + { + skip(); + continue; + } + holes.sort_by(|a, b| a.0.total_cmp(&b.0)); + for (lo, hi) in less(&outline, &holes) { + if hi - lo >= round { + return (Vec::new(), left); + } + let from = lo.rem_euclid(round); + let to = from + (hi - lo); + inside.push((from, to.min(round))); + if to > round { + inside.push((0.0, to - round)); + } + } + } + // What lies between them, west to east, from the turn's start. + inside.sort_by(|a, b| a.0.total_cmp(&b.0)); + let mut outside = Vec::new(); + let mut reached = 0.0; + for (from, to) in inside { + if from > reached { + outside.push((reached, from)); + } + reached = f64::max(reached, to); + } + if reached < round { + outside.push((reached, round)); + } + // A stretch running across the turn's end and on from its start is one. + if let [first, .., last] = outside.as_mut_slice() + && first.0 <= 0.0 + && last.1 >= round + { + last.1 = round + first.1; + outside.remove(0); + } + (outside, left) +} + +/// Adds to `out` the stretches of the line at `y` inside `ring`, between the places it crosses +/// the line, paired west to east; false when a crossing can't be worked out (a coordinate missing +/// or infinite). +fn ring_stretches(ring: &geo::LineString, y: f64, out: &mut Vec<(f64, f64)>) -> bool { + let mut crossings: Vec = ring + .lines() + .filter(|e| (e.start.y > y) != (e.end.y > y)) + .map(|e| e.start.x + (y - e.start.y) * (e.end.x - e.start.x) / (e.end.y - e.start.y)) + .collect(); + if crossings.iter().any(|x| !x.is_finite()) { + return false; + } + crossings.sort_by(f64::total_cmp); + out.extend(crossings.as_chunks::<2>().0.iter().filter(|[lo, hi]| lo < hi).map(|&[lo, hi]| (lo, hi))); + true +} + +/// The parts of the stretches `of` (apart, west to east, as [`ring_stretches`] gives them) outside +/// every one of `holes` (sorted by where they start), in one pass along both: a ring crossing the +/// line a million times costs no more than its crossings. +fn less(of: &[(f64, f64)], holes: &[(f64, f64)]) -> Vec<(f64, f64)> { + // The holes' stretches joined where they overlap, so they too lie apart, west to east. + let mut apart: Vec<(f64, f64)> = Vec::with_capacity(holes.len()); + for &(lo, hi) in holes { + match apart.last_mut() { + Some(last) if lo <= last.1 => last.1 = last.1.max(hi), + _ => apart.push((lo, hi)), + } + } + let mut out = Vec::new(); + let mut next = 0; + for &(mut lo, hi) in of { + // Holes ending before this stretch end before every later one too. + while apart.get(next).is_some_and(|h| h.1 <= lo) { + next += 1; + } + // Those overlapping it: the last of them may overlap the next stretch, so it's kept. + for &(hole_lo, hole_hi) in apart.get(next..).unwrap_or_default() { + if hole_lo >= hi { + break; + } + if hole_lo > lo { + out.push((lo, hole_lo)); + } + lo = lo.max(hole_hi); + if lo >= hi { + break; + } + } + if lo < hi { + out.push((lo, hi)); + } + } + out +} + +/// Points of `ds` that fall in no polygon of `by`. A point at a pole (to within +/// [`rustgis_geom::POLE_TOLERANCE`]) also falls in the polygons of `by` when together they surround +/// the pole (`surrounds_pole`); ones that only touch it, at a corner, hold it on their boundary, +/// as everywhere else (#1033). Features of `by` the test at a pole had to leave out (a coordinate +/// missing or infinite, or a million turns around) when the rest don't surround it, or might not +/// have with them, are counted in `skipped`, each once at either pole, for the check to say its +/// result there may be wrong. pub fn not_covered(ds: &Dataset, rows: &[usize], by: &Dataset) -> Checked { let all: Vec = (0..by.features.len()).collect(); let polys = polygons(by, &all); + // Whether the polygons together surround each pole, found at the first point there, with the + // features left out. + type Holds = Option<(bool, Vec)>; + let (mut holds_north, mut holds_south): (Holds, Holds) = (None, None); + let with_skipped = |c: Checked, north: &Holds, south: &Holds| { + let mut left_out: Vec = north.iter().chain(south).flat_map(|(_, rows)| rows.iter().copied()).collect(); + left_out.sort_unstable(); + left_out.dedup(); + Checked { skipped: left_out.len(), ..c } + }; + let parts = || polys.iter().flat_map(|(r, mp, _)| mp.0.iter().map(move |p| (*r, p))); let mut out = Vec::new(); crate::control::begin(rows.len()); for (i, r) in rows.iter().enumerate() { if crate::control::cancelled() { - return Checked::partial(out); + return with_skipped(Checked::partial(out), &holds_north, &holds_south); } crate::control::reached(i); let pts: Vec = match ds.features.get(*r).and_then(|f| f.geometry.as_ref()) { @@ -351,17 +1027,23 @@ pub fn not_covered(ds: &Dataset, rows: &[usize], by: &Dataset) -> Checked { _ => continue, }; for c in pts { - let m = lonlat_to_merc(c); + let m = lonlat_to_plane(c); let p = geo::Point::new(m.x, m.y); - if !polys.iter().any(|(_, mp, bb)| rects_touch(bb, &geo::Rect::new(p.0, p.0)) && mp.contains(&p)) { + let surrounded = match crate::convert::reaches_poles(c.y, c.y) { + (true, _) => holds_north.get_or_insert_with(|| surrounds_pole(parts(), true)).0, + (_, true) => holds_south.get_or_insert_with(|| surrounds_pole(parts(), false)).0, + _ => false, + }; + let covered = surrounded || polys.iter().any(|(_, mp, bb)| rects_touch(bb, &geo::Rect::new(p.0, p.0)) && mp.contains(&p)); + if !covered { out.push(TopologyError { rows: vec![*r], at: c }); } } if out.len() >= MAX_ERRORS { - return Checked::partial(out); + return with_skipped(Checked::partial(out), &holds_north, &holds_south); } } - Checked::done(out) + with_skipped(Checked::done(out), &holds_north, &holds_south) } #[cfg(test)] @@ -369,6 +1051,73 @@ mod tests { use super::*; use rustgis_data::GeometryType; + #[test] + fn holes_are_taken_from_an_outline_in_one_pass() { + // Holes overlapping each other, one running across two stretches, one covering a stretch. + let got = less(&[(0.0, 10.0), (20.0, 30.0), (40.0, 50.0)], &[(1.0, 5.0), (2.0, 3.0), (4.0, 8.0), (9.0, 25.0), (40.0, 50.0)]); + assert_eq!(got, [(0.0, 1.0), (8.0, 9.0), (25.0, 30.0)]); + // A million stretches and a million holes between them, as a ring zigzagging across the + // line would make: one pass, not a million times a million. + let at = |i: u32, d: f64| f64::from(i) * 4.0 + d; + let of: Vec<(f64, f64)> = (0..1_000_000).map(|i| (at(i, 0.0), at(i, 2.0))).collect(); + let holes: Vec<(f64, f64)> = (0..1_000_000).map(|i| (at(i, 1.0), at(i, 3.0))).collect(); + let got = less(&of, &holes); + assert!(got.len() == 1_000_000 && got.iter().all(|(lo, hi)| hi - lo == 1.0)); + } + + #[test] + fn a_polygon_left_out_at_a_pole_is_counted() { + // Around the South Pole, a cap from 20°E round to 10°E, and the 10° between them written a + // million turns and more on (1e9° + 90° is 10°E): that one can't be used at the pole, so a + // point there is reported, and the check counts the polygon it left out. A whole cap + // covers the point without looking further, and counts none. + let c = Coord::new; + let ring = |w: f64, e: f64| Geometry::Polygon(vec![vec![vec![c(w, -85.0), c(e, -85.0), c(e, -90.0), c(w, -90.0), c(w, -85.0)]]]); + let cap = || ring(-180.0, 180.0); + let (open, far) = (ring(20.0, 370.0), ring(1e9 + 90.0, 1e9 + 100.0)); + let mut pts = Dataset::new("pts", GeometryType::Point, vec![]); + pts.push(Some(Geometry::Point(c(0.0, -90.0))), vec![]); + let layer = |gs: Vec| { + let mut by = Dataset::new("by", GeometryType::Polygon, vec![]); + for g in gs { + by.push(Some(g), vec![]); + } + by + }; + let checked = not_covered(&pts, &[0], &layer(vec![open, far])); + assert!(checked.errors.len() == 1 && checked.skipped == 1, "{checked:?}"); + let checked = not_covered(&pts, &[0], &layer(vec![cap()])); + assert!(checked.errors.is_empty() && checked.skipped == 0, "{checked:?}"); + // A polygon whose latitudes miss the pole's lines plays no part, however far around it's + // written: with a whole cap, in either order, the point is covered and none is counted. + let at = |w: f64, e: f64, y0: f64, y1: f64| vec![c(w, y0), c(e, y0), c(e, y1), c(w, y1), c(w, y0)]; + let band = Geometry::Polygon(vec![vec![at(1e9, 1e9 + 10.0, 0.0, 10.0)]]); + for gs in [vec![band.clone(), cap()], vec![cap(), band.clone()]] { + let checked = not_covered(&pts, &[0], &layer(gs)); + assert!(checked.errors.is_empty() && checked.skipped == 0, "{checked:?}"); + } + // Left out at both poles, a feature of two parts at one and a feature at the other, beside + // that band: two features, each once. + let open = |y0: f64, y1: f64| Geometry::Polygon(vec![vec![at(20.0, 370.0, y0, y1)]]); + let halves = + |y0: f64, y1: f64| Geometry::Polygon(vec![vec![at(1e9 + 90.0, 1e9 + 95.0, y0, y1)], vec![at(1e9 + 95.0, 1e9 + 100.0, y0, y1)]]); + let mut poles = Dataset::new("pts", GeometryType::Point, vec![]); + poles.push(Some(Geometry::Point(c(0.0, -90.0))), vec![]); + poles.push(Some(Geometry::Point(c(0.0, 90.0))), vec![]); + let by = layer(vec![open(-85.0, -90.0), halves(-85.0, -90.0), band, open(85.0, 90.0), halves(85.0, 90.0)]); + let checked = not_covered(&poles, &[0, 1], &by); + assert!(checked.errors.len() == 2 && checked.skipped == 2, "{checked:?}"); + // A cap 500 m across about the South Pole covers it out to its lines within 87 m; the line + // 860 m off lies wholly outside it, where the check ends. A polygon left out that reaches + // that line might have changed the answer: the point is still covered, but it's counted. + let small = Geometry::Polygon(vec![vec![at(-180.0, 180.0, -90.0, -89.9955)]]); + let checked = not_covered(&pts, &[0], &layer(vec![small.clone()])); + assert!(checked.errors.is_empty() && checked.skipped == 0, "{checked:?}"); + let across = Geometry::Polygon(vec![vec![at(1e9, 1e9 + 10.0, -89.995, -89.99)]]); + let checked = not_covered(&pts, &[0], &layer(vec![small, across])); + assert!(checked.errors.is_empty() && checked.skipped == 1, "{checked:?}"); + } + fn square(x: f64, y: f64, s: f64) -> Geometry { Geometry::Polygon(vec![vec![vec![ Coord::new(x, y), @@ -379,6 +1128,310 @@ mod tests { ]]]) } + #[test] + fn a_polar_end_counts_its_work_and_stops_at_the_limit() { + /// Whether line 0's end at 89°N meets a line, all in one cell, with the work it took + /// (`None`: it passed the limit, starting from `work`). + fn run(lines: &[(usize, Vec>)], work: usize) -> (Option, usize) { + let end = lonlat_to_plane(Coord::new(0.0, 89.0)); + let mut grid = std::collections::HashMap::new(); + grid.insert((0, 0), vec![(0, 0), (1, 0)]); + let cells = Cells::new((end.x - 1e9, end.y - 1e9), 1e10, (end.x, end.y)); + let mut work = work; + let met = PolarEnds::new(lines, &grid, cells).meets((0, 0, false), (end, plane_to_lonlat(end)), 10.0, &mut work); + (met.ok(), work) + } + let p = |x: f64, y: f64| lonlat_to_plane(Coord::new(x, y)); + let own = vec![p(0.0, 88.0), p(0.0, 89.0)]; + // Line 1 meets it there. + assert_eq!(run(&[(0, vec![own.clone()]), (1, vec![vec![p(1.0, 88.0), p(0.0, 89.0)]])], 0).0, Some(true)); + // Line 1 is 1,000 points at 86°N, far off but in the cell: the end measures its own edge's + // length (one piece, 111 km, and 26 halvings placing twice the tolerance in it to a + // millimeter) and the edge from there on (one unit), reads the cell and its two entries, + // opens line 1 (one unit) and finds its 999 edges out of reach (one unit per READS_A_WORK). + let apart = [(0, vec![own]), (1, vec![(0..1_000).map(|i| p(90.0 + f64::from(i) * 1e-4, 86.0)).collect()])]; + assert_eq!(run(&apart, 0), (Some(false), 1 + 26 + 1 + 1 + 2 + 1 + 999 / READS_A_WORK)); + // Starting short of the limit, it stops on the first unit past it: its own edge, a + // halving, the cell, an entry, opening line 1, or reading line 1's edges. + for start in [MAX_WORK, MAX_WORK - 1, MAX_WORK - 2, MAX_WORK - 4, MAX_WORK - 10] { + assert_eq!(run(&apart, start), (None, MAX_WORK + 1), "{start}"); + } + // Reads count across the ends: a second end reads line 1's edges on from the first's 999. + let end = p(0.0, 89.0); + let cells = Cells::new((end.x - 1e9, end.y - 1e9), 1e10, (end.x, end.y)); + let mut grid = std::collections::HashMap::new(); + grid.insert((0, 0), vec![(0, 0), (1, 0)]); + let mut polar = PolarEnds::new(&apart, &grid, cells); + let mut ends = [(None, 0), (None, 0)]; + for (met, work) in &mut ends { + *met = polar.meets((0, 0, false), (end, plane_to_lonlat(end)), 10.0, work).ok(); + } + // Each end: 6 units and its 26 halvings, as above, and the reads. + let first = 6 + 26 + 999 / READS_A_WORK; + assert_eq!(ends, [(Some(false), first), (Some(false), 6 + 26 + (2 * 999 / READS_A_WORK - 999 / READS_A_WORK))]); + // The end's own part's edges out of reach, past twice the tolerance, are read too: nine + // toward and along 86°N, with no other line, after its own first edge's piece (and its 26 + // halvings) and the edge from twice the tolerance on. + let own: Vec = (0..9).map(|i| p(90.0 + f64::from(i) * 1e-4, 86.0)).chain([p(0.0, 88.0), end]).collect(); + let (lone, mut work) = ([(0, vec![own])], 0); + let met = PolarEnds::new(&lone, &std::collections::HashMap::new(), cells).meets( + (0, 0, false), + (end, plane_to_lonlat(end)), + 10.0, + &mut work, + ); + assert_eq!((met.ok(), work), (Some(false), 1 + 26 + 1 + 9 / READS_A_WORK)); + } + + #[test] + fn a_small_layer_is_checked_at_a_large_tolerance() { + // 2,000 lines 100 m long side by side across 89 m, checked at 1 km: the grid's cells are as + // wide as the layer, and each segment goes only into those within it, so the check + // finishes. It wouldn't if each segment filled every cell its widened box covers: about + // 54,000 each, twice the work limit in all. + let mut ds = Dataset::new("l", GeometryType::Polyline, vec![]); + for i in 0..2_000 { + let x = f64::from(i) * 4e-7; + ds.push(Some(Geometry::Line(vec![vec![Coord::new(x, 0.0), Coord::new(x, 0.0009)]])), vec![]); + } + let checked = dangles(&ds, &(0..2_000).collect::>(), 1_000.0); + assert!(checked.complete && checked.errors.is_empty(), "complete {}, {} errors", checked.complete, checked.errors.len()); + } + + #[test] + fn an_edge_with_a_coordinate_out_of_range_is_one_piece() { + // Each is one piece, and coarse; so is a stretch of one in edge_from, a zero-length + // stretch of an edge whose length overflows included (∞ × 0 is no number). + let b = lonlat_to_plane(Coord::new(0.0, 89.0)); + for a in [Coord::new(f64::INFINITY, b.y), Coord::new(f64::NAN, b.y)] { + let (mut work, mut coarse) = (0, false); + let _ = walk_edge((a, b), true, 0.0, &mut work, &mut coarse); + assert!(work == 1 && coarse, "{a:?}"); + let (mut work, mut coarse, ends) = (0, false, (plane_to_lonlat(a), plane_to_lonlat(b))); + let e = about_pole(Coord::new(0.0, 89.0), true); + assert!(edge_from(e, (a, b), ends, (0.0, 1.0), true, 1.0, &mut work, &mut coarse).is_ok() && coarse && work == 1, "{a:?}"); + } + let (a, b) = (Coord::new(-f64::MAX, b.y), Coord::new(f64::MAX, b.y)); + let (mut work, mut coarse, ends) = (0, false, (plane_to_lonlat(a), plane_to_lonlat(b))); + let e = about_pole(Coord::new(0.0, 89.0), true); + assert!(edge_from(e, (a, b), ends, (0.5, 0.5), true, 1.0, &mut work, &mut coarse).is_ok() && coarse && work == 1); + } + + #[test] + fn a_long_edge_is_walked_only_as_far_as_it_needs() { + // An edge winding round the North Pole about 194 times (70,000° of longitude), more than + // POLAR_PIECES_MAX pieces of POLAR_LENGTH_STEP: it's walked in that many, and says so. The + // walk from its start passes 20 m within its first piece, and stops there. A stretch of + // it needing more than POLAR_PIECES_MAX pieces is measured in that many, and says so too. + let (a, b) = (lonlat_to_plane(Coord::new(0.0, 88.0)), lonlat_to_plane(Coord::new(70_000.0, 88.0))); + assert!(70_000f64.to_radians() / POLAR_LENGTH_STEP > POLAR_PIECES_MAX); + let (mut work, mut coarse) = (0, false); + let Ok(Walked::Past(t)) = walk_edge((a, b), true, 20.0, &mut work, &mut coarse) else { panic!("not past") }; + let passed = about_pole(plane_to_lonlat(Coord::new(a.x + (b.x - a.x) * t, a.y)), true); + let start = about_pole(plane_to_lonlat(a), true); + let past = (passed.x - start.x).hypot(passed.y - start.y) - 20.0; + assert!(t < 1.0 / POLAR_PIECES_MAX && (0.0..=POLAR_STRAY_MIN_M).contains(&past), "{t} {past}"); + // One piece, and the halvings placing 20 m within it. + let walked = work; + assert!(coarse && walked > 1 && walked <= 1 + PLACE_HALVINGS, "{walked}"); + let ends = (plane_to_lonlat(a), plane_to_lonlat(b)); + coarse = false; + assert!(edge_from(start, (a, b), ends, (0.0, 1.0), true, 1.0, &mut work, &mut coarse).is_ok()); + assert!(coarse && work == walked + POLAR_PIECES_MAX as usize, "{work}"); + // Under the cap (600° of longitude), the walk is as fine as it needs, and not coarse. + let (mut work, mut coarse) = (0, false); + let short = lonlat_to_plane(Coord::new(600.0, 88.0)); + assert!(matches!(walk_edge((a, short), true, 20.0, &mut work, &mut coarse), Ok(Walked::Past(_))) && !coarse); + } + + #[test] + fn an_edge_whose_piece_count_overflows_is_measured_as_one_coarse_piece() { + // An edge at 88°N running 1e300 m across the plane (coordinates far out of range): its + // piece count overflows. It's measured as one piece, for one unit of work, and is coarse, + // not cut into POLAR_PIECES_MAX pieces that measure nothing better. + let a = lonlat_to_plane(Coord::new(0.0, 88.0)); + let b = Coord::new(1e300, a.y); + let (ends, mut work, mut coarse) = ((plane_to_lonlat(a), plane_to_lonlat(b)), 0, false); + let e = about_pole(Coord::new(0.0, 88.0), true); + assert!(edge_from(e, (a, b), ends, (0.0, 1.0), true, 1.0, &mut work, &mut coarse).is_ok()); + assert!(coarse && work == 1, "{work}"); + } + + #[test] + fn a_stretch_measured_coarsely_leaves_the_check_incomplete_not_stopped() { + // A line 11 m from the North Pole winding round it about 19,000 times (7,000,000° of + // longitude): each end, within 20 m of the pole, reaches the whole of it, which needs about + // a million pieces. It's measured in POLAR_PIECES_MAX, and the check says it's incomplete; + // the other line's ends are still checked, and dangle. + let c = Coord::new; + let mut ds = Dataset::new("l", GeometryType::Polyline, vec![]); + ds.push(Some(Geometry::Line(vec![vec![c(0.0, 89.9999), c(7_000_000.0, 89.9999)]])), vec![]); + ds.push(Some(Geometry::Line(vec![vec![c(0.0, 86.0), c(0.001, 86.0)]])), vec![]); + let checked = dangles(&ds, &[0, 1], 20.0); + let other: Vec = checked.errors.iter().filter(|e| e.rows == [1]).map(|e| e.at).collect(); + assert!(!checked.complete && other.len() == 2, "complete {}, {other:?}", checked.complete); + // Measured as closely as it needs, a line winding a few times is complete. + ds.features.clear(); + ds.push(Some(Geometry::Line(vec![vec![c(0.0, 89.9999), c(2_000.0, 89.9999)]])), vec![]); + assert!(dangles(&ds, &[0], 20.0).complete); + } + + #[test] + fn a_polar_end_meets_its_own_line_past_twice_the_tolerance_within_an_edge() { + // A hairpin at 86°N: 14.8 m east from E, then back west through E for 3°. The line passes + // E 29.6 m along it, past twice the 10 m tolerance, inside the second edge's first piece + // (7.8 km): E doesn't dangle, only the far end does. + let mut ds = Dataset::new("l", GeometryType::Polyline, vec![]); + ds.push(Some(Geometry::Line(vec![vec![Coord::new(0.0, 86.0), Coord::new(0.0019, 86.0), Coord::new(-3.0, 86.0)]])), vec![]); + let checked = dangles(&ds, &[0], 10.0); + let at: Vec = checked.errors.iter().map(|e| e.at).collect(); + assert!(checked.complete && at.len() == 1 && at.iter().all(|c| (c.x + 3.0).abs() < 1e-9), "{at:?}"); + } + + #[test] + fn a_polar_end_on_a_line_reaching_far_into_the_band_dangles() { + // Lines down a meridian, or nearly, from beyond 85.05° far into the band: each is one + // piece, along which the plane stretches about 11.6 times more per meter beyond 85.05° + // than within it. Twice the tolerance is found on the ground, not in proportion to the + // plane, so the polar end doesn't meet its own line a few meters on, and both ends dangle. + let ends = [((0.0, 86.0), (0.0, 60.0)), ((0.0, -86.0), (0.0, -60.0)), ((0.0, 89.0), (0.0, 0.0)), ((0.0, 86.0), (0.5, 60.0))]; + for ((x0, y0), (x1, y1)) in ends { + let mut ds = Dataset::new("l", GeometryType::Polyline, vec![]); + ds.push(Some(Geometry::Line(vec![vec![Coord::new(x0, y0), Coord::new(x1, y1)]])), vec![]); + let checked = dangles(&ds, &[0], 10.0); + assert!(checked.complete && checked.errors.len() == 2, "({x0}, {y0}) to ({x1}, {y1}): {:?}", checked.errors); + } + // The walk passes 20 m from the end on the ground, in the frame about the pole. + let (a, b) = (lonlat_to_plane(Coord::new(0.0, 86.0)), lonlat_to_plane(Coord::new(0.0, 60.0))); + let (mut work, mut coarse) = (0, false); + let Ok(Walked::Past(t)) = walk_edge((a, b), true, 20.0, &mut work, &mut coarse) else { panic!("not past") }; + let at = about_pole(plane_to_lonlat(Coord::new(a.x, a.y + (b.y - a.y) * t)), true); + let start = about_pole(plane_to_lonlat(a), true); + let past = (at.x - start.x).hypot(at.y - start.y) - 20.0; + assert!((0.0..=POLAR_STRAY_MIN_M).contains(&past), "{t} {past}"); + } + + #[test] + fn an_edge_whose_length_overflows_leaves_the_check_incomplete() { + // An edge along 88°N from −9.2e302° to 9.2e302° of longitude: its length in the plane + // overflows. It's walked as one piece between its own vertices, to a length that's a + // number, and the check says it's incomplete rather than measuring it as nothing. + let (w, e) = (Coord::new(-9.2e302, 88.0), Coord::new(9.2e302, 88.0)); + let (a, b) = (lonlat_to_plane(w), lonlat_to_plane(e)); + assert!(a.is_finite() && b.is_finite() && (b.x - a.x).is_infinite()); + let (mut work, mut coarse) = (0, false); + let walked = walk_edge((a, b), true, f64::INFINITY, &mut work, &mut coarse); + assert!(coarse && work == 1 && matches!(walked, Ok(Walked::Short(l)) if l.is_finite()), "{:?}", walked.ok()); + let mut ds = Dataset::new("l", GeometryType::Polyline, vec![]); + ds.push(Some(Geometry::Line(vec![vec![w, e]])), vec![]); + assert!(!dangles(&ds, &[0], 10.0).complete); + } + + #[test] + fn a_polar_end_meets_its_own_line_from_exactly_twice_the_tolerance() { + // The hairpin above at a tolerance of half its first edge's length: the walk comes to + // twice the tolerance exactly at that edge's end, so all of the second edge counts, and it + // passes back through E within its first piece. E doesn't dangle, only the far end does. + let (e, f, g) = (Coord::new(0.0, 86.0), Coord::new(0.0019, 86.0), Coord::new(-3.0, 86.0)); + let (mut work, mut coarse) = (0, false); + let first = (lonlat_to_plane(e), lonlat_to_plane(f)); + let Ok(Walked::Short(length)) = walk_edge(first, true, f64::INFINITY, &mut work, &mut coarse) else { panic!("not short") }; + let mut ds = Dataset::new("l", GeometryType::Polyline, vec![]); + ds.push(Some(Geometry::Line(vec![vec![e, f, g]])), vec![]); + let checked = dangles(&ds, &[0], length / 2.0); + let at: Vec = checked.errors.iter().map(|e| e.at).collect(); + assert!(checked.complete && at.len() == 1 && at.iter().all(|c| (c.x + 3.0).abs() < 1e-9), "{at:?}"); + } + + #[test] + fn an_edge_winding_round_a_pole_meets_its_own_end() { + // One edge along 88°N through 400° of longitude passes back over both its ends (0° and, + // a turn on, 40°): neither dangles, though each is far, along the edge, from where it + // passes. Through 300° it doesn't come back, and both do. + let line = |to: f64| { + let mut ds = Dataset::new("l", GeometryType::Polyline, vec![]); + ds.push(Some(Geometry::Line(vec![vec![Coord::new(0.0, 88.0), Coord::new(to, 88.0)]])), vec![]); + dangles(&ds, &[0], 10.0) + }; + let wound = line(400.0); + assert!(wound.complete && wound.errors.is_empty(), "{:?}", wound.errors); + let open = line(300.0); + assert!(open.complete && open.errors.len() == 2, "{:?}", open.errors); + } + + #[test] + fn a_polar_end_reads_only_the_cells_holding_parts() { + // An end 1.1 m from the North Pole reaches every longitude: a whole row of 512 cells, of + // which only the last holds a part (line 1, far off). The end measures its own edge's + // length (one piece, 11 km, and 24 halvings placing twice the tolerance in it to a + // millimeter) and the edge from there on, and reads that cell, its entry and line 1 (a + // unit each); the empty cells cost nothing, and cells beyond the grid's rows and columns + // aren't looked in. + let p = |x: f64, y: f64| lonlat_to_plane(Coord::new(x, y)); + let end = p(0.0, 90.0 - 1e-5); + let lines = [(0, vec![vec![p(0.0, 89.9), end]]), (1, vec![vec![p(90.0, 86.0), p(90.1, 86.0)]])]; + let mut grid = std::collections::HashMap::new(); + for cell in [(511, 0), (-1, 0), (512, 0), (511, 1)] { + grid.insert(cell, vec![(1, 0)]); + } + let cells = Cells::new((end.x - 1.0, end.y - 0.5e9), 1e9, (end.x - 1.0 + 511.5e9, end.y)); + assert_eq!(cells.last, (511, 0)); + let mut work = 0; + let met = PolarEnds::new(&lines, &grid, cells).meets((0, 0, false), (end, plane_to_lonlat(end)), 10.0, &mut work); + assert_eq!((met.ok(), work), (Some(false), 5 + 24)); + } + + #[test] + fn a_polar_edge_is_measured_along_its_spiral() { + // Straight in the plane from 86°N to 89°N, turning 120° about the pole: in the frame about + // the pole it's a spiral, far from the chord between its ends. + let (la, lb) = (Coord::new(0.0, 86.0), Coord::new(120.0, 89.0)); + let (a, b) = (lonlat_to_plane(la), lonlat_to_plane(lb)); + let mid = about_pole(plane_to_lonlat(Coord::new((a.x + b.x) / 2.0, (a.y + b.y) / 2.0)), true); + let chord = rustgis_data::geometry::seg_dist2(&[about_pole(la, true), about_pole(lb, true)], mid).sqrt(); + assert!(chord > 10_000.0, "{chord}"); + // A 10 m tolerance allows pieces to stray 1 cm. + let (mut work, mut coarse, ends) = (0, false, (plane_to_lonlat(a), plane_to_lonlat(b))); + let d = edge_from(mid, (a, b), ends, (0.0, 1.0), true, 10.0, &mut work, &mut coarse).map(f64::sqrt).unwrap_or(f64::NAN); + assert!(d < 0.01 && work > 1 && work <= POLAR_PIECES_MAX as usize && !coarse, "{d} m, {work} pieces"); + // A vertex is on its edge (within a micrometer), and a stretch is measured on its own. + let at_b = edge_from(about_pole(ends.1, true), (a, b), ends, (0.0, 1.0), true, 10.0, &mut work, &mut coarse).unwrap_or(f64::NAN); + assert!(at_b < 1e-12, "{at_b}"); + let first_half = edge_from(about_pole(ends.1, true), (a, b), ends, (0.0, 0.5), true, 10.0, &mut work, &mut coarse).unwrap_or(0.0); + assert!(first_half.sqrt() > 10_000.0, "{first_half}"); + } + + #[test] + fn reach_finds_the_stretches_of_an_edge_within_it() { + let round = std::f64::consts::TAU * rustgis_geom::MERCATOR_R; + let c = Coord::new; + let reach = Reach { end: c(0.0, 0.0), wx: 10.0, wy: 5.0 }; + // Each edge, and its stretches within the window, a turn east of it and a turn west of it. + let cases = [ + (c(-20.0, 0.0), c(20.0, 0.0), [Some((0.25, 0.75)), None, None]), + (c(0.0, -10.0), c(20.0, 10.0), [Some((0.25, 0.5)), None, None]), + (c(-20.0, 10.0), c(20.0, 10.0), [None, None, None]), + (c(round - 20.0, 0.0), c(round, 0.0), [None, Some((0.5, 1.0)), None]), + (c(-round, 0.0), c(round, 0.0), [Some((0.5, 0.5)), Some((1.0, 1.0)), Some((0.0, 0.0))]), + (c(f64::NAN, 0.0), c(20.0, 0.0), [None, None, None]), + (c(0.0, f64::NAN), c(20.0, 0.0), [None, None, None]), + ]; + for (a, b, want) in cases { + let got = reach.spans(&a, &b); + let close = got.iter().zip(&want).all(|(g, w)| match (g, w) { + (Some(g), Some(w)) => (g.0 - w.0).abs() < 1e-6 && (g.1 - w.1).abs() < 1e-6, + (g, w) => g.is_none() && w.is_none(), + }); + assert!(close, "{a:?}–{b:?}: {got:?}"); + // The box test never turns away an edge with a stretch in reach. + assert!(got.iter().all(Option::is_none) || reach.boxed(&a, &b), "{a:?}–{b:?}"); + } + // A tolerance reaching the pole reaches every longitude. + let around = Reach { end: c(0.0, 0.0), wx: f64::INFINITY, wy: 5.0 }; + assert_eq!(around.spans(&c(-1e9, -10.0), &c(1e9, 10.0)), [Some((0.25, 0.75)), None, None]); + assert!(Reach::new(c(0.0, 0.0), 89.999_99, 10.0).wx.is_infinite() && Reach::new(c(0.0, 0.0), 89.0, 10.0).wx.is_finite()); + } + #[test] fn checks_find_problems() { let mut ds = Dataset::new("p", GeometryType::Polygon, vec![]); diff --git a/crates/analysis/tests/polar.rs b/crates/analysis/tests/polar.rs new file mode 100644 index 000000000..933ab94f8 --- /dev/null +++ b/crates/analysis/tests/polar.rs @@ -0,0 +1,765 @@ +//! Polar features keep their shape and their relations through the analysis tools (#1033): +//! nothing is clamped onto 85.05°, where Web Mercator ends. +use rustgis_analysis::{Inputs, Output, ToolResult}; +use rustgis_data::{Coord, Dataset, Field, FieldType, Geometry, GeometryType, Value}; +use std::collections::HashMap; + +struct Map(HashMap); +impl Inputs for Map { + fn dataset(&self, layer: &str) -> Option<&Dataset> { + self.0.get(layer) + } + fn selection(&self, _: &str) -> &[usize] { + &[] + } +} + +fn square(x0: f64, y0: f64, x1: f64, y1: f64) -> Geometry { + let c = Coord::new; + Geometry::Polygon(vec![vec![vec![c(x0, y0), c(x1, y0), c(x1, y1), c(x0, y1), c(x0, y0)]]]) +} + +fn layer(name: &str, geometries: Vec) -> (String, Dataset) { + let mut ds = Dataset::new(name, GeometryType::Polygon, vec![Field::new("k", FieldType::Text)]); + for g in geometries { + ds.push(Some(g), vec![Value::Text("a".into())]); + } + (name.to_string(), ds) +} + +/// Two overlapping polygons between 87°S and 89°S, a mask over them, and polygons from 86°S to 87°S +/// and 88°S to 89°S, about 111 km apart. +fn polar() -> Map { + Map(HashMap::from([ + layer("polar", vec![square(0.0, -89.0, 1.0, -87.0), square(0.5, -89.0, 1.5, -87.0)]), + layer("mask", vec![square(-10.0, -89.9, 10.0, -86.0)]), + layer("upper", vec![square(0.0, -87.0, 1.0, -86.0)]), + layer("lower", vec![square(0.0, -89.0, 1.0, -88.0)]), + ])) +} + +fn dataset(r: ToolResult) -> Dataset { + match r.output { + Output::Layer(d) | Output::Table(d) => d, + _ => panic!("expected a layer"), + } +} + +/// Every latitude of every feature's geometry, with the feature count. +fn latitudes(ds: &Dataset) -> (usize, Vec) { + let mut lats = Vec::new(); + for f in ds.features.iter() { + if let Some(g) = &f.geometry { + g.for_each_coord(|c| lats.push(c.y)); + } + } + (ds.features.len(), lats) +} + +fn run(m: &Map, tool: &str, params: serde_json::Value) -> Dataset { + dataset(rustgis_analysis::run(tool, m, ¶ms).unwrap_or_else(|e| panic!("{tool}: {e}"))) +} + +#[test] +fn polar_polygons_keep_their_shape_through_the_tools() { + let m = polar(); + let between = |lats: &[f64]| !lats.is_empty() && lats.iter().all(|y| (-89.0 - 1e-9..=-87.0 + 1e-9).contains(y)); + let reach = |lats: &[f64], y: f64| lats.iter().any(|v| (v - y).abs() < 1e-9); + + // Clip by a mask that covers them: both come through whole. + let (n, lats) = latitudes(&run(&m, "clip", serde_json::json!({"input": "polar", "clip": "mask"}))); + assert!(n == 2 && between(&lats) && reach(&lats, -89.0) && reach(&lats, -87.0), "clip: {n} {lats:?}"); + // Dissolve unions them into one polygon spanning both (to the overlay's grid, under a + // millimeter: 1e-7° here). + let d = run(&m, "dissolve", serde_json::json!({"input": "polar", "field": "k"})); + let (n, lats) = latitudes(&d); + let e = d.features.first().and_then(|f| f.geometry.as_ref()).map(|g| g.extent()).unwrap(); + assert!(n == 1 && between(&lats) && e.min_x.abs() < 1e-7 && (e.max_x - 1.5).abs() < 1e-7, "dissolve: {lats:?} {e:?}"); + // Simplify keeps a square's corners, hull and points stay inside. + let (_, lats) = latitudes(&run(&m, "simplify", serde_json::json!({"input": "polar", "tolerance": 0.001}))); + assert!(between(&lats) && reach(&lats, -89.0) && reach(&lats, -87.0), "simplify: {lats:?}"); + let (_, lats) = latitudes(&run(&m, "convexHull", serde_json::json!({"input": "polar"}))); + assert!(between(&lats) && reach(&lats, -89.0) && reach(&lats, -87.0), "convexHull: {lats:?}"); + let (n, lats) = latitudes(&run(&m, "featureToPoint", serde_json::json!({"input": "polar"}))); + assert!(n == 2 && lats.iter().all(|y| (-89.0..-87.0).contains(y)), "featureToPoint: {lats:?}"); +} + +#[test] +fn a_cap_around_a_pole_keeps_its_pole_through_the_overlays() { + let c = Coord::new; + // As files cut a cap around the South Pole: along 85°S, down to the pole at 180°, and back. + let cap = Geometry::Polygon(vec![vec![vec![c(-180.0, -85.0), c(180.0, -85.0), c(180.0, -90.0), c(-180.0, -90.0), c(-180.0, -85.0)]]]); + let m = Map(HashMap::from([layer("cap", vec![cap]), layer("mask", vec![square(-180.0, -90.0, 180.0, -80.0)])])); + // Within the overlays' grid of the pole and of 85°S, and nowhere in between them. + let pole = |lats: &[f64]| lats.iter().any(|y| (y + 90.0).abs() <= rustgis_geom::POLE_TOLERANCE); + let edge = |lats: &[f64]| lats.iter().any(|y| (y + 85.0).abs() < 1e-7); + for (tool, params) in [ + ("clip", serde_json::json!({"input": "cap", "clip": "mask"})), + ("dissolve", serde_json::json!({"input": "cap", "field": "k"})), + ("simplify", serde_json::json!({"input": "cap", "tolerance": 0.001})), + ] { + let (n, lats) = latitudes(&run(&m, tool, params)); + assert!(n == 1 && pole(&lats) && edge(&lats) && lats.iter().all(|y| (-90.0 - 1e-9..=-85.0 + 1e-7).contains(y)), "{tool}: {lats:?}"); + } +} + +#[test] +fn polar_polygons_apart_are_apart() { + let m = polar(); + // 86–87°S and 88–89°S don't meet: they did when both were clamped onto 85.05°S. + let r = rustgis_analysis::run("selectByLocation", &m, &serde_json::json!({"input": "upper", "select": "lower"})).unwrap(); + let Output::Selection { rows, .. } = r.output else { panic!("expected a selection") }; + assert!(rows.is_empty(), "{rows:?}"); + // The closest one is 1° of latitude away near 88°S: 111.69 km on the WGS 84 meridian, not 0. + let joined = + run(&m, "spatialJoin", serde_json::json!({"target": "upper", "join": "lower", "match": "Closest", "distanceField": "gap"})); + let gap = joined.field_index("gap").map(|i| joined.value(0, i).clone()); + let km = match gap { + Some(Value::Float(v)) => v, + other => panic!("gap: {other:?}"), + }; + assert!((km - 111.69).abs() < 0.05, "{km}"); +} + +#[test] +fn line_ends_near_a_pole_meet_within_the_tolerance() { + // Two short lines at 89.6°N, 500 m apart in the plane, 3.5 m on the ground: no dangles at 10 m. + let c = Coord::new; + let dx = (500.0 / rustgis_geom::MERCATOR_R).to_degrees(); + let mut ds = Dataset::new("lines", GeometryType::Polyline, vec![]); + for x in [0.0, dx] { + ds.push(Some(Geometry::Line(vec![vec![c(x, 89.6), c(x, 89.6001)]])), vec![]); + } + let checked = rustgis_analysis::topology::dangles(&ds, &[0, 1], 10.0); + assert!(checked.complete && checked.errors.is_empty(), "{:?}", checked.errors.iter().map(|e| e.at).collect::>()); + + // The same lines 5 km apart in the plane, about 35 m on the ground: all four ends dangle. + // (Clamped onto 85.05°, each line was a single point, and none was found.) + let dx = (5_000.0 / rustgis_geom::MERCATOR_R).to_degrees(); + let mut ds = Dataset::new("lines", GeometryType::Polyline, vec![]); + for x in [0.0, dx] { + ds.push(Some(Geometry::Line(vec![vec![c(x, 89.6), c(x, 89.6001)]])), vec![]); + } + let checked = rustgis_analysis::topology::dangles(&ds, &[0, 1], 10.0); + assert!(checked.complete && checked.errors.len() == 4, "{:?}", checked.errors.iter().map(|e| e.at).collect::>()); + + // Lines end to end along 89.9°N, where a degree of longitude is 194 m: 5 m apart on the ground + // they meet, and only their outer ends dangle; 30 m apart, all four do. + let per_degree = rustgis_geom::MERCATOR_R * 89.9f64.to_radians().cos() * std::f64::consts::PI / 180.0; + for (gap_m, dangling) in [(5.0, 2), (30.0, 4)] { + let gap = gap_m / per_degree; + let mut ds = Dataset::new("lines", GeometryType::Polyline, vec![]); + ds.push(Some(Geometry::Line(vec![vec![c(0.0, 89.9), c(1.0, 89.9)]])), vec![]); + ds.push(Some(Geometry::Line(vec![vec![c(1.0 + gap, 89.9), c(2.0, 89.9)]])), vec![]); + let checked = rustgis_analysis::topology::dangles(&ds, &[0, 1], 10.0); + let at: Vec = checked.errors.iter().map(|e| e.at).collect(); + assert!(checked.complete && at.len() == dangling, "{gap_m} m: {at:?}"); + } + + // Meridians at 0°E and 90°E from 80°N meet at the North Pole, where every longitude is one + // place: only their ends at 80°N dangle. + let mut ds = Dataset::new("meridians", GeometryType::Polyline, vec![]); + for x in [0.0, 90.0] { + ds.push(Some(Geometry::Line(vec![vec![c(x, 80.0), c(x, 90.0)]])), vec![]); + } + let checked = rustgis_analysis::topology::dangles(&ds, &[0, 1], 10.0); + let at: Vec = checked.errors.iter().map(|e| e.at.y).collect(); + assert!(checked.complete && at.len() == 2 && at.iter().all(|y| (y - 80.0).abs() < 1e-9), "{at:?}"); +} + +#[test] +fn line_ends_near_a_pole_are_looked_up_not_scanned() { + // 7,000 lines end to end around 88°S, 200 m each: every end meets its neighbor's, the chain's + // first and last too, across ±180°. Each end is measured only against the lines near it, so the + // check finishes (against every line, it would pass the work limit). + let c = Coord::new; + let n = 7_000; + let lon = |i: i32| -180.0 + 360.0 * f64::from(i) / f64::from(n); + let mut ds = Dataset::new("chain", GeometryType::Polyline, vec![]); + for i in 0..n { + ds.push(Some(Geometry::Line(vec![vec![c(lon(i), -88.0), c(lon(i + 1), -88.0)]])), vec![]); + } + let rows: Vec = (0..ds.features.len()).collect(); + let checked = rustgis_analysis::topology::dangles(&ds, &rows, 10.0); + let at: Vec = checked.errors.iter().map(|e| e.at).collect(); + assert!(checked.complete && at.is_empty(), "{} found: {:?}", at.len(), at.iter().take(4).collect::>()); + + // 8,000 short lines on top of one another at 89°N: each end meets the first other line it + // looks at, so the check finishes (measuring every line against every end, it wouldn't). + let mut ds = Dataset::new("pile", GeometryType::Polyline, vec![]); + for _ in 0..8_000 { + ds.push(Some(Geometry::Line(vec![vec![c(0.0, 89.0), c(0.001, 89.0)]])), vec![]); + } + let rows: Vec = (0..ds.features.len()).collect(); + let checked = rustgis_analysis::topology::dangles(&ds, &rows, 10.0); + assert!(checked.complete && checked.errors.is_empty(), "{} found", checked.errors.len()); +} + +/// How many line ends dangle at `tol_m` among `lines`, each a list of (longitude, latitude). +fn dangling(lines: &[&[(f64, f64)]], tol_m: f64) -> usize { + let mut ds = Dataset::new("lines", GeometryType::Polyline, vec![]); + for line in lines { + ds.push(Some(Geometry::Line(vec![line.iter().map(|&(x, y)| Coord::new(x, y)).collect()])), vec![]); + } + let rows: Vec = (0..ds.features.len()).collect(); + let checked = rustgis_analysis::topology::dangles(&ds, &rows, tol_m); + assert!(checked.complete); + checked.errors.len() +} + +#[test] +fn line_ends_close_to_a_pole_are_measured_across_it() { + // Meridian lines from 89.9°, ending `rho` meters from the pole (the North Pole, or the South + // Pole when `s` is -1): their ends there meet (only the two at 89.9° dangle) when within the + // tolerance of the other line, through the pole. + let from_pole = |rho_m: f64| 90.0 - degrees(rho_m); + for (lon_a, rho_a, lon_b, rho_b, s, tol, want) in [ + (0.0, 11.0, 40.0, 2.0, 1.0, 10.0, 2), // a's end 7.07 m from b's line (9.55 m from its end) + (0.0, 11.0, 60.0, 5.5, 1.0, 10.0, 2), // 9.53 m from b's line + (0.0, 2.0, 90.0, 11.0, 1.0, 10.0, 4), // 11.18 m from end to end + (0.0, 4.5, 180.0, 4.5, 1.0, 10.0, 2), // 9.0 m, through the pole + (0.0, 4.5, 180.0, 4.5, -1.0, 10.0, 2), // the same through the South Pole + (0.0, 1.0, 180.0, 10.0, 1.0, 10.0, 4), // 11.0 m + (0.0, 11.1, 180.0, 0.0, 1.0, 15.0, 2), // one at the pole, 11.1 m away + ] { + let a = [(lon_a, s * 89.9), (lon_a, s * from_pole(rho_a))]; + let b = [(lon_b, s * 89.9), (lon_b, s * from_pole(rho_b))]; + let got = dangling(&[&a, &b], tol); + assert_eq!(got, want, "{lon_a}° {rho_a} m and {lon_b}° {rho_b} m at {tol} m ({s})"); + } + // A lone line reaching the pole still dangles there. + assert_eq!(dangling(&[&[(0.0, 89.0), (0.0, 90.0)]], 10.0), 2); +} + +/// Degrees of latitude in `m` meters near a pole (WGS 84's radius of curvature there, a²/b). +fn degrees(m: f64) -> f64 { + (m / rustgis_geom::WGS84_POLAR_RADIUS_M).to_degrees() +} + +#[test] +fn line_ends_near_a_pole_are_measured_to_the_edge_as_it_runs() { + // Beyond 85.05° an edge straight in the plane curves in the frame about the pole that ends are + // measured in: along a parallel it bows away from the pole (by up to 0.2 m over 0.09° at 85°), + // and running toward the pole it spirals, so it's measured in pieces short enough to follow it + // to a millimeter. An end 1.1 m poleward of an edge's middle dangles at 1 m (the default + // tolerance), and an end on an edge along a parallel meets it at 0.1 m. (At 84.9°, within Web + // Mercator's square, ends are checked as before.) + for lat in [84.9, 85.06, 87.0, 89.0] { + let edge = [(0.0, lat), (0.09, lat + degrees(11.0))]; + let near = lat + degrees(5.5 + 1.1); + assert_eq!(dangling(&[&edge, &[(0.045, near), (0.045, near + degrees(100.0))]], 1.0), 4, "{lat}"); + let parallel = [(0.0, lat), (0.09, lat)]; + assert_eq!(dangling(&[¶llel, &[(0.045, lat), (0.045, lat - degrees(100.0))]], 0.1), 3, "{lat}"); + } + // An end on the middle (in the plane) of an edge running mostly toward the pole meets it, at + // each tolerance: only the edge's ends and the other line's far end dangle. + for ((x0, y0), (x1, y1), tol) in [ + ((0.0, 86.0), (0.0069, 89.9), 1.0), + ((0.0, 86.0), (0.0069, 89.9), 10.0), + ((0.0, 85.1), (0.01, 89.99), 10.0), + ((0.0, 85.1), (0.05, 89.99), 10.0), + ((0.0, 85.06), (0.02, 89.9999), 0.5), + ] { + let (a, b) = (rustgis_geom::lonlat_to_plane(Coord::new(x0, y0)), rustgis_geom::lonlat_to_plane(Coord::new(x1, y1))); + let mid = rustgis_geom::plane_to_lonlat(Coord::new((a.x + b.x) / 2.0, (a.y + b.y) / 2.0)); + let edge = [(x0, y0), (x1, y1)]; + assert_eq!(dangling(&[&edge, &[(mid.x, mid.y), (mid.x + 0.01, mid.y)]], tol), 3, "{edge:?} at {tol} m"); + } +} + +#[test] +fn a_graticule_near_a_pole_is_checked_in_full() { + // Meridians every 0.25° from 80°N up to 86°N, ending on that parallel, drawn as one edge right + // round: only their 1,440 ends at 80°N dangle (the parallel's own ends meet across ±180°), and + // the check finishes, since only the stretch of the long edge within each end's reach is + // measured. + let mut lines: Vec> = (0..1_440) + .map(|i| { + let x = -180.0 + 0.25 * f64::from(i); + vec![(x, 80.0), (x, 86.0)] + }) + .collect(); + lines.push(vec![(-180.0, 86.0), (180.0, 86.0)]); + let refs: Vec<&[(f64, f64)]> = lines.iter().map(Vec::as_slice).collect(); + assert_eq!(dangling(&refs, 1.0), 1_440); +} + +#[test] +fn line_ends_near_a_pole_meet_north_south_and_across_180() { + // Lines end to end along the meridian at 10°E around 89°N: 5 m apart they meet, 30 m apart they + // don't (north–south the plane stretches by sec 85.05° beyond it, not by sec(latitude)). + let gap = degrees; + for (gap_m, want) in [(5.0, 2), (30.0, 4)] { + let a = [(10.0, 88.9), (10.0, 89.0)]; + let b = [(10.0, 89.0 + gap(gap_m)), (10.0, 89.1)]; + assert_eq!(dangling(&[&a, &b], 10.0), want, "{gap_m} m"); + } + // Lines along 89.9°N meeting across ±180°, 0.02° (3.9 m) apart: only their outer ends dangle. + assert_eq!(dangling(&[&[(170.0, 89.9), (179.99, 89.9)], &[(-179.99, 89.9), (-170.0, 89.9)]], 10.0), 2); + // One line all the way around 89.9°N ends where it starts: no dangle. + assert_eq!(dangling(&[&[(-180.0, 89.9), (180.0, 89.9)]], 10.0), 0); + // A line curling back onto itself at 88°S, like a 6: only its first end dangles. + assert_eq!(dangling(&[&[(0.0, -88.1), (0.0, -88.0), (1.0, -88.0), (1.0, -88.05), (0.0, -88.05)]], 10.0), 1); +} + +#[test] +fn an_edge_crossing_85_degrees_is_straight_in_the_plane() { + // A mask whose edge runs from (0°, 80°N) to (10°E, 89°N), across 85.05°: the tools follow it + // straight in the plane, so a line along 84.5°N is cut at 3.99°E. (Web Mercator, clamping the + // edge's far end onto 85.05°, cut it at 8.50°E.) + let c = Coord::new; + let mut line = Dataset::new("line", GeometryType::Polyline, vec![Field::new("k", FieldType::Text)]); + line.push(Some(Geometry::Line(vec![vec![c(0.0, 84.5), c(10.0, 84.5)]])), vec![Value::Text("a".into())]); + let mask = Geometry::Polygon(vec![vec![vec![c(0.0, 80.0), c(10.0, 89.0), c(10.0, 80.0), c(0.0, 80.0)]]]); + let m = Map(HashMap::from([("line".to_string(), line), layer("mask", vec![mask])])); + let clipped = run(&m, "clip", serde_json::json!({"input": "line", "clip": "mask"})); + let e = clipped.features.first().and_then(|f| f.geometry.as_ref()).map(|g| g.extent()).unwrap(); + assert!((e.min_x - 3.986).abs() < 0.002 && (e.max_x - 10.0).abs() < 1e-6, "{e:?}"); +} + +#[test] +fn simplify_near_a_pole_drops_nothing_farther_than_its_tolerance() { + // Beyond 85.05° the plane stretches east–west by sec(latitude) but north–south by a constant + // sec(85.05°), so a tolerance in meters can't be right both ways: for a feature wholly beyond + // 85°, Simplify converts it (at the feature's label point, as everywhere) with a scale no + // larger than the smaller of the two (sec 85°), and nothing it drops is farther than the + // tolerance on the ground, either way. + let c = Coord::new; + // A box from 0° to 10°E and 89°S to 86°S. Its northern edge zigzags 0.05° (5.6 km) and 0.02° + // (2.2 km) north and south; its western edge 1° (3.4 to 7.3 km) and 0.3° (0.7 to 1.9 km) east + // and west. + let mut ring = Vec::new(); + for i in 0..=20 { + let bump = match i % 4 { + 1 => 0.05, + 3 => -0.02, + _ => 0.0, + }; + ring.push(c(10.0 - 0.5 * f64::from(i), -86.0 + bump)); + } + for j in 1..12 { + let bump = match j % 4 { + 1 => 1.0, + 3 => -0.3, + _ => 0.0, + }; + ring.push(c(bump, -86.0 - 0.25 * f64::from(j))); + } + ring.extend([c(0.0, -89.0), c(10.0, -89.0), c(10.0, -86.0)]); + let original = ring.clone(); + let m = Map(HashMap::from([layer("box", vec![Geometry::Polygon(vec![vec![ring]])])])); + let tolerance_km = 3.0; + let simplified = run(&m, "simplify", serde_json::json!({"input": "box", "tolerance": tolerance_km})); + let kept: Vec = match simplified.features.first().and_then(|f| f.geometry.as_ref()) { + Some(Geometry::Polygon(parts)) => parts.first().and_then(|rings| rings.first()).cloned().unwrap_or_default(), + other => panic!("expected a polygon: {other:?}"), + }; + assert!(kept.len() > 4 && kept.len() < original.len(), "{} of {} kept", kept.len(), original.len()); + // How far a point is from what's kept, on the ground: each kept edge followed as the plane + // draws it, in 2,000 steps (under 170 m each here: hence 100 m allowed past the tolerance). + let (to_plane, back) = (rustgis_geom::lonlat_to_plane, rustgis_geom::plane_to_lonlat); + let from_kept = |v: Coord| { + kept.windows(2) + .flat_map(|w| { + let (a, b) = (to_plane(w[0]), to_plane(w[1])); + let steps = 2_000; + (0..=steps).map(move |k| { + let t = f64::from(k) / f64::from(steps); + back(Coord::new(a.x + (b.x - a.x) * t, a.y + (b.y - a.y) * t)) + }) + }) + .map(|q| rustgis_geom::wgs84_distance_m(v, q)) + .fold(f64::INFINITY, f64::min) + }; + for v in original { + let d = from_kept(v); + assert!(d <= tolerance_km * 1000.0 + 100.0, "{v:?} is {d:.0} m from the simplified box"); + } +} + +#[test] +fn the_plane_keeps_longitudes_past_180_as_they_are() { + // A polygon across the antimeridian, its longitudes running on from 179.6° to 180.2° (as moving + // a feature across it leaves them): its convex hull is 0.6° wide, not folded around the world, + // at 10°N as at 86°S, beyond 85.05°. + let c = Coord::new; + for lat in [10.0, -86.0] { + let ring = vec![c(179.6, lat), c(180.2, lat), c(180.2, lat + 0.4), c(179.6, lat + 0.4), c(179.6, lat)]; + let m = Map(HashMap::from([layer("across", vec![Geometry::Polygon(vec![vec![ring]])])])); + let hull = run(&m, "convexHull", serde_json::json!({"input": "across"})); + let e = hull.features.first().and_then(|f| f.geometry.as_ref()).map(|g| g.extent()).unwrap(); + assert!((e.max_x - e.min_x - 0.6).abs() < 1e-9 && (e.min_x - 179.6).abs() < 1e-9, "{lat}: {e:?}"); + assert!((e.min_y - lat).abs() < 1e-7 && (e.max_y - lat - 0.4).abs() < 1e-7, "{lat}: {e:?}"); + } +} + +#[test] +fn a_radius_near_a_pole_reaches_as_far_as_it_says() { + let c = Coord::new; + let mut pts = Dataset::new("pts", GeometryType::Point, vec![]); + // 33.5 km and 55.8 km north of 89°N, and 447 km away at 85°N, on the same meridian. + for y in [89.3, 89.5, 85.0] { + pts.push(Some(Geometry::Point(c(10.0, y))), vec![]); + } + let found = rustgis_analysis::network::within(&pts, &[0, 1, 2], c(10.0, 89.0), 50_000.0); + assert_eq!(found.iter().map(|(r, _)| *r).collect::>(), [0], "{found:?}"); + let d = found.first().map(|(_, d)| *d).unwrap_or_default(); + assert!((d - 33_507.0).abs() < 50.0, "{d}"); + + // 120 km around the South Pole reaches 89°S (111.7 km) and not 88.9°S (122.9 km). + let mut pts = Dataset::new("pts", GeometryType::Point, vec![]); + for (x, y) in [(0.0, -89.0), (10.0, -88.9)] { + pts.push(Some(Geometry::Point(c(x, y))), vec![]); + } + let found = rustgis_analysis::network::within(&pts, &[0, 1], c(0.0, -90.0), 120_000.0); + assert_eq!(found.iter().map(|(r, _)| *r).collect::>(), [0], "{found:?}"); +} + +#[test] +fn a_radius_reaches_polar_shapes_along_their_edges() { + let c = Coord::new; + let ground = rustgis_geom::wgs84_distance_m; + let shapes = |g: Geometry| { + let mut ds = Dataset::new("shapes", GeometryType::Polygon, vec![]); + ds.push(Some(g), vec![]); + ds + }; + let reach = |ds: &Dataset, at: Coord, r: f64| rustgis_analysis::network::within(ds, &[0], at, r).first().map(|(_, d)| *d); + // Inside a box from 84°N to 89°N: no distance at all. From a center that's no place, nothing. + assert_eq!(reach(&shapes(square(0.0, 84.0, 90.0, 89.0)), c(45.0, 85.0), 50_000.0), Some(0.0)); + for at in [c(45.0, 1000.0), c(45.0, f64::NAN), c(f64::INFINITY, 85.0)] { + assert_eq!(reach(&shapes(square(0.0, 84.0, 90.0, 89.0)), at, 50_000_000.0), None, "{at:?}"); + } + // A cap written from 0° to 360° holds a center at 10°W, inside it on the ground, however its + // longitude is written. A radius that's no distance reaches nothing. + let cap360 = shapes(Geometry::Polygon(vec![vec![vec![c(0.0, 85.0), c(360.0, 85.0), c(360.0, 90.0), c(0.0, 90.0), c(0.0, 85.0)]]])); + for x in [-10.0, 710.0, -370.0] { + assert_eq!(reach(&cap360, c(x, 88.0), 100_000.0), Some(0.0), "{x}"); + } + for r in [f64::NAN, -1.0] { + assert!(rustgis_analysis::network::within(&cap360, &[0], c(-10.0, 88.0), r).is_empty(), "{r}"); + } + // A polygon whose ring is one point, 10 km off, written once or closed: as far as that point. + for ring in [vec![c(10.0, 86.0)], vec![c(10.0, 86.0), c(10.0, 86.0)]] { + let d = reach(&shapes(Geometry::Polygon(vec![vec![ring]])), c(10.0, 86.09), 50_000.0); + assert!(d.is_some_and(|d| (d - ground(c(10.0, 86.09), c(10.0, 86.0))).abs() < 10.0), "{d:?}"); + } + // A cap cut at ±180°, its edge along 85.5°N half a degree north: that far, not as far as the + // pole. + let cap = Geometry::Polygon(vec![vec![vec![c(-180.0, 85.5), c(180.0, 85.5), c(180.0, 90.0), c(-180.0, 90.0), c(-180.0, 85.5)]]]); + let d = reach(&shapes(cap), c(0.0, 85.0), 100_000.0); + assert!(d.is_some_and(|d| (d - ground(c(0.0, 85.0), c(0.0, 85.5))).abs() < 10.0), "{d:?}"); + // A meridian from pole to pole, a degree west at 45°N: about 79 km off, not 0 by a chord. + let mut line = Dataset::new("meridian", GeometryType::Polyline, vec![]); + line.push(Some(Geometry::Line(vec![vec![c(0.0, -90.0), c(0.0, 90.0)]])), vec![]); + let d = reach(&line, c(1.0, 45.0), 100_000.0); + assert!(d.is_some_and(|d| (d - ground(c(1.0, 45.0), c(0.0, 45.0))).abs() < 100.0), "{d:?}"); + // From 80°N, 1,600 km reaches a point beyond 85.05° on the far side of the pole, 1,563 km away + // over it. One within ±85.05°, from a center within it, is measured in local Web Mercator as + // before, which doesn't go over the pole: 2,300 km from 80°N doesn't reach 80°N across it. + let mut pts = Dataset::new("pts", GeometryType::Point, vec![]); + pts.push(Some(Geometry::Point(c(180.0, 86.0))), vec![]); + let d = reach(&pts, c(0.0, 80.0), 1_600_000.0); + assert!(d.is_some_and(|d| (d - ground(c(0.0, 80.0), c(180.0, 86.0))).abs() < 100.0), "{d:?}"); + let mut pts = Dataset::new("pts", GeometryType::Point, vec![]); + pts.push(Some(Geometry::Point(c(180.0, 80.0))), vec![]); + assert_eq!(reach(&pts, c(0.0, 80.0), 2_300_000.0), None); + // A line along 180° from 88.1°N to 89.4°S, through the antipode (0°, 89°S) of a center at 89°N, + // where the plane at the center tears: 400 km reaches it 2.9° away over the pole (324 km), not + // at a chord across the tear; 200 km doesn't. + let mut line = Dataset::new("through", GeometryType::Polyline, vec![]); + line.push(Some(Geometry::Line(vec![vec![c(180.0, 88.1), c(180.0, -89.4)]])), vec![]); + let d = reach(&line, c(0.0, 89.0), 400_000.0); + assert!(d.is_some_and(|d| (d - ground(c(0.0, 89.0), c(180.0, 88.1))).abs() < 100.0), "{d:?}"); + assert_eq!(reach(&line, c(0.0, 89.0), 200_000.0), None); + // A line from near the center (180°, 89°N; 55 km off at its nearest) down to 89°S and along it + // past the center's antipode (0°, 89°S), where the plane at the center draws neighboring steps + // on opposite sides of it, their chord 40 km from the center: not within 50 km. + let mut line = Dataset::new("far", GeometryType::Polyline, vec![]); + // (Steps from 150.25°E fall either side of 0°, not on it.) + line.push(Some(Geometry::Line(vec![vec![c(150.25, 89.5), c(150.25, -89.0), c(-170.0, -89.0)]])), vec![]); + assert_eq!(reach(&line, c(180.0, 89.0), 50_000.0), None); + let d = reach(&line, c(180.0, 89.0), 60_000.0); + assert!(d.is_some_and(|d| (52_000.0..58_000.0).contains(&d)), "{d:?}"); +} + +#[test] +fn a_radius_follows_hostile_lines_near_a_pole_in_bounded_steps() { + // 200 lines, each going 348 times nearly all the way around along 89.95°N, 5.6 km from a center + // at 89.9°N: about 50 million half-degree steps, far past the search's budget. The edges + // followed before it runs out find the first line where it is; every line still counts, the + // rest at their boxes' distance (5.5 km), never more than they are. + let c = Coord::new; + let path: Vec = (0..=348).map(|i| c(if i % 2 == 0 { -179.5 } else { 179.5 }, 89.95)).collect(); + let mut lines = Dataset::new("hostile", GeometryType::Polyline, vec![]); + for _ in 0..200 { + lines.push(Some(Geometry::Line(vec![path.clone()])), vec![]); + } + let rows: Vec = (0..200).collect(); + let found = rustgis_analysis::network::within(&lines, &rows, c(0.0, 89.9), 10_000.0); + assert_eq!(found.len(), 200); + assert!(found.iter().all(|(_, d)| (5_500.0..5_600.0).contains(d)), "{:?}", found.first()); +} + +#[test] +fn a_larger_radius_never_loses_a_feature() { + // Each shape is measured one way whatever the radius. Within ±85.05°, in local Web Mercator: a + // point at 75°N from Paris, whether the radius reaches 85.05° or not, and one at 82.95°N from + // 84°N the same. + let c = Coord::new; + let found = |at: Coord, p: Coord, r: f64| { + let mut pts = Dataset::new("pts", GeometryType::Point, vec![]); + pts.push(Some(Geometry::Point(p)), vec![]); + rustgis_analysis::network::within(&pts, &[0], at, r).first().map(|(_, d)| *d) + }; + let paris = c(2.35, 48.86); + assert_eq!(found(paris, c(2.35, 75.0), 4_000_000.0), found(paris, c(2.35, 75.0), 4_100_000.0)); + let (at, p) = (c(0.0, 84.0), c(3.9277, 82.9505)); + let (near, far) = (found(at, p, 117_000.0), found(at, p, 117_020.0)); + assert!(near.is_some() && near == far, "{near:?} {far:?}"); + // On the ground: a line along 89.99°N, 10 km from a center at 89.9°N, found at 20 km, where 400 + // lines around 89.6°N (33 km off) are out of reach, and at 50 km, where following those lines + // first leaves too few of the search's steps for it, so it counts as near as its box (9.95 km): + // never farther than it is. + let at = c(0.0, 89.9); + let mut ds = Dataset::new("polar", GeometryType::Polyline, vec![]); + for _ in 0..400 { + ds.push(Some(Geometry::Line(vec![vec![c(-179.5, 89.6), c(179.5, 89.6)]])), vec![]); + } + ds.push(Some(Geometry::Line(vec![vec![c(-150.0, 89.99), c(150.0, 89.99)]])), vec![]); + let rows: Vec = (0..=400).collect(); + let point = |r: f64| rustgis_analysis::network::within(&ds, &rows, at, r).into_iter().find(|(row, _)| *row == 400).map(|(_, d)| d); + let (near, far) = (point(20_000.0), point(50_000.0)); + let fell_back = far.is_some_and(|f| (9_940.0..9_960.0).contains(&f) && near.is_some_and(|n| f < n)); + assert!(near.is_some_and(|d| (9_900.0..10_200.0).contains(&d)) && fell_back, "{near:?} {far:?}"); +} + +#[test] +fn features_reaching_the_same_pole_meet_there() { + // Wedges from 89°S to the South Pole at 0°E and 90°E: apart in the plane, where the pole is a + // line, but one place on the ground. One at 0°E reaching the North Pole meets neither. + let wedge = |x: f64, s: f64| square(x, s * 90.0, x + 10.0, s * 89.0); + let m = Map(HashMap::from([ + layer("south0", vec![wedge(0.0, -1.0)]), + layer("south90", vec![wedge(90.0, -1.0)]), + layer("north0", vec![wedge(0.0, 1.0)]), + ])); + let selected = |input: &str, select: &str| { + let r = rustgis_analysis::run("selectByLocation", &m, &serde_json::json!({"input": input, "select": select})).unwrap(); + let Output::Selection { rows, .. } = r.output else { panic!("expected a selection") }; + rows + }; + assert_eq!((selected("south0", "south90"), selected("north0", "south90")), (vec![0], vec![])); + let r = rustgis_analysis::run("selectByLocation", &m, &serde_json::json!({"input": "south0", "select": "south90", "invert": "Yes"})) + .unwrap(); + let Output::Selection { rows, .. } = r.output else { panic!("expected a selection") }; + assert!(rows.is_empty(), "{rows:?}"); + let joined = |target: &str, join: &str, how: &str| { + let mut params = serde_json::json!({"target": target, "join": join, "match": how}); + if how == "Closest" { + params["distanceField"] = serde_json::json!("gap"); + } + run(&m, "spatialJoin", params) + }; + let sj = joined("south0", "south90", "Intersects"); + assert_eq!(sj.value(0, sj.field_index("join_count").unwrap()), &Value::Int(1)); + let sj = joined("north0", "south90", "Intersects"); + assert_eq!(sj.value(0, sj.field_index("join_count").unwrap()), &Value::Int(0)); + let sj = joined("south0", "south90", "Closest"); + assert_eq!(sj.value(0, sj.field_index("gap").unwrap()), &Value::Float(0.0)); +} + +#[test] +fn a_clipped_wedge_still_meets_others_at_the_pole() { + // A wedge from 0°E to 90°E, clipped, still reaches the pole, where it meets one at 120°E: the + // overlay's grid moves its pole vertices, but not out of the pole's tolerance. + let m = Map(HashMap::from([ + layer("wedge", vec![square(0.0, -90.0, 90.0, -89.0)]), + layer("mask", vec![square(-180.0, -90.0, 180.0, -80.0)]), + ])); + let clipped = run(&m, "clip", serde_json::json!({"input": "wedge", "clip": "mask"})); + let off = latitudes(&clipped).1.iter().map(|y| (y + 90.0).abs()).fold(f64::INFINITY, f64::min); + let m = Map(HashMap::from([("clipped".to_string(), clipped), layer("other", vec![square(120.0, -90.0, 130.0, -89.0)])])); + let r = rustgis_analysis::run("selectByLocation", &m, &serde_json::json!({"input": "clipped", "select": "other"})).unwrap(); + let Output::Selection { rows, .. } = r.output else { panic!("expected a selection") }; + assert_eq!(rows, [0], "pole vertices {off:e}° off"); +} + +#[test] +fn a_spatial_join_finds_every_feature_at_the_pole_through_its_index() { + // Twenty wedges reaching the South Pole all around it, and five squares far off: enough join + // features (over 16) for Spatial Join to look them up in its grid, where a wedge's extent + // meets only its neighbors'. + let wedges = (0..20).map(|i| square(-180.0 + 18.0 * f64::from(i), -90.0, -172.0 + 18.0 * f64::from(i), -89.0)); + let far = (0..5).map(|i| square(f64::from(i), 0.0, f64::from(i) + 0.5, 0.5)); + let m = Map(HashMap::from([layer("target", vec![square(0.0, -90.0, 10.0, -89.0)]), layer("join", wedges.chain(far).collect())])); + let sj = run(&m, "spatialJoin", serde_json::json!({"target": "target", "join": "join", "match": "Intersects"})); + assert_eq!(sj.value(0, sj.field_index("join_count").unwrap()), &Value::Int(20)); +} + +#[test] +fn a_point_at_a_pole_is_found_through_each_index() { + // The South Pole points of the test below, with twenty points and twenty squares far off: more + // than 16 features on each side, so the tools look them up in their grids, where the pole + // points' extents don't meet the wedge's. + let c = Coord::new; + let mut pts = Dataset::new("pts", GeometryType::Point, vec![Field::new("k", FieldType::Text)]); + let polar = [c(0.0, -90.0), c(120.0, -90.0), c(0.0, -90.0 + 5e-8), c(45.0, -90.0 + 5e-8)]; + for p in polar.into_iter().chain((0..20).map(|i| c(100.0 + f64::from(i), 10.0))) { + pts.push(Some(Geometry::Point(p)), vec![Value::Text("a".into())]); + } + let far = (0..20).map(|i| square(f64::from(i), 0.0, f64::from(i) + 0.5, 0.5)); + let (_, polys) = layer("polys", std::iter::once(square(0.0, -90.0, 90.0, -89.0)).chain(far).collect()); + let m = Map(HashMap::from([("pts".to_string(), pts), ("polys".to_string(), polys)])); + for (tool, params, want) in [ + ("clip", serde_json::json!({"input": "pts", "clip": "polys"}), 4), + ("intersect", serde_json::json!({"input": "pts", "overlay": "polys"}), 4), + ("erase", serde_json::json!({"input": "pts", "erase": "polys"}), 20), + ] { + assert_eq!(run(&m, tool, params).features.len(), want, "{tool}"); + } + let counted = run(&m, "countPoints", serde_json::json!({"polygons": "polys", "points": "pts"})); + assert_eq!(counted.value(0, counted.field_index("point_count").unwrap()), &Value::Int(4)); + let split = rustgis_analysis::run("split", &m, &serde_json::json!({"input": "pts", "split": "polys"})).unwrap(); + let Output::Layers(layers) = split.output else { panic!("expected layers") }; + let [rustgis_analysis::OutputLayer::Vector(part)] = layers.as_slice() else { panic!("expected one layer: {layers:?}") }; + assert_eq!(part.features.len(), 4); + let joined = run(&m, "spatialJoin", serde_json::json!({"target": "pts", "join": "polys", "match": "Have their center in"})); + let count = joined.field_index("join_count").unwrap(); + let found: Vec = (0..24).filter(|r| joined.value(*r, count) == &Value::Int(1)).collect(); + assert_eq!(found, [0, 1, 2, 3]); +} + +#[test] +fn a_point_at_a_pole_meets_what_reaches_it_and_lies_in_what_surrounds_it() { + // In the plane a pole is a line, on the ground one place. Points at the South Pole, written at + // 0°E and 120°E and 5e-8° off it, touch a wedge from 0°E to 90°E reaching the pole: Clip and + // Intersect keep them and Summarize Within counts them, as they do points on an edge. Only a + // polygon that surrounds the pole holds them inside, as Must Be Covered By counts; nothing + // reaches the North Pole point. A point 5e-8° off the pole within the wedge, at 45°E, lies + // inside it as anywhere else. + let c = Coord::new; + let mut pts = Dataset::new("pts", GeometryType::Point, vec![Field::new("k", FieldType::Text)]); + for p in [c(0.0, -90.0), c(120.0, -90.0), c(0.0, -90.0 + 5e-8), c(0.0, 90.0), c(45.0, -90.0 + 5e-8)] { + pts.push(Some(Geometry::Point(p)), vec![Value::Text("a".into())]); + } + let (_, wedge) = layer("wedge", vec![square(0.0, -90.0, 90.0, -89.0)]); + let m = Map(HashMap::from([("pts".to_string(), pts.clone()), ("wedge".to_string(), wedge.clone())])); + for (tool, params) in [ + ("clip", serde_json::json!({"input": "pts", "clip": "wedge"})), + ("intersect", serde_json::json!({"input": "pts", "overlay": "wedge"})), + ] { + assert_eq!(latitudes(&run(&m, tool, params)).0, 4, "{tool}"); + } + let counted = run(&m, "countPoints", serde_json::json!({"polygons": "wedge", "points": "pts"})); + assert_eq!(counted.value(0, counted.field_index("point_count").unwrap()), &Value::Int(4)); + // Erase takes them away, leaving the North Pole point; Split puts them in the wedge's layer; + // Spatial Join finds each in the wedge when it has its center there. + let erased = run(&m, "erase", serde_json::json!({"input": "pts", "erase": "wedge"})); + assert_eq!(latitudes(&erased), (1, vec![90.0])); + let split = rustgis_analysis::run("split", &m, &serde_json::json!({"input": "pts", "split": "wedge"})).unwrap(); + let Output::Layers(layers) = split.output else { panic!("expected layers") }; + let [rustgis_analysis::OutputLayer::Vector(part)] = layers.as_slice() else { panic!("expected one layer: {layers:?}") }; + assert_eq!(part.features.len(), 4); + let joined = run(&m, "spatialJoin", serde_json::json!({"target": "pts", "join": "wedge", "match": "Have their center in"})); + let count = joined.field_index("join_count").unwrap(); + let counts: Vec<&Value> = (0..5).map(|r| joined.value(r, count)).collect(); + assert_eq!(counts, [&Value::Int(1), &Value::Int(1), &Value::Int(1), &Value::Int(0), &Value::Int(1)]); + let uncovered = |by: &Dataset| -> Vec { + rustgis_analysis::topology::not_covered(&pts, &[0, 1, 2, 3, 4], by).errors.iter().flat_map(|e| e.rows.clone()).collect() + }; + assert_eq!(uncovered(&wedge), [0, 1, 2, 3]); + // Surrounding the pole: a cap written from −180° to 180° or from 0° to 360°, one in two halves + // (meeting exactly, or 1e-9° apart, as an overlay's rounding can leave them), one in parts + // written in both ways, one with an island far west; not one leaving a gap of 10° or of 1° (its + // pole points lie on its edge), nor one with a hole reaching the pole (the point beside the + // pole, at 0°, lies inside it, as anywhere else), even one narrowing to a corner there, under a + // meter wide in the plane a meter off the pole, unless another polygon fills it; nor halves + // with a gap narrowing to the pole between them (along a meridian or slanting, 0.5° wide at + // 85°S, millimeters a meter off the pole in the plane), or one with a notch narrowing to it; + // nor one whose longitudes can't be taken (one infinite, or so far around that a meter is lost + // in rounding). A cap 110 m across holds it, in halves a seam apart too; a hole near the pole + // that doesn't reach it is reported as if it did, as the check leans to reporting. + let ring = |w: f64, e: f64| vec![c(w, -85.0), c(e, -85.0), c(e, -90.0), c(w, -90.0), c(w, -85.0)]; + let island = vec![c(-100.0, 0.0), c(-95.0, 0.0), c(-95.0, 5.0), c(-100.0, 5.0), c(-100.0, 0.0)]; + let hole = vec![c(-60.0, -90.0), c(-50.0, -90.0), c(-50.0, -89.0), c(-60.0, -89.0), c(-60.0, -90.0)]; + let narrowing = vec![c(-60.0, -89.0), c(-50.0, -89.0), c(-55.0, -90.0), c(-60.0, -89.0)]; + let tapering = vec![c(0.0, -85.0), c(180.0, -85.0), c(180.0, -90.0), c(1.0, -90.0), c(0.0, -89.0), c(0.0, -85.0)]; + let notched = vec![ + c(-180.0, -85.0), + c(180.0, -85.0), + c(180.0, -90.0), + c(2.0, -90.0), + c(0.0, -89.5), + c(-2.0, -90.0), + c(-180.0, -90.0), + c(-180.0, -85.0), + ]; + let infinite = vec![c(-180.0, -85.0), c(f64::INFINITY, -85.0), c(180.0, -90.0), c(-180.0, -90.0), c(-180.0, -85.0)]; + let half_wedge = vec![c(0.5, -85.0), c(180.0, -85.0), c(180.0, -90.0), c(0.0, -90.0), c(0.5, -85.0)]; + let near = vec![c(10.0, -89.999), c(20.0, -89.999), c(20.0, -89.99), c(10.0, -89.99), c(10.0, -89.999)]; + let small = vec![c(-180.0, -89.9995), c(180.0, -89.9995), c(180.0, -90.0), c(-180.0, -90.0), c(-180.0, -89.9995)]; + let slanted_west = vec![c(-180.0, -85.0), c(10.0, -85.0), c(0.0, -90.0), c(-180.0, -90.0), c(-180.0, -85.0)]; + let slanted_east = vec![c(10.5, -85.0), c(180.0, -85.0), c(180.0, -90.0), c(0.0, -90.0), c(10.5, -85.0)]; + let small_west = vec![c(-180.0, -89.9995), c(0.0, -89.9995), c(0.0, -90.0), c(-180.0, -90.0), c(-180.0, -89.9995)]; + let small_east = vec![c(1e-9, -89.9995), c(180.0, -89.9995), c(180.0, -90.0), c(1e-9, -90.0), c(1e-9, -89.9995)]; + let notch_outline = vec![ + c(0.0, -85.0), + c(360.0, -85.0), + c(360.0, -90.0), + c(200.0, -90.0), + c(200.0, -86.0), + c(100.0, -86.0), + c(100.0, -90.0), + c(0.0, -90.0), + c(0.0, -85.0), + ]; + let notch_hole = vec![c(100.0, -90.0), c(200.0, -90.0), c(200.0, -86.5), c(100.0, -86.5), c(100.0, -90.0)]; + let short_of = vec![c(-180.0, -85.0), c(180.0, -85.0), c(180.0, -89.9999995), c(-180.0, -89.9999995), c(-180.0, -85.0)]; + let nearly = vec![ + c(-180.0, -85.0), + c(180.0, -85.0), + c(180.0, -90.0), + c(2.0, -90.0), + c(0.0, -89.9999996), + c(-2.0, -90.0), + c(-180.0, -90.0), + c(-180.0, -85.0), + ]; + let one = |rings: Vec>>| vec![Geometry::Polygon(rings)]; + let each = |rings: Vec>>| rings.into_iter().map(|r| Geometry::Polygon(vec![r])).collect::>(); + for (name, features, want) in [ + ("cap", one(vec![vec![ring(-180.0, 180.0)]]), vec![3]), + ("cap from 0°", one(vec![vec![ring(0.0, 360.0)]]), vec![3]), + ("halves", one(vec![vec![ring(-180.0, 0.0)], vec![ring(0.0, 180.0)]]), vec![3]), + ("halves a seam apart", one(vec![vec![ring(-180.0, 0.0)], vec![ring(1e-9, 180.0)]]), vec![3]), + ("halves overlapping", one(vec![vec![ring(-180.0, 10.0)], vec![ring(-10.0, 180.0)]]), vec![3]), + ("halves as two features", each(vec![vec![ring(-180.0, 0.0)], vec![ring(0.0, 180.0)]]), vec![3]), + ("parts both ways", one(vec![vec![ring(-180.0, -10.0)], vec![ring(0.0, 180.0)], vec![ring(350.0, 360.0)]]), vec![3]), + ("island", one(vec![vec![ring(0.0, 360.0)], vec![island.clone()]]), vec![3]), + ("gap", one(vec![vec![ring(0.0, 350.0)]]), vec![0, 1, 2, 3]), + ("halves a degree apart", each(vec![vec![ring(-180.0, 0.0)], vec![ring(1.0, 180.0)]]), vec![0, 1, 2, 3]), + ("hole", one(vec![vec![ring(-180.0, 180.0), hole.clone()]]), vec![0, 1, 3]), + ("hole narrowing to the pole", one(vec![vec![ring(-180.0, 180.0), narrowing.clone()]]), vec![0, 1, 3]), + ( + "that hole filled", + vec![Geometry::Polygon(vec![vec![ring(-180.0, 180.0), narrowing.clone()]]), Geometry::Polygon(vec![vec![narrowing.clone()]])], + vec![3], + ), + ("halves a gap narrowing to the pole apart", each(vec![vec![ring(-180.0, 0.0)], vec![tapering.clone()]]), vec![0, 1, 2, 3]), + ("a notch to the pole", one(vec![vec![notched.clone()]]), vec![0, 1, 2, 3]), + ("halves a 0.5° wedge apart", each(vec![vec![ring(-180.0, 0.0)], vec![half_wedge.clone()]]), vec![0, 1, 2, 3]), + ("a hole near the pole, not reaching it", one(vec![vec![ring(-180.0, 180.0), near.clone()]]), vec![0, 1, 3]), + ("halves a slanted 0.5° wedge apart", each(vec![vec![slanted_west.clone()], vec![slanted_east.clone()]]), vec![0, 1, 3]), + ("a small cap in halves a seam apart", each(vec![vec![small_west.clone()], vec![small_east.clone()]]), vec![3]), + ("a cap stopping 5e-7° short of the pole", one(vec![vec![short_of.clone()]]), vec![0, 1, 2, 3, 4]), + ("a notch ending 4e-7° short of the pole", one(vec![vec![nearly.clone()]]), vec![0, 1, 2, 3]), + ("a hole beside its outline, in its gap", one(vec![vec![ring(0.0, 350.0), ring(350.0, 360.0)]]), vec![0, 1, 2, 3]), + ("a hole filling a notch in its outline", one(vec![vec![notch_outline.clone(), notch_hole.clone()]]), vec![0, 1, 2, 3]), + ("a cap 110 m across", one(vec![vec![small.clone()]]), vec![3]), + ("a longitude infinite", one(vec![vec![infinite.clone()]]), vec![0, 1, 3]), + ("a turn written 1e15° away", one(vec![vec![ring(1e15, 1e15 + 360.0)]]), vec![0, 1, 2, 3, 4]), + ] { + assert_eq!(uncovered(&layer("cap", features).1), want, "{name}"); + } +} diff --git a/crates/analysis/tests/reference.rs b/crates/analysis/tests/reference.rs index 5120f44ae..76f206801 100644 --- a/crates/analysis/tests/reference.rs +++ b/crates/analysis/tests/reference.rs @@ -104,7 +104,7 @@ fn polygons(gs: &[Geometry]) -> MultiPolygon { MultiPolygon::new(gs.iter().filter_map(to_multipolygon).flat_map(|m| m.0).collect()) } -/// The lines among `gs`, together, in Web Mercator metres. +/// The lines among `gs`, together, in the analysis plane. fn lines(gs: &[Geometry]) -> MultiLineString { MultiLineString::new(gs.iter().filter_map(to_multiline).flat_map(|m| m.0).collect()) } diff --git a/crates/analysis/tests/tabulate_intersection.rs b/crates/analysis/tests/tabulate_intersection.rs index 0325934e1..962b43395 100644 --- a/crates/analysis/tests/tabulate_intersection.rs +++ b/crates/analysis/tests/tabulate_intersection.rs @@ -255,10 +255,10 @@ fn overlapping_zones_still_say_what_lies_outside() { assert_eq!(got.iter().map(|(_, n)| n[2]).collect::>(), [50.0, 50.0]); } -/// A class feature reaching beyond the map's latitude limit, and features lacking values in two sum -/// fields: one notice for each kind. +/// Features lacking values in two sum fields: one notice for each field, and none for a zone and a +/// class feature reaching past 85.05°, which are measured there like anywhere else (#1033). #[test] -fn one_notice_per_kind_for_missing_values_and_the_latitude_limit() { +fn one_notice_per_field_for_missing_values() { let polar = |x0: f64, x1: f64, y0: f64, y1: f64| { let ring = [(x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0)].map(|(x, y)| Coord::new(x, y)).to_vec(); Some(Geometry::Polygon(vec![vec![ring]])) @@ -285,12 +285,55 @@ fn one_notice_per_kind_for_missing_values_and_the_latitude_limit() { "Tabulate Intersection: 1 row", "1 zone overlaps no class feature, so it has no rows", "class features with no value in a sum field are left out of its shares: 2 class features in 'pop', 1 class feature in 'income'", - "1 zone feature and 1 class feature reach beyond 85.05° latitude, where the map ends: parts there are taken as lying along it, with no area", + "1 zone feature and 1 class feature have an edge beyond 85.05° running east or west: the analysis plane takes it straight in longitude and latitude beyond 85.05°, and their percentages and shares follow the plane, while area_km2 takes the shortest way on the ground, which near a pole can differ widely", "1 class feature lies partly or wholly outside the zones: its value outside isn't shared out", ] ); } +/// Polygons reaching past 85.05°, where the overlay goes on to the poles (#1033): a class feature +/// from 84°N to 89°N shares out 91.42% of its value to a zone from 80°N to 88°N (on the sphere, +/// the part from 84° to 88° over the whole), not all of it as when its part beyond 85.05° had no +/// area; and one wholly beyond, from 86°N to 90°N, has its row in a zone from 80°N to the pole, +/// 16.03% of it, rather than none. +#[test] +fn polygons_past_85_degrees_are_measured_to_the_pole() { + let polar = |x0: f64, x1: f64, y0: f64, y1: f64| { + let ring = [(x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0)].map(|(x, y)| Coord::new(x, y)).to_vec(); + Some(Geometry::Polygon(vec![vec![ring]])) + }; + let mut z = Dataset::new("zones", GeometryType::Polygon, vec![Field::new("zone", FieldType::Text)]); + z.push_next(polar(0.0, 10.0, 80.0, 88.0), vec![Value::Text("A".into())]); + z.push_next(polar(20.0, 30.0, 80.0, 90.0), vec![Value::Text("B".into())]); + let mut c = Dataset::new("classes", GeometryType::Polygon, vec![Field::new("pop", FieldType::Double)]); + c.push_next(polar(0.0, 10.0, 84.0, 89.0), vec![Value::Float(1000.0)]); + c.push_next(polar(20.0, 30.0, 86.0, 90.0), vec![Value::Float(500.0)]); + let l = layers(z, c); + let params = json!({"zones": "zones", "zoneFields": "zone", "classes": "classes", "sumFields": ["pop"]}); + let res = rustgis_analysis::run("tabulateIntersection", &l, ¶ms).unwrap(); + // Each has an edge beyond 85.05° running east or west, which the map and the ground take + // differently: the tool says so. + let said = "2 zone features and 2 class features have an edge beyond 85.05° running east or west"; + assert!(res.messages.iter().any(|m| m.starts_with(said)), "{:?}", res.messages); + let got = rows(&table(res)); + // zone, then area_km2, percentage, pop. + assert_eq!(got.len(), 2, "{got:?}"); + let (a, b) = (&got[0], &got[1]); + assert_eq!((a.0[0].as_str(), b.0[0].as_str()), ("A", "B")); + // The parts' areas on WGS 84, as polygons of geodesic edges (GeographicLib 2.1 gives 34,625.271 + // and 17,323.559 km²), not 0 for the stretch beyond 85.05°. + assert!((a.1[0] - 34_625.271).abs() < 0.01 && (b.1[0] - 17_323.559).abs() < 0.01, "{a:?} {b:?}"); + assert!((a.1[1] - 33.388).abs() < 0.01 && (a.1[2] - 914.216).abs() < 0.1, "{a:?}"); + assert!((b.1[1] - 16.034).abs() < 0.01 && (b.1[2] - 500.0).abs() < 1e-6, "{b:?}"); + // Within ±85.05°, or beyond it only along meridians, nothing of the kind. + let mut z = Dataset::new("zones", GeometryType::Polygon, vec![Field::new("zone", FieldType::Text)]); + z.push_next(polar(0.0, 10.0, 70.0, 85.0), vec![Value::Text("A".into())]); + let mut c = Dataset::new("classes", GeometryType::Polygon, vec![Field::new("pop", FieldType::Double)]); + c.push_next(polar(0.0, 10.0, 75.0, 85.0), vec![Value::Float(1000.0)]); + let res = rustgis_analysis::run("tabulateIntersection", &layers(z, c), ¶ms).unwrap(); + assert!(res.messages.iter().all(|m| !m.contains("beyond 85.05°")), "{:?}", res.messages); +} + /// Features of about 10 m in zones of half a degree: the overlay rounds their parts to a grid /// scaled to the zones, so a feature wholly inside a zone still counts whole, and one across the /// border between two zones shares its value out in full within that rounding, with no notice. diff --git a/crates/engine/src/cmd/design.rs b/crates/engine/src/cmd/design.rs index a6dea3609..8b0967f5a 100644 --- a/crates/engine/src/cmd/design.rs +++ b/crates/engine/src/cmd/design.rs @@ -1289,8 +1289,11 @@ fn topology_validate(s: &mut Session, p: &Json) -> CmdResult { use rustgis_analysis::topology as t; let tol = f(p, "toleranceM").unwrap_or(1.0).clamp(0.0, 10_000.0); let mut errors = Vec::new(); - // Rules whose check stopped early (an error or work limit): their counts are partial. + // Rules whose check stopped early (an error or work limit), or measured some part more coarsely + // than it's meant to (a line winding round a pole many times): their counts may be off. let mut incomplete = Vec::new(); + // Polygons a check had to leave out, by rule. + let mut skipped = Vec::new(); for rule in s.map().editing.topology.clone() { let Some(l) = s.map().layer_by_name(rule.layer()) else { continue }; let rows: Vec = l.visible_rows().iter().enumerate().filter(|(_, v)| **v).map(|(i, _)| i).collect(); @@ -1307,6 +1310,14 @@ fn topology_validate(s: &mut Session, p: &Json) -> CmdResult { if !found.complete { incomplete.push(rule.label().to_string()); } + if found.skipped > 0 { + // The features left out are the covering layer's. + let of = match &rule { + TopologyRule::MustBeCoveredBy { by, .. } => by.as_str(), + _ => l.name.as_str(), + }; + skipped.push(format!("{} on '{}': {} feature(s) of '{}'", rule.label(), l.name, found.skipped, of)); + } errors.extend(found.errors.into_iter().map(|e| crate::TopologyIssue { rule: rule.label().to_string(), layer: l.name.clone(), @@ -1315,6 +1326,18 @@ fn topology_validate(s: &mut Session, p: &Json) -> CmdResult { lat: e.at.y, })); } + if !incomplete.is_empty() { + s.warn(format!( + "Checks that stopped early (an error or work limit) or measured a line round a pole coarsely, so their counts may be off: {}", + incomplete.join(", ") + )); + } + if !skipped.is_empty() { + s.warn(format!( + "Polygons left out of the check at a pole (a coordinate missing or infinite, or a million turns around), so its result there may be wrong: {}", + skipped.join("; ") + )); + } let n = errors.len(); let out: Vec = errors.iter().take(1000).map(|e| json!({"rule": e.rule, "layer": e.layer, "rows": e.rows, "lon": e.lon, "lat": e.lat})).collect(); diff --git a/crates/engine/src/cmd/select.rs b/crates/engine/src/cmd/select.rs index c1f80a88f..a5764bd5a 100644 --- a/crates/engine/src/cmd/select.rs +++ b/crates/engine/src/cmd/select.rs @@ -239,17 +239,38 @@ fn by_radius(s: &mut Session, p: &Value) -> CmdResult { flat.validate_selection_path(&boundary)?; let scale = s.scale(); let mut hits = Vec::new(); + let mut bounded = Vec::new(); for layer in s.map_mut().selectable_layers(scale) { let drawn = layer.drawn_rows(scale); let rows: Vec = drawn.iter().enumerate().filter_map(|(row, shown)| shown.then_some(row)).collect(); - let rows: Vec = network::within(&layer.dataset, &rows, center, radius).into_iter().map(|(row, _)| row).collect(); + let (near, n) = network::within_counted(&layer.dataset, &rows, center, radius); + if n > 0 { + bounded.push((layer.name.clone(), n)); + } + let rows: Vec = near.into_iter().map(|(row, _)| row).collect(); if !rows.is_empty() { hits.push((layer.id, rows)); } } + warn_bounded(s, &bounded); apply_hits(s, hits, m) } +/// One warning for the features of each layer a radius search ([`network::within_counted`]) +/// counted at a bound on their distance near a pole, with their counts (#1033); none when there are +/// none. +pub(crate) fn warn_bounded(s: &mut Session, bounded: &[(String, usize)]) { + if bounded.is_empty() { + return; + } + let total: usize = bounded.iter().map(|(_, n)| n).sum(); + let layers: Vec = bounded.iter().map(|(layer, n)| format!("{layer} {n}")).collect(); + s.warn(format!( + "{total} feature(s) near a pole ({}) couldn't be measured in full (past a search's 250,000 steps a layer, or along an edge over 512° long) and were counted at the least distance they could be: they may lie farther, even beyond the radius", + layers.join(", ") + )); +} + /// The hit count describes the query, regardless of whether its rows are added or removed. fn apply_hits(s: &mut Session, hits: Vec<(u32, Vec)>, m: Mode) -> CmdResult { if m == Mode::New { diff --git a/crates/engine/src/cmd/workflows.rs b/crates/engine/src/cmd/workflows.rs index d0331a5da..aa3788176 100644 --- a/crates/engine/src/cmd/workflows.rs +++ b/crates/engine/src/cmd/workflows.rs @@ -723,17 +723,24 @@ fn route(s: &mut Session, p: &Value) -> CmdResult { fn neighborhood(s: &mut Session, p: &Value) -> CmdResult { let at = Coord::new(f(p, "lon").ok_or("lon is required")?, f(p, "lat").ok_or("lat is required")?); + if !(at.x.is_finite() && at.y.abs() <= 90.0) { + return Err("lon must be a finite number and lat within -90 to 90".into()); + } let radius = f(p, "radius").filter(|r| *r > 0.0).ok_or("radius is required (metres)")?.min(5e6); let only: Option> = p.get("layers").and_then(Value::as_array).map(|a| a.iter().filter_map(|v| v.as_str().map(str::to_ascii_lowercase)).collect()); let mut out = Vec::new(); let mut found: Vec<(u32, Vec)> = Vec::new(); + let mut bounded = Vec::new(); for l in s.map().layers.iter().filter(|l| s.map().shows(l) && !l.is_raster() && !l.graphics) { if only.as_ref().is_some_and(|o| !o.contains(&l.name.to_ascii_lowercase())) { continue; } let rows: Vec = l.visible_rows().iter().enumerate().filter(|(_, v)| **v).map(|(i, _)| i).collect(); - let near = network::within(&l.dataset, &rows, at, radius); + let (near, n) = network::within_counted(&l.dataset, &rows, at, radius); + if n > 0 { + bounded.push((l.name.clone(), n)); + } if near.is_empty() { continue; } @@ -749,6 +756,7 @@ fn neighborhood(s: &mut Session, p: &Value) -> CmdResult { out.push(json!({"layer": l.name, "count": near.len(), "nearestKm": near.first().map(|x| x.1 / 1000.0), "sums": sums})); found.push((l.id, near.into_iter().map(|x| x.0).collect())); } + super::select::warn_bounded(s, &bounded); if b(p, "select").unwrap_or(false) { s.map_mut().clear_selection(); for (id, rows) in found { diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index a5cdfe5e3..fa1adb1e7 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -6635,6 +6635,25 @@ fn analysis_workflows_models_scripts_and_raster_functions() { .run("analysis.neighborhood", &json!({"lon": 2.35, "lat": 48.85, "radius": 300000, "layers": ["Cities"], "select": true})) .unwrap(); assert!(n["layers"][0]["count"].as_u64().unwrap() >= 2); + // A center that's no place is refused, not taken as a pole (#1033). + assert!(s.run("analysis.neighborhood", &json!({"lon": 2.35, "lat": 1000, "radius": 300000})).is_err()); + // Lines near a pole with more vertices than one search measures: those past its steps are + // counted at the least distance they could be, and it says how many (#1033); a search that + // measures every one says nothing. + let mut ds = rustgis_data::Dataset::new("Spirals", rustgis_data::GeometryType::Polyline, vec![]); + for _ in 0..300 { + let line = (0..=348).map(|i| rustgis_data::Coord::new(if i % 2 == 0 { -179.5 } else { 179.5 }, 89.95)).collect(); + ds.push(Some(rustgis_data::Geometry::Line(vec![line])), vec![]); + } + s.map_mut().add_dataset(ds); + let mark = s.notice_mark(); + let n = s.run("analysis.neighborhood", &json!({"lon": 0, "lat": 89.9, "radius": 10000, "layers": ["Spirals"]})).unwrap(); + assert_eq!(n["layers"][0]["count"], 300); + let said: Vec<&str> = s.notices_since(mark).iter().map(|n| n.text.as_str()).collect(); + assert!(said.iter().any(|t| t.starts_with("299 feature(s) near a pole (Spirals 299)") && t.contains("beyond the radius")), "{said:?}"); + let mark = s.notice_mark(); + s.run("analysis.neighborhood", &json!({"lon": 0, "lat": 89.9, "radius": 10000, "layers": ["Cities"]})).unwrap(); + assert!(s.notices_since(mark).iter().all(|n| !n.text.contains("near a pole")), "{:?}", s.notices_since(mark)); s.run("select.byAttributes", &json!({"layer": "Countries", "where": "name = 'France'"})).unwrap(); let sim = s.run("analysis.simulate", &json!({"layer": "Countries", "count": 20, "seed": 3})).unwrap(); assert_eq!(sim["points"], 20); @@ -25517,6 +25536,58 @@ fn background_runs_follow_a_renamed_map_a_moved_window_and_a_loaded_layer() { let _ = std::fs::remove_file(path); } +#[test] +fn validate_topology_says_which_checks_are_incomplete() { + // A line 11 m from the North Pole winding round it some 19,000 times is measured coarsely, so + // Must Not Have Dangles is incomplete and the validation says so; a plain line, nothing (#1033). + use rustgis_data::{Coord, Dataset, Geometry, GeometryType}; + let c = Coord::new; + for (to, warned) in [(7_000_000.0, true), (1.0, false)] { + let mut s = Session::new(Project::default()); + let mut paths = Dataset::new("Paths", GeometryType::Polyline, vec![]); + paths.push(Some(Geometry::Line(vec![vec![c(0.0, 89.9999), c(to, 89.9999)]])), vec![]); + s.map_mut().add_dataset(paths); + s.run("topology.rules", &json!({"add": {"rule": "mustNotHaveDangles", "layer": "Paths"}})).unwrap(); + let mark = s.notice_mark(); + let v = s.run("topology.validate", &Value::Null).unwrap(); + assert_eq!(v["incomplete"].as_array().map(Vec::len), Some(usize::from(warned)), "{v}"); + let said = s + .notices_since(mark) + .iter() + .any(|n| n.text.starts_with("Checks that stopped early") && n.text.ends_with("Must Not Have Dangles")); + assert_eq!(said, warned, "{:?}", s.notices_since(mark)); + } +} + +#[test] +fn must_be_covered_by_says_which_polygons_it_left_out_at_a_pole() { + // Around the South Pole, a cap from 20°E round to 10°E and the 10° between written a + // billion degrees on: the check at the pole can't use that one, reports the point there and + // warns why; with a whole cap it says nothing of the kind (#1033). + use rustgis_data::{Coord, Dataset, Geometry, GeometryType}; + let c = Coord::new; + let ring = |w: f64, e: f64| Geometry::Polygon(vec![vec![vec![c(w, -85.0), c(e, -85.0), c(e, -90.0), c(w, -90.0), c(w, -85.0)]]]); + for (caps, warned) in [(vec![ring(20.0, 370.0), ring(1e9 + 90.0, 1e9 + 100.0)], true), (vec![ring(-180.0, 180.0)], false)] { + let mut s = Session::new(Project::default()); + let mut by = Dataset::new("Caps", GeometryType::Polygon, vec![]); + for g in caps { + by.push(Some(g), vec![]); + } + let mut pts = Dataset::new("Pole", GeometryType::Point, vec![]); + pts.push(Some(Geometry::Point(c(0.0, -90.0))), vec![]); + s.map_mut().add_dataset(by); + s.map_mut().add_dataset(pts); + s.run("topology.rules", &json!({"add": {"rule": "mustBeCoveredBy", "layer": "Pole", "by": "Caps"}})).unwrap(); + let mark = s.notice_mark(); + let v = s.run("topology.validate", &Value::Null).unwrap(); + assert_eq!(v["errors"], u64::from(warned)); + let said = s.notices_since(mark).iter().any(|n| { + n.text.contains("left out of the check at a pole") && n.text.contains("Must Be Covered By on 'Pole': 1 feature(s) of 'Caps'") + }); + assert_eq!(said, warned, "{:?}", s.notices_since(mark)); + } +} + #[test] fn a_pdf_says_which_frames_and_pages_it_draws_as_pictures() { let mut s = Session::new(Project::default()); diff --git a/crates/geom/src/lib.rs b/crates/geom/src/lib.rs index b6423d51a..4151844cb 100644 --- a/crates/geom/src/lib.rs +++ b/crates/geom/src/lib.rs @@ -1,7 +1,8 @@ //! L0: coordinates, extents, the Web Mercator projection, geodesic measurement and map scale. //! //! Everything here is plain math with no I/O. Longitude/latitude are WGS 84 degrees; projected -//! coordinates are Web Mercator metres (EPSG:3857). +//! coordinates are Web Mercator meters (EPSG:3857), which [`lonlat_to_plane`] continues beyond +//! 85.05° to the poles for analysis. #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] pub mod body; @@ -22,10 +23,17 @@ pub const EARTH_RADIUS_M: f64 = 6_371_008.8; /// 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; +/// The WGS 84 ellipsoid's radius of curvature at the poles, a²/b = a/(1 − f) (meters): near a pole, +/// distances along and across meridians are about this radius times the angle. +pub const WGS84_POLAR_RADIUS_M: f64 = body::EARTH.a / (1.0 - 1.0 / body::EARTH.inv_f); /// Web Mercator sphere radius (metres). pub const MERCATOR_R: f64 = 6_378_137.0; /// Latitude limit of Web Mercator. pub const MAX_LAT: f64 = 85.051_128_78; +/// How near ±90° (degrees) a latitude counts as the pole: about a centimeter, more than geometry +/// operations' rounding moves a pole vertex (around 1e-8°). For the polar edits and the features +/// meeting at a pole that build on this (#1033); here only the tests use it. +pub const POLE_TOLERANCE: f64 = 1e-7; /// Screen pixel size assumed when converting resolution to a map scale (0.28 mm, OGC). pub const PIXEL_M: f64 = 0.000_28; @@ -173,6 +181,49 @@ pub fn merc_to_lonlat(c: Coord) -> Coord { Coord::new(lon, lat) } +/// Web Mercator's northing at [`MAX_LAT`] north and south, each as [`lonlat_to_merc`] rounds it +/// (the two differ by a rounding), and how fast it grows there (meters a radian of latitude, +/// R·sec [`MAX_LAT`]): where [`lonlat_to_plane`] carries it on to the poles. +static MERC_EDGE: LazyLock<(f64, f64, f64)> = LazyLock::new(|| { + let edge = |lat: f64| lonlat_to_merc(Coord::new(0.0, lat)).y; + (edge(MAX_LAT), edge(-MAX_LAT), MERCATOR_R / MAX_LAT.to_radians().cos()) +}); + +/// Longitude/latitude (degrees) → the plane the analysis tools work in (meters): Web Mercator +/// within ±[`MAX_LAT`], exactly as [`lonlat_to_merc`] gives it, and beyond, where Web Mercator +/// stops, continued straight in latitude at the rate it reaches there, so the two meet smoothly and +/// the poles lie at ±26,424 km. Nothing within ±90° is clamped (a latitude past one is taken as +/// the pole), so a polar feature keeps its shape through an overlay. Beyond 85.05°, where the plane +/// is linear in latitude, an edge straight in the plane is straight in longitude and latitude too +/// (Web Mercator's view still draws it clamped to its edge). +pub fn lonlat_to_plane(c: Coord) -> Coord { + let lat = c.y.clamp(-90.0, 90.0); + // False for NaN, which then stays NaN below. + if lat.abs() <= MAX_LAT { + return lonlat_to_merc(c); + } + let (north, south, rate) = *MERC_EDGE; + let beyond = (lat.abs() - MAX_LAT).to_radians() * rate; + Coord::new(MERCATOR_R * c.x.to_radians(), if lat > 0.0 { north + beyond } else { south - beyond }) +} + +/// The analysis plane's meters → longitude/latitude (degrees), undoing [`lonlat_to_plane`]: +/// [`merc_to_lonlat`] within Web Mercator's square, and beyond it back along the continuation, to +/// ±90° at most. +pub fn plane_to_lonlat(p: Coord) -> Coord { + let (north, south, rate) = *MERC_EDGE; + let beyond = if p.y > north { + p.y - north + } else if p.y < south { + south - p.y + } else { + // Within the square, or NaN, which merc_to_lonlat keeps NaN. + return merc_to_lonlat(p); + }; + let lat = (MAX_LAT + (beyond / rate).to_degrees()).min(90.0); + Coord::new((p.x / MERCATOR_R).to_degrees(), lat.copysign(p.y)) +} + /// Projects a lon/lat extent to Web Mercator, preserving unwrapped longitudes for spans below /// 360°. A span of a full turn or more maps to the standard world bounds. pub fn extent_to_merc(e: &Extent) -> Extent { @@ -775,6 +826,43 @@ mod tests { assert!(spherical_centre(&[]).is_none()); } + #[test] + fn the_analysis_plane_is_web_mercator_within_its_square_and_reaches_the_poles() { + // Within ±MAX_LAT, to the bit what the tools used before, both ways, past ±180° too. + for i in -40..=40 { + let lat = MAX_LAT * f64::from(i) / 40.0; + for lon in [-540.0, -180.0, -0.5, 0.0, 33.3, 179.999, 180.0, 541.25] { + let c = Coord::new(lon, lat); + let (m, p) = (lonlat_to_merc(c), lonlat_to_plane(c)); + assert!(m.x.to_bits() == p.x.to_bits() && m.y.to_bits() == p.y.to_bits(), "{c:?}"); + let (a, b) = (merc_to_lonlat(m), plane_to_lonlat(m)); + assert!(a.x.to_bits() == b.x.to_bits() && a.y.to_bits() == b.y.to_bits(), "{c:?}"); + } + } + // Beyond, nothing is clamped: latitude climbs on to the poles, about 26,424 km out, meeting + // Web Mercator at its edge with the same slope (R·sec 85.05° a radian), and comes back. + let y = |lat: f64| lonlat_to_plane(Coord::new(10.0, lat)).y; + assert!((y(90.0) - 26_423_582.9).abs() < 1.0 && (y(-90.0) + y(90.0)).abs() < 1e-6, "{}", y(90.0)); + let h = 1e-6; + let (below, above) = ((y(MAX_LAT) - y(MAX_LAT - h)) / h, (y(MAX_LAT + h) - y(MAX_LAT)) / h); + assert!((above / below - 1.0).abs() < 1e-5, "{below} {above}"); + let mut last = y(MAX_LAT); + for i in 1..=100 { + let lat = MAX_LAT + (90.0 - MAX_LAT) * f64::from(i) / 100.0; + assert!(y(lat) > last, "{lat}"); + last = y(lat); + for sign in [1.0, -1.0] { + let back = plane_to_lonlat(lonlat_to_plane(Coord::new(-170.25, sign * lat))); + assert!((back.y - sign * lat).abs() < 1e-9 && (back.x + 170.25).abs() < 1e-12, "{lat}: {back:?}"); + } + } + // Junk stays junk or lands on a pole, never panics or goes past one. + assert!(lonlat_to_plane(Coord::new(0.0, f64::NAN)).y.is_nan() && plane_to_lonlat(Coord::new(0.0, f64::NAN)).y.is_nan()); + assert_eq!(lonlat_to_plane(Coord::new(0.0, f64::INFINITY)).y, y(90.0)); + assert_eq!(lonlat_to_plane(Coord::new(0.0, -1e300)).y, y(-90.0)); + assert_eq!((plane_to_lonlat(Coord::new(0.0, 1e300)).y, plane_to_lonlat(Coord::new(0.0, f64::NEG_INFINITY)).y), (90.0, -90.0)); + } + #[test] fn wgs84_direct_and_inverse_agree_with_geographiclib() { // GeographicLib's destinations (pyproj.Geod(ellps="WGS84").fwd), as in the test above. diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 00a28fd47..0c6b1c904 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -8,8 +8,9 @@ Lasso** (drag a freehand outline), or **Select by Radius** (drag from the centre selection; Ctrl/Cmd removes; Esc cancels an unfinished outline. Area tools use even-odd fill, so a self-crossing outline selects its alternating filled regions. Radius uses the Neighborhood Explorer's local Mercator distance and retained geographic prefilter; its outline and readout are approximate -and some features near the outline may be excluded. The main Select button remains click/rectangle -selection. +and some features near the outline may be excluded. A feature reaching beyond 85.05°, where Web +Mercator stops, is measured on the ground instead, so it may be selected outside the outline, which +stops at 85.05°, or left out inside it. The main Select button remains click/rectangle selection. All selection tools take only features that are drawn, in selectable layers: a definition query or a display filter hiding a feature keeps it out, and a click picks the feature drawn on top among the selectable layers. **Select By Attributes** takes a SQL clause and checks it as you type; **Select By @@ -105,6 +106,46 @@ written some other way (a token typed into a SQL query, say) can still be there: before sharing it. A point buffered where the map was clicked is copied as the command that buffered it, with the point. +**Near the poles.** The geometry tools and the topology checks work on Web Mercator's coordinates, +carried on beyond 85.05°, where Web Mercator stops, to the poles instead of clamped there, so polar +features keep their shape. An edge wholly beyond 85.05° is straight in longitude and latitude. One +that crosses 85.05° is straight in the tools, but the Web Mercator map clamps its far end onto +85.05° and draws it straight to there, so a cut along such an edge can fall well east or west of +where that map shows it: about 48 km at 84.5°N for an edge from (0°, 80°N) to (10°E, 89°N). A +pole is a line in these coordinates, not a point, so Minimum Bounding Geometry of points spread +around a pole spans their longitudes (a sliver, for points on one parallel), not the cap they +enclose on the ground. The globe draws a long edge along the great circle between its ends, so near +a pole it shows a little poleward of where the tools take it: about 850 m at the middle of a 10° +edge along 88°N (as everywhere, the tools take an edge as straight on the map, not on the globe). +Beyond 85.05°, the dangle check measures line ends about their pole (on WGS 84's curvature there, +within 0.13% out to 85°), right up to the pole and across ±180°; an edge winding round a pole so +many times that it can't be measured closely is measured coarsely, and the check says it's +incomplete. Features that both reach a pole +(to within 1e-7°, about a centimeter) meet there, one place on the ground, in Select By Location +and Spatial Join (Intersects), whatever their longitudes; in Clip, the overlays, Erase, Split, +Summarize Within and Spatial Join ("Have their center in") a point at a pole falls in a polygon +reaching it (a line or polygon meeting another only at a pole has no part in it, as at a corner +anywhere); and a point at a pole is covered (Must Be Covered By) by polygons that together surround +the pole (elsewhere, a point on the boundary between two polygons is still reported, as before), +not by ones that only touch it at a corner or leave a gap reaching it (a hole no other polygon +fills included). The check looks at circles around the pole from 1.1 cm to 8.7 km from it on the +ground, out to the first the polygons don't reach at all: a gap more than 9e-6° of longitude wide +where it crosses one counts, and narrower ones are taken for the seams the overlays' rounding +leaves. It leans to reporting: such a gap counts even if it doesn't reach the pole, and so do +polygons reaching less far from the pole on one side than on another, within 8.7 km. +Spatial Join's Within, Contains and Largest overlap still see a pole as a line. Radius selections +and the Neighborhood Explorer measure a feature reaching beyond 85.05°, or any feature from a center +beyond it, on the ground, following each edge as the tools take it (straight on the map) to within +about 125 m out to 5,000 km (farther, more, but only ever nearer); the +Neighborhood Explorer also takes a center beyond 85.05°, while radius selections take centers within +±85.05°. A feature within ±85.05° searched from a center within it is still measured in local +Mercator, which doesn't go over the pole, so from a high latitude a feature just within 85.05° on +the far side of the pole can be left out while one just beyond it, nearer on the ground, is found. +A search measures at most 250,000 points and steps along edges in each layer, in steps of at most +half a degree; past them, or along an edge over 512° long (as only a hostile file has), a feature is +counted at the least distance it could be, so it may lie beyond the radius, and the tool warns how +many were, layer by layer. + These are the built-in tools; GeoLibre and plugin tools join them when available. Use these ids in `gp.run`'s `tool` parameter (`gp.tools` lists the available tools and their parameters): @@ -209,8 +250,11 @@ gp.run={"tool":"vectorToPmtiles","params":{"input":"Roads","minZoom":0,"maxZoom" least three different, non-collinear corners (plus the closing vertex), so a tolerance larger than a polygon produces a small ring instead of restoring its original vertices. Holes and multipart rings are retained. This is not a topology-preserving simplifier: check the result -where holes or neighbouring polygons must stay aligned. The default tolerance remains 25 km; -use a smaller value for local boundaries. +where holes or neighboring polygons must stay aligned. The default tolerance remains 25 km; +use a smaller value for local boundaries. A line or polygon has its tolerance taken east–west at +the scale of its middle (its label point's latitude), at most 85°'s, so beyond 85.05°, where the +ground is far narrower, less is simplified away than asked (a 3 km tolerance at 89.9°N takes out +only what lies within about 60 m east–west), and the tool says how many reach that far. **Kernel Density** (`kernelDensity`) calculates a continuous magnitude-per-unit-area raster surface from point features using a Quartic (biweight, default) or Gaussian kernel (truncated at 3σ and @@ -224,7 +268,8 @@ e.g. `EPSG:32618`). Because bandwidth and cell size are measured in coordinate-space linear units, global cylindrical projections such as Web Mercator (EPSG:3857) introduce scale distortion that increases with latitude; use a local projected coordinate system (such as UTM or State Plane) when accurate ground distances -are required. Leaving `radius` empty or setting it to `0` automatically estimates the bandwidth +are required. A point the system can't hold, such as one beyond 85.05° in Web Mercator, is left +out, and the tool says how many were. Leaving `radius` empty or setting it to `0` automatically estimates the bandwidth using Silverman's rule of thumb based on the points' spatial distribution (falling back to extent span / 30, or 1000 coordinate-space linear units for a single point or coincident points). When `pixelSize` is omitted, an adaptive cell size is calculated to span approximately 250 cells along the longer side @@ -415,7 +460,8 @@ along it are kept, but polygons that only share an edge or a corner, and a line one, give nothing. The input can be points, lines or polygons and keeps its geometry type; the overlay must be polygons. An overlay field whose name the input already has gets a suffix (`id_2`). Overlapping features within one layer are not split against each other. Like Clip and Erase, the -overlay works on the map's Web Mercator coordinates, so the parts' edges are straight on the map. +overlay works on the map's Web Mercator coordinates, so the parts' edges are straight on the map +within ±85.05° (beyond it, see **Near the poles** above). Selections in both layers are honoured. To measure each part, run Calculate Geometry Attributes on the result: @@ -726,15 +772,19 @@ billion), so a part that leaves nothing outside its zone wider than that grid is and only parts outside every zone wider than the grid are reported. A feature so narrow next to its zones that the rounding can move its area by more than 1% is measured only roughly, unless it lies wholly inside a zone, and the tool says how many there are. The percentages of classes that overlap -each other can add up to more than 100. Like the other overlay tools, it works on the map, which -ends at 85.05° north and south: parts of polygons beyond that are taken as lying along it, so -overlaps there have no area ([#1033](https://github.com/opengeos/RustGIS/issues/1033)); the tool -says how many features reach beyond it. +each other can add up to more than 100. Like the other overlay tools, it works on the map's +coordinates carried on past 85.05° to the poles (see **Near the poles** above), so polygons there +keep their area. Near a pole an edge along a parallel runs far from the shortest way between its +ends, which `area_km2` follows, so the two measures can differ widely: for a quarter cap from 88°N +to the pole, written with one edge along 88°N, `area_km2` is about 64% of the area its +percentage is taken from. Add vertices along such edges to bring them together. The tool says how many zone and class features have no polygon (no geometry, an empty one, a line or point, or a ring collapsed to no area; they are left out), how many cross the antimeridian without being split there (the map, and with it the percentages and shares, takes them the long way round, -while area_km2 takes their ground shape), how many zones overlap no class feature and so have no +while area_km2 takes their ground shape), how many have an edge beyond 85.05° running east or west +(taken as the map draws it, so their percentages and area_km2 can differ widely, as above), how +many zones overlap no class feature and so have no rows, how many class features have no value in a sum field (they add nothing to its shares) and, with sum fields, how many class features lie partly or wholly outside the zones. Only selected features are used when a layer has a selection, and the table is named `_Tabulate` unless diff --git a/docs/user-guide/editing.md b/docs/user-guide/editing.md index 8977e77a6..856146aff 100644 --- a/docs/user-guide/editing.md +++ b/docs/user-guide/editing.md @@ -95,6 +95,8 @@ you edit the value or choose units; click **Apply elevation** to apply the new i (Move, and `edit.traverse`, which draws a line or polygon from distances and bearings). - **Topology**: rules (must not overlap, must not have gaps, must not self-intersect, must not have dangles, must be covered by); Validate lists the errors, and clicking one zooms to it. + Beyond 85.05° they follow the geometry tools (see + [Geoprocessing tools](analysis.md#geoprocessing-tools), **Near the poles**). - **Status**: features added, changed and deleted per layer since the last save. ### Moving shared polygon vertices