fix(analysis): keep buffer outlines at their distance between vertices

crates/analysis/src/geodesic.rs

@@ -5,12 +5,20 @@
 //! from the centre are true and the error elsewhere grows with the square of the distance from it.
 //! So a small feature is buffered in one piece, and a long line or a large polygon's outline in
 //! stretches of a few hundred kilometres, each in its own projection, and the pieces are joined.
+//!
+//! 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
+//! within 0.01% of the distance (or 1 cm). Arcs are drawn as `geo` draws them by default, close to
+//! GEOS and PostGIS's 8 edges to a quarter circle: a point's circle with 32 edges, a line's caps and
+//! a polygon's corners with 15 to a half circle. Between an arc's vertices the outline runs inside
+//! it, by up to 0.48% of the distance round a point and 0.55% at a cap or corner.
 
 use geo::{
     BooleanOps, Buffer, Contains, Coord as GC, LineString, MultiLineString, MultiPolygon, Orient, Point, Polygon, orient::Direction,
 };
 use rustgis_data::{Coord, Geometry};
-use rustgis_geom::{Aeqd, geodesic_interpolate, haversine_m, spherical_centre};
+use rustgis_geom::{Aeqd, geodesic_interpolate, haversine_m, lonlat_to_merc, merc_to_lonlat, spherical_centre};
 
 /// A feature within this distance of its centre is buffered in one projection (the error is
 /// under 0.01%).
@@ -22,6 +30,20 @@ const STRETCH_M: f64 = 300_000.0;
 const DENSIFY_M: f64 = 10_000.0;
 /// At most this many vertices are added to one line or ring (four times round the Earth at 10 km).
 const MAX_ADDED: f64 = 16_000.0;
+/// Back in lon/lat, an outline's straight edges stay within this fraction of the distance of the
+/// edges drawn in the plane, or [`MIN_TOLERANCE_M`], whichever is more.
+const OUTLINE_TOLERANCE: f64 = 2e-5;
+/// The smallest tolerance (metres): finer only adds vertices.
+const MIN_TOLERANCE_M: f64 = 0.01;
+/// An edge is halved at most this many times on its way back to lon/lat (4096 pieces): enough for
+/// the longest edges a buffer has (about 300 km) to come within 1 cm where they stray most.
+const MAX_SPLITS: u32 = 12;
+/// How far, per square metre of an edge's length, a straight edge in lon/lat can stray from the
+/// same edge straight in the plane: measured at 2e-8 in mid-latitudes (the long side of a 2 km
+/// buffer of a 96 km line strayed 190 m, 2e-8 × 96,000²) and 1.9e-7 at worst, along a 298 km
+/// line at 84°N, with a margin. An edge too short to stray as far as the tolerance isn't
+/// checked, so an arc's short edges cost nothing.
+const STRAY_PER_M2: f64 = 3e-7;
 
 /// A feature's buffer at `metres` (negative shrinks a polygon), in lon/lat; `None` when nothing
 /// is left.
@@ -53,15 +75,23 @@ impl Plane {
         LineString::new(pts.iter().map(|c| self.0.forward(*c)).map(|p| GC { x: p.x, y: p.y }).collect())
     }
 
-    /// A buffer drawn in the plane, in lon/lat degrees.
-    fn back(&self, mp: &MultiPolygon<f64>) -> MultiPolygon<f64> {
-        MultiPolygon::new(
-            mp.0.iter().map(|p| Polygon::new(self.ring(p.exterior()), p.interiors().iter().map(|r| self.ring(r)).collect())).collect(),
-        )
+    /// A buffer at `metres` drawn in the plane, in lon/lat degrees.
+    fn back(&self, mp: &MultiPolygon<f64>, metres: f64) -> MultiPolygon<f64> {
+        let tolerance = (metres.abs() * OUTLINE_TOLERANCE).max(MIN_TOLERANCE_M);
+        let ring = |r: &LineString<f64>| self.ring(r, tolerance);
+        MultiPolygon::new(mp.0.iter().map(|p| Polygon::new(ring(p.exterior()), p.interiors().iter().map(ring).collect())).collect())
     }
 
-    /// A ring drawn in the plane, in lon/lat degrees.
-    fn ring(&self, r: &LineString<f64>) -> LineString<f64> {
+    /// A ring drawn in the plane, in lon/lat degrees, its edges split as [`Plane::split`] needs.
+    fn ring(&self, r: &LineString<f64>, tolerance: f64) -> LineString<f64> {
+        let mut dense = Vec::with_capacity(r.0.len());
+        for w in r.0.windows(2) {
+            let [a, b] = w else { continue };
+            dense.push(*a);
+            self.split(*a, *b, tolerance, 0, &mut dense);
+        }
+        dense.extend(r.0.last());
+        let r = &LineString::new(dense);
         let lonlat: Vec<Coord> = r.0.iter().map(|p| self.0.inverse(Coord::new(p.x, p.y))).collect();
         let around = |lat: f64| {
             let pole = self.0.forward(Coord::new(0.0, lat));
@@ -82,6 +112,44 @@ impl Plane {
         LineString::new(out)
     }
 
+    /// 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
+    /// line crossing the equator); an edge that strays is halved, and each half checked again, up
+    /// to [`MAX_SPLITS`] times.
+    fn split(&self, a: GC<f64>, b: GC<f64>, tolerance: f64, depth: u32, out: &mut Vec<GC<f64>>) {
+        if depth >= MAX_SPLITS || STRAY_PER_M2 * ((b.x - a.x).powi(2) + (b.y - a.y).powi(2)) <= tolerance {
+            return;
+        }
+        let at = |t: f64| self.0.inverse(Coord::new(a.x + t * (b.x - a.x), a.y + t * (b.y - a.y)));
+        let (la, lb) = (at(0.0), at(1.0));
+        let beyond = |c: Coord| c.y.abs() >= rustgis_geom::MAX_LAT;
+        let strays = 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.
+            (true, false) | (false, true) => true,
+            (false, false) => {
+                // `b` on `a`'s side of the antimeridian, so the straight edge is the short one.
+                let lb = Coord::new(la.x + wrap(lb.x - la.x), lb.y);
+                let (ma, mb) = (lonlat_to_merc(la), lonlat_to_merc(lb));
+                [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
+                })
+            }
+        };
+        if !strays {
+            return;
+        }
+        let mid = GC { x: (a.x + b.x) / 2.0, y: (a.y + b.y) / 2.0 };
+        self.split(a, mid, tolerance, depth + 1, out);
+        out.push(mid);
+        self.split(mid, b, tolerance, depth + 1, out);
+    }
+
     /// A ring round a pole, which goes once round the globe in longitude, as a lon/lat ring: cut at
     /// the antimeridian and closed along the pole, so in lon/lat it holds the cap between the ring and
     /// the top (or bottom) of the map.
@@ -123,24 +191,28 @@ fn wrap(d: f64) -> f64 {
     (d + 180.0).rem_euclid(360.0) - 180.0
 }
 
+/// A point's buffer: a circle drawn in a projection centred on the point.
 fn point(c: Coord, metres: f64) -> MultiPolygon<f64> {
     if !c.is_finite() {
         return MultiPolygon::new(vec![]);
     }
     let plane = Plane::new(c, c.x);
-    plane.back(&Point::new(0.0, 0.0).buffer(metres))
+    plane.back(&Point::new(0.0, 0.0).buffer(metres), metres)
 }
 
+/// A line's buffer: in one projection centred on it when it fits, else in stretches.
 fn lines(parts: &[Vec<Coord>], metres: f64) -> MultiPolygon<f64> {
     let all: Vec<Coord> = parts.iter().flatten().copied().collect();
     let Some(centre) = spherical_centre(&all) else { return MultiPolygon::new(vec![]) };
     if fits(&all, centre) {
         let plane = Plane::new(centre, all.first().map_or(centre.x, |c| c.x));
-        return plane.back(&MultiLineString::new(parts.iter().map(|p| plane.line(p)).collect()).buffer(metres));
+        return plane.back(&MultiLineString::new(parts.iter().map(|p| plane.line(p)).collect()).buffer(metres), metres);
     }
     union(stretched(parts, metres))
 }
 
+/// A polygon's buffer: in one projection centred on it when it fits, else the polygon with its
+/// outline's buffer (in stretches) added or, inward, taken away.
 fn polygons(polys: &[Vec<Vec<Coord>>], metres: f64) -> MultiPolygon<f64> {
     let all: Vec<Coord> = polys.iter().flatten().flatten().copied().collect();
     let Some(centre) = spherical_centre(&all) else { return MultiPolygon::new(vec![]) };
@@ -158,7 +230,7 @@ fn polygons(polys: &[Vec<Vec<Coord>>], metres: f64) -> MultiPolygon<f64> {
     };
     if fits(&all, centre) {
         let plane = Plane::new(centre, all.first().map_or(centre.x, |c| c.x));
-        return plane.back(&shape(&|r| plane.line(r)).buffer(metres));
+        return plane.back(&shape(&|r| plane.line(r)).buffer(metres), metres);
     }
     // A buffer distributes over a union: the polygon's buffer is the polygon with its outline's
     // buffer added (or, inward, taken away). The polygon stays in lon/lat, along the same densified
@@ -193,7 +265,7 @@ fn stretched(parts: &[Vec<Coord>], metres: f64) -> Vec<MultiPolygon<f64>> {
             let Some(stretch) = d.get(start..=end.min(d.len().saturating_sub(1))) else { break };
             if let (Some(centre), Some(first)) = (spherical_centre(stretch), stretch.first()) {
                 let plane = Plane::new(centre, first.x);
-                out.push(plane.back(&plane.line(stretch).buffer(metres)));
+                out.push(plane.back(&plane.line(stretch).buffer(metres), metres));
             }
             // Consecutive stretches share a vertex, so their buffers overlap.
             if end + 1 >= d.len() {
@@ -205,20 +277,36 @@ fn stretched(parts: &[Vec<Coord>], metres: f64) -> Vec<MultiPolygon<f64>> {
     out
 }
 
-/// A line with every edge split along its great circle into pieces of at most [`DENSIFY_M`], or
+/// A line with every edge split along its geodesic on WGS 84 into pieces of at most [`DENSIFY_M`], or
 /// longer ones for a line so long it would gain more than [`MAX_ADDED`] vertices.
 fn densify(pts: &[Coord]) -> Vec<Coord> {
     let step = DENSIFY_M.max(rustgis_geom::geodesic_length_m(pts) / MAX_ADDED);
     let mut out: Vec<Coord> = pts.first().copied().into_iter().collect();
     for w in pts.windows(2) {
         let [a, b] = w else { continue };
         let n = (haversine_m(*a, *b) / step).ceil().clamp(1.0, MAX_ADDED) as usize;
-        out.extend((1..n).map(|i| geodesic_interpolate(*a, *b, i as f64 / n as f64)));
+        out.extend((1..n).map(|i| on_wgs84(*a, *b, i as f64 / n as f64)));
         out.push(*b);
     }
     out
 }
 
+/// The point a fraction `t` of the way from `a` to `b` along the geodesic between them on WGS 84,
+/// the path the rest of the buffer follows (a great circle runs tens of metres from it over a
+/// few hundred kilometres). Coordinates that aren't on the Earth (a latitude beyond ±90°, say)
+/// have no geodesic, so they get the great circle, which stays finite.
+fn on_wgs84(a: Coord, b: Coord, t: f64) -> Coord {
+    let on_earth = |c: Coord| c.x.is_finite() && (-90.0..=90.0).contains(&c.y);
+    if on_earth(a) && on_earth(b) {
+        let p = geo::InterpolatePoint::point_at_ratio_between(&geo::Geodesic, geo::Point::new(a.x, a.y), geo::Point::new(b.x, b.y), t);
+        let p = Coord::new(a.x + wrap(p.x() - a.x), p.y());
+        if p.is_finite() {
+            return p;
+        }
+    }
+    geodesic_interpolate(a, b, t)
+}
+
 fn union(pieces: Vec<MultiPolygon<f64>>) -> MultiPolygon<f64> {
     if pieces.len() == 1 {
         return pieces.into_iter().next().unwrap_or_else(|| MultiPolygon::new(vec![]));

crates/analysis/src/tools.rs

@@ -4853,6 +4853,115 @@ mod tests {
         }
     }
 
+    /// Points along each edge of each ring at the fractions `ts`, as the tools draw edges: straight
+    /// in Web Mercator.
+    fn along_edges(polys: &[Vec<Vec<Coord>>], ts: &[f64]) -> Vec<Coord> {
+        use rustgis_geom::{lonlat_to_merc, merc_to_lonlat};
+        polys
+            .iter()
+            .flatten()
+            .flat_map(|ring| ring.windows(2).map(|w| (lonlat_to_merc(w[0]), lonlat_to_merc(w[1]))).collect::<Vec<_>>())
+            .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)))))
+            .collect()
+    }
+
+    /// Points along the geodesic from `a` to `b` on WGS 84, at most 20 m apart.
+    fn geodesic_samples(a: Coord, b: Coord) -> Vec<Coord> {
+        use geo::InterpolatePoint;
+        let n = (rustgis_geom::haversine_m(a, b) / 20.0).ceil().max(1.0) as usize;
+        (0..=n)
+            .map(|k| geo::Geodesic.point_at_ratio_between(geo::Point::new(a.x, a.y), geo::Point::new(b.x, b.y), k as f64 / n as f64))
+            .map(|p| Coord::new(p.x(), p.y()))
+            .collect()
+    }
+
+    /// The distance (metres) from `c` to a geodesic given by its [`geodesic_samples`], and whether
+    /// its nearest point is between the ends: the nearest sample on the sphere (among every 50th,
+    /// then its neighbours), then the nearest point on the ellipsoid along the geodesic either side
+    /// of it, by golden-section search to well under a millimetre.
+    fn m_to_geodesic(c: Coord, samples: &[Coord]) -> (f64, bool) {
+        use geo::InterpolatePoint;
+        let near = |range: std::ops::Range<usize>, step: usize| {
+            range
+                .step_by(step)
+                .min_by(|&i, &j| rustgis_geom::haversine_m(c, samples[i]).total_cmp(&rustgis_geom::haversine_m(c, samples[j])))
+        };
+        let coarse = near(0..samples.len(), 50).unwrap_or(0);
+        let i = near(coarse.saturating_sub(50)..(coarse + 51).min(samples.len()), 1).unwrap_or(coarse);
+        let (lo, hi) = (samples[i.saturating_sub(1)], samples[(i + 1).min(samples.len() - 1)]);
+        let at = |t: f64| {
+            let p = geo::Geodesic.point_at_ratio_between(geo::Point::new(lo.x, lo.y), geo::Point::new(hi.x, hi.y), t);
+            on_wgs84_m(c, Coord::new(p.x(), p.y()))
+        };
+        let (mut a, mut b) = (0.0, 1.0);
+        let g = (5f64.sqrt() - 1.0) / 2.0;
+        for _ in 0..60 {
+            let (x1, x2) = (b - g * (b - a), a + g * (b - a));
+            if at(x1) <= at(x2) { b = x2 } else { a = x1 }
+        }
+        (at((a + b) / 2.0).min(at(0.0)).min(at(1.0)), i > 0 && i + 1 < samples.len())
+    }
+
+    /// 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
+    /// within 0.01% of the distance (or 1 cm) along a line's sides, and no further inside an arc
+    /// than its edges cut: 32 to a point's circle and 15 to a half circle at a line's cap, as `geo`
+    /// draws them.
+    #[test]
+    fn buffer_outlines_keep_their_distance_between_vertices() {
+        let fractions = [0.25, 0.5, 0.75];
+        let points = |polys: &[Vec<Vec<Coord>>]| -> Vec<Coord> {
+            polys.iter().flatten().flatten().copied().chain(along_edges(polys, &fractions)).collect()
+        };
+        // An arc of n edges to a circle cuts d(1 − cos(π/n)) inside it between vertices.
+        let sag = |d: f64, n: f64| d * (1.0 - (std::f64::consts::PI / n).cos());
+
+        let centre = Coord::new(10.125, 46.125);
+        for c in points(&buffered(Geometry::Point(centre), 1.0)) {
+            let d = on_wgs84_m(centre, c);
+            assert!(d < 1_000.1 && d > 1_000.0 - sag(1_000.0, 32.0) - 0.1, "{c:?} is {d} m from the centre");
+        }
+
+        // A long east–west line; a diagonal one across the equator, where the stray changes side
+        // halfway along an edge (a 10 m buffer of it strayed 0.65 m at the quarter points); and a
+        // 298 km one at 84°N, whose sides needed more halvings than the limit first allowed.
+        for (a, b, km) in [
+            (Coord::new(9.875, 46.25), Coord::new(11.125, 46.25), 2.0),
+            (Coord::new(-0.6, -0.6), Coord::new(0.6, 0.6), 0.01),
+            (Coord::new(-12.559_750_946_199_074, 83.853_999_752_232_98), Coord::new(12.559_750_946_199_074, 83.853_999_752_232_98), 0.01),
+        ] {
+            let (d0, samples) = (km * 1000.0, geodesic_samples(a, b));
+            for c in points(&buffered(Geometry::Line(vec![vec![a, b]]), km)) {
+                let (d, beside) = m_to_geodesic(c, &samples);
+                // Along the sides, within 0.01% or 1 cm (and the sampler's millimetre); round the
+                // caps, within the sag of their 15 edges to a half circle.
+                let side = (d0 * 1e-4).max(0.01) + 1e-3;
+                let bound = if beside { side } else { sag(d0, 30.0) + side };
+                assert!(d < d0 + side && d > d0 - bound, "{km} km: {c:?} is {d} m from the line");
+            }
+        }
+
+        // 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.
+        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 d = on_wgs84_m(north, c);
+            assert!(d < 1_000_000.0 * (1.0 + 1e-4) && d > 1_000_000.0 - sag(1_000_000.0, 32.0) - 100.0, "{c:?} is {d} m from {north:?}");
+        }
+
+        // Coordinates off the Earth have no geodesic: buffering them doesn't panic.
+        let off = [Coord::new(0.0, 100.0), Coord::new(20.0, 100.0)];
+        let mut ds = Dataset::new("f", GeometryType::Polyline, vec![]);
+        ds.push(Some(Geometry::Line(vec![off.to_vec()])), vec![]);
+        let m = Map(HashMap::from([("f".to_string(), ds)]));
+        if let Err(e) = crate::run("buffer", &m, &serde_json::json!({"input": "f", "distance": 1})) {
+            assert!(!e.contains("unexpectedly"), "{e}");
+        }
+    }
+
     /// Buffers are geodesic (#394): Web Mercator's scale changes with latitude, so one scale for
     /// the whole feature made a 1000 km circle at 60°N reach 872 km north and 1142 km south, and a
     /// 50 km buffer of a line from 5°S to 33°S 59 km wide at one end and 50 km at the other.

docs/user-guide/analysis.md

@@ -211,8 +211,12 @@ is the layer name; `vectorToPmtiles` uses it as a file destination.
 **Buffer** distances, and the distance of **Select Layer By Location**, are measured on the Earth
 (the WGS 84 ellipsoid, as PROJ and GeographicLib measure it), not in the map's display
 projection. Each feature is buffered in an azimuthal equidistant projection on the ellipsoid
-centred on it, so a point's buffer is exact; long lines and large polygon outlines are
-buffered in stretches of a few hundred kilometres to limit projection error. Negative Buffer
+centred on it, so distances from a point are true; long lines and large polygon outlines are
+buffered in stretches of a few hundred kilometres to limit projection error. Up to 85° of
+latitude, a buffer's sides are within 0.01 % of the distance (or 1 cm). Its curves are drawn with
+straight edges (32 to a point's circle, 15 to a half circle at a line's end or a polygon's corner,
+close to GEOS and PostGIS's 8 to a quarter circle), whose corners are at the distance and whose
+middles run up to about 0.5 % of it inside the true curve. Negative Buffer
 distances shrink polygons. Buffer and Select Layer By Location limit distances to 5,000 km.
 **Calculate Geometry Attributes** measures polygon areas and perimeters, and line lengths,
 on the WGS 84 ellipsoid, independently of the map's display projection. Holes are subtracted