From 5a15379571491403e97cb456a6632362045a2ce0 Mon Sep 17 00:00:00 2001 From: mavrick Date: Sat, 31 Oct 2026 23:38:24 -0500 Subject: [PATCH 1/2] fix(analysis): refine buffer arcs to meet reference tolerance --- crates/analysis/src/geodesic.rs | 30 ++++++++++++++------- crates/analysis/src/tools.rs | 16 +++++------ crates/analysis/tests/reference/cases.json | 12 --------- crates/analysis/tests/reference/generate.py | 10 +++---- docs/user-guide/analysis.md | 8 +++--- 5 files changed, 36 insertions(+), 40 deletions(-) diff --git a/crates/analysis/src/geodesic.rs b/crates/analysis/src/geodesic.rs index e84645240..9443f4e44 100644 --- a/crates/analysis/src/geodesic.rs +++ b/crates/analysis/src/geodesic.rs @@ -15,13 +15,15 @@ //! stretch's centre, where the sides are within 0.05%. A ring gains at most [`MAX_ADDED`] points, //! which only a ring of thousands of long edges uses up; past that, its edges are left whole. //! -//! 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. +//! Arcs use a 1.8° step: 200 edges to a point's circle and about 100 to a semicircular cap. +//! Between vertices, those edges cut at most 0.013% of the distance inside the circle in the +//! projection plane. The buffer backend rounds a corner's segment count down, so a short corner +//! can have a step approaching 3.6° and cut up to 0.05% inside, on top of the projection error. use geo::{ - BooleanOps, Buffer, Contains, Coord as GC, LineString, MultiLineString, MultiPolygon, Orient, Point, Polygon, orient::Direction, + BooleanOps, Buffer, Contains, Coord as GC, LineString, MultiLineString, MultiPolygon, Orient, Point, Polygon, + algorithm::buffer::{BufferStyle, LineCap, LineJoin}, + orient::Direction, }; use rustgis_data::{Coord, Geometry}; use rustgis_geom::{Aeqd, geodesic_interpolate, haversine_m, lonlat_to_merc, merc_to_lonlat, spherical_centre}; @@ -52,6 +54,15 @@ const MAX_SPLITS: u32 = 12; /// checked, so an arc's short edges cost nothing. const STRAY_PER_M2: f64 = 3e-7; +/// The finest round-cap/join step supported by the buffer backend (1.8°). A circle's area +/// deficit is under 0.017%, leaving room for projection error in the reference's 0.1% bound. +const ARC_STEP: f64 = std::f64::consts::PI / 100.0; + +/// Round arcs at the same density for points, line caps and polygon corners, including stretches. +fn buffer_style(metres: f64) -> BufferStyle { + BufferStyle::new(metres).line_cap(LineCap::Round(ARC_STEP)).line_join(LineJoin::Round(ARC_STEP)) +} + /// A feature's buffer at `metres` (negative shrinks a polygon), in lon/lat; `None` when nothing /// is left. pub fn geodesic_buffer(g: &Geometry, metres: f64) -> Option { @@ -210,7 +221,7 @@ fn point(c: Coord, metres: f64) -> MultiPolygon { return MultiPolygon::new(vec![]); } let plane = Plane::new(c, c.x); - plane.back(&Point::new(0.0, 0.0).buffer(metres), metres) + plane.back(&Point::new(0.0, 0.0).buffer_with_style(buffer_style(metres)), metres) } /// A line's buffer: in one projection centred on it when it fits, else in stretches. @@ -219,7 +230,8 @@ fn lines(parts: &[Vec], metres: f64) -> MultiPolygon { 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), metres); + return plane + .back(&MultiLineString::new(parts.iter().map(|p| plane.line(p)).collect()).buffer_with_style(buffer_style(metres)), metres); } union(stretched(parts, metres)) } @@ -243,7 +255,7 @@ fn polygons(polys: &[Vec>], metres: f64) -> MultiPolygon { }; 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), metres); + return plane.back(&shape(&|r| plane.line(r)).buffer_with_style(buffer_style(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 @@ -278,7 +290,7 @@ fn stretched(parts: &[Vec], metres: f64) -> Vec> { 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), metres)); + out.push(plane.back(&plane.line(stretch).buffer_with_style(buffer_style(metres)), metres)); } // Consecutive stretches share a vertex, so their buffers overlap. if end + 1 >= d.len() { diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index a18d41843..61402dd45 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -4991,21 +4991,21 @@ mod tests { /// 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.012% 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. + /// than 0.013% 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. #[test] fn buffer_outlines_keep_their_distance_between_vertices() { let fractions = [0.25, 0.5, 0.75]; let points = |polys: &[Vec>]| -> Vec { 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()); + // Bound the geometric error independently of how many vertices the buffer emits. + let arc_error = |d: f64| d * 1.3e-4; 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"); + assert!(d < 1_000.1 && d > 1_000.0 - arc_error(1_000.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 @@ -5034,9 +5034,9 @@ mod tests { 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.012% or 1 cm (and the sampler's millimetre); round the - // caps, within the sag of their 15 edges to a half circle. + // caps, add the 0.013% arc approximation in the projection plane. let side = (d0 * 1.2e-4).max(0.01) + 1e-3; - let bound = if beside { side } else { sag(d0, 30.0) + side }; + let bound = if beside { side } else { arc_error(d0) + side }; assert!(d < d0 + side && d > d0 - bound, "{km} km: {c:?} is {d} m from the line"); } } @@ -5047,7 +5047,7 @@ mod tests { 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 + 1.2e-4) && d > 1_000_000.0 - sag(1_000_000.0, 32.0) - 100.0, "{c:?} is {d} m from {north:?}"); + 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:?}"); } // Coordinates off the Earth have no geodesic: buffering them doesn't panic. diff --git a/crates/analysis/tests/reference/cases.json b/crates/analysis/tests/reference/cases.json index c725501bf..59f8f5f8c 100644 --- a/crates/analysis/tests/reference/cases.json +++ b/crates/analysis/tests/reference/cases.json @@ -912,10 +912,6 @@ "expect": { "features": "expected/buffer_points_1km.geojson", "tolerance": 0.001 - }, - "known": { - "issue": "#1159: Buffer outlines sit up to 10 % of the distance off", - "tolerance": 0.1 } }, { @@ -929,10 +925,6 @@ "expect": { "features": "expected/buffer_stations_120km.geojson", "tolerance": 0.001 - }, - "known": { - "issue": "#1159: Buffer outlines sit up to 10 % of the distance off", - "tolerance": 0.1 } }, { @@ -946,10 +938,6 @@ "expect": { "features": "expected/buffer_roads_2km.geojson", "tolerance": 0.001 - }, - "known": { - "issue": "#1159: Buffer outlines sit up to 10 % of the distance off", - "tolerance": 0.1 } }, { diff --git a/crates/analysis/tests/reference/generate.py b/crates/analysis/tests/reference/generate.py index ffa23df51..5417d16b9 100644 --- a/crates/analysis/tests/reference/generate.py +++ b/crates/analysis/tests/reference/generate.py @@ -504,16 +504,12 @@ def buffered(name, rows, d): return write_expected(name, out) tol = 1e-3 - # Arcs are drawn with 32 edges, and long edges come back from the plane as straight lines in - # lon/lat: a 2 km buffer of a 96 km road is up to 190 m off. - # Until then each must still agree within 10 %, so a buffer that breaks further fails the test. - known = {"issue": "#1159: Buffer outlines sit up to 10 % of the distance off", "tolerance": 0.1} case("buffer points", "buffer", {"input": "points", "distance": 1, "units": "Kilometers"}, - {"features": buffered("buffer_points_1km", POINTS, 1000.0), "tolerance": tol}, known) + {"features": buffered("buffer_points_1km", POINTS, 1000.0), "tolerance": tol}) case("buffer points far", "buffer", {"input": "stations", "distance": 120, "units": "Kilometers"}, - {"features": buffered("buffer_stations_120km", STATIONS, 120_000.0), "tolerance": tol}, known) + {"features": buffered("buffer_stations_120km", STATIONS, 120_000.0), "tolerance": tol}) case("buffer lines", "buffer", {"input": "roads", "distance": 2000, "units": "Meters"}, - {"features": buffered("buffer_roads_2km", ROADS, 2000.0), "tolerance": tol}, known) + {"features": buffered("buffer_roads_2km", ROADS, 2000.0), "tolerance": tol}) case("buffer a long line", "buffer", {"input": "long", "distance": 5, "units": "Kilometers"}, {"features": buffered("buffer_long_5km", LONG, 5000.0), "tolerance": tol}) case("buffer polygons", "buffer", {"input": "parcels", "distance": 1, "units": "Kilometers"}, diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 914d5621b..8dd1067ac 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -214,10 +214,10 @@ projection. Each feature is buffered in an azimuthal equidistant projection on t 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.012 % of the distance (or 1 cm), and within 0.05 % along -a winding line or outline buffered in stretches. 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 +a winding line or outline buffered in stretches. Its curves are drawn with straight edges (200 +to a point's circle, about 100 to a semicircular line end). Between vertices, those edges cut up +to 0.013 % of the distance inside the true curve in the projection plane, and up to 0.05 % at +polygon corners and line joins, in addition to the projection error. 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 From 7a5c7b7a363727fce6a532418e058220762941b5 Mon Sep 17 00:00:00 2001 From: mavrick Date: Sun, 1 Nov 2026 01:09:17 -0600 Subject: [PATCH 2/2] fix(analysis): bound buffer output and reduce arc overhead --- crates/analysis/src/geodesic.rs | 73 ++++++++++++++++++++++++++++----- crates/analysis/src/tools.rs | 51 ++++++++++++++++++++--- docs/user-guide/analysis.md | 10 +++-- 3 files changed, 115 insertions(+), 19 deletions(-) diff --git a/crates/analysis/src/geodesic.rs b/crates/analysis/src/geodesic.rs index 9443f4e44..d4bd2e0a3 100644 --- a/crates/analysis/src/geodesic.rs +++ b/crates/analysis/src/geodesic.rs @@ -15,13 +15,13 @@ //! stretch's centre, where the sides are within 0.05%. A ring gains at most [`MAX_ADDED`] points, //! which only a ring of thousands of long edges uses up; past that, its edges are left whole. //! -//! Arcs use a 1.8° step: 200 edges to a point's circle and about 100 to a semicircular cap. -//! Between vertices, those edges cut at most 0.013% of the distance inside the circle in the -//! projection plane. The buffer backend rounds a corner's segment count down, so a short corner -//! can have a step approaching 3.6° and cut up to 0.05% inside, on top of the projection error. +//! Arcs use a 2.8125° step: 128 edges to a point's circle and 64 to a semicircular cap. +//! Between vertices, those edges cut at most 0.031% of the distance inside the circle in the +//! projection plane. A short corner can have a step approaching 5.625° when the backend rounds +//! its segment count down, cutting at most 0.125% inside, on top of the projection error. use geo::{ - BooleanOps, Buffer, Contains, Coord as GC, LineString, MultiLineString, MultiPolygon, Orient, Point, Polygon, + BooleanOps, Buffer, Contains, Coord as GC, CoordsIter, LineString, MultiLineString, MultiPolygon, Orient, Point, Polygon, algorithm::buffer::{BufferStyle, LineCap, LineJoin}, orient::Direction, }; @@ -54,9 +54,11 @@ const MAX_SPLITS: u32 = 12; /// checked, so an arc's short edges cost nothing. const STRAY_PER_M2: f64 = 3e-7; -/// The finest round-cap/join step supported by the buffer backend (1.8°). A circle's area -/// deficit is under 0.017%, leaving room for projection error in the reference's 0.1% bound. -const ARC_STEP: f64 = std::f64::consts::PI / 100.0; +/// A circle's area deficit is under 0.041% with 128 edges: less than half the reference's +/// 0.1% bound, leaving room for projection error without the cost of the backend's finest step. +const ARC_STEP: f64 = std::f64::consts::PI / 64.0; + +pub(crate) const VERTEX_BUDGET_ERROR: &str = "Buffer exceeds its vertex budget; select fewer features"; /// Round arcs at the same density for points, line caps and polygon corners, including stretches. fn buffer_style(metres: f64) -> BufferStyle { @@ -66,16 +68,39 @@ fn buffer_style(metres: f64) -> BufferStyle { /// A feature's buffer at `metres` (negative shrinks a polygon), in lon/lat; `None` when nothing /// is left. pub fn geodesic_buffer(g: &Geometry, metres: f64) -> Option { + geodesic_buffer_with_limit(g, metres, usize::MAX).ok().flatten() +} + +/// A buffer that must fit the caller's remaining vertex budget. Multipoint discs are counted +/// before collecting them for union, so even overlapping discs cannot grow an unbounded batch. +/// Other geometries are checked before conversion to the returned geometry; this bounds retained +/// outlines, not all temporary allocations inside the polygon overlay backend. +pub(crate) fn geodesic_buffer_with_limit(g: &Geometry, metres: f64, max_vertices: usize) -> Result, &'static str> { if !metres.is_finite() { - return None; + return Ok(None); } let mp = match g { Geometry::Point(c) => point(*c, metres), - Geometry::MultiPoint(v) => union(v.iter().map(|c| point(*c, metres)).collect()), + Geometry::MultiPoint(v) => { + let mut remaining = max_vertices; + let mut discs = Vec::new(); + for c in v { + let disc = point(*c, metres); + if disc.0.is_empty() { + continue; + } + remaining = remaining.checked_sub(disc.coords_count()).ok_or(VERTEX_BUDGET_ERROR)?; + discs.push(disc); + } + union(discs) + } Geometry::Line(parts) => lines(&finite(parts), metres), Geometry::Polygon(polys) => polygons(&polys.iter().map(|rings| finite(rings)).collect::>(), metres), }; - to_geometry(&mp) + if mp.coords_count() > max_vertices { + return Err(VERTEX_BUDGET_ERROR); + } + Ok(to_geometry(&mp)) } /// An azimuthal equidistant plane whose centre's longitude is within 180° of `reference`, a @@ -351,6 +376,32 @@ fn to_geometry(mp: &MultiPolygon) -> Option { mod tests { use super::*; + /// Samples the actual projected outlines: shallow turns exercise a rounded-down join, + /// while larger turns would fail this bound at the old default arc density. + #[test] + fn round_polygon_corners_and_line_joins_keep_their_error_bound() { + for degrees in [0.25_f64, 2.8, 2.9, 5.5, 8.3, 89.0, 179.0] { + let angle = degrees.to_radians(); + let line = LineString::from(vec![(-100.0, 0.0), (0.0, 0.0), (100.0 * angle.cos(), 100.0 * angle.sin())]); + for shape in [geo::Geometry::LineString(line.clone()), geo::Geometry::Polygon(Polygon::new(line, vec![]))] { + let buffered = shape.buffer_with_style(buffer_style(1.0)); + let mut probes = 0; + for edge in buffered.0.iter().flat_map(|p| p.exterior().0.windows(2)) { + let [a, b] = edge else { continue }; + if (a.x.hypot(a.y) - 1.0).abs() > 1e-5 || (b.x.hypot(b.y) - 1.0).abs() > 1e-5 { + continue; + } + for t in [0.0, 0.25, 0.5, 0.75, 1.0] { + let distance = (a.x + t * (b.x - a.x)).hypot(a.y + t * (b.y - a.y)); + assert!((1.0 - 0.00125 - 1e-5..=1.0 + 1e-5).contains(&distance), "{degrees}°: {distance}"); + probes += 1; + } + } + assert!(probes > 0, "no corner samples at {degrees}°"); + } + } + } + #[test] fn a_ring_gains_at_most_max_added_points() { // 2,000 edges of about 100 km zigzagging across the plane, each needing dozens of splits. diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index 61402dd45..89a38257a 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -9,7 +9,7 @@ use std::borrow::Cow; use std::collections::{HashMap, HashSet, TryReserveError}; use crate::convert::{from_line, from_multipolygon, from_point, to_geo, to_multiline, to_multipolygon}; -use crate::geodesic::geodesic_buffer; +use crate::geodesic::{VERTEX_BUDGET_ERROR, geodesic_buffer, geodesic_buffer_with_limit}; use crate::{Inputs, Output, OutputLayer, ParamKind as K, ParamSpec as P, ToolResult, ToolSpec}; const UNITS: &[&str] = &["Kilometers", "Meters", "Miles", "Feet"]; @@ -1137,7 +1137,18 @@ fn union_all<'a>(polys: impl IntoIterator>) -> Mult // --------------------------------------------------------------------------------------------- // tools +/// At most 32 MB of retained outline coordinates, before dataset/overlay overhead. The same +/// budget applies before dissolve, so a large dissolved request cannot accumulate every disc. +const MAX_BUFFER_VERTICES: usize = 2_000_000; + +/// Buffer with the output budget shared by all selected input features. fn buffer(inp: &dyn Inputs, p: &Json) -> Result { + buffer_with_vertex_limit(inp, p, MAX_BUFFER_VERTICES) +} + +/// Runs the tool with a vertex budget; the explicit limit also lets the regression exercise +/// exhaustion with a few features rather than allocating millions of coordinates. +fn buffer_with_vertex_limit(inp: &dyn Inputs, p: &Json, limit: usize) -> Result { let ds = input(inp, p, "input")?; let d = num(p, "distance").ok_or("Distance must be a number")?; let metres = to_metres(d, s(p, "units"))?; @@ -1146,10 +1157,12 @@ fn buffer(inp: &dyn Inputs, p: &Json) -> Result { } let mut out = copy_schema(ds, out_name(p, format!("{}_Buffer", ds.name)), GeometryType::Polygon); let mut all: Vec> = Vec::new(); + let mut remaining = limit; for r in rows(inp, p, "input", ds) { let Some(f) = ds.features.get(r) else { continue }; let Some(g) = &f.geometry else { continue }; - let Some(geom) = geodesic_buffer(g, metres) else { continue }; + let Some(geom) = geodesic_buffer_with_limit(g, metres, remaining)? else { continue }; + remaining = remaining.saturating_sub(geom.coord_count()); if s(p, "dissolve") == "All" { all.extend(to_multipolygon(&geom)); } else { @@ -1160,6 +1173,9 @@ fn buffer(inp: &dyn Inputs, p: &Json) -> Result { out.fields.clear(); let u = union_all(all.iter().flat_map(|m| m.0.iter())); if let Some(geom) = from_multipolygon(&u) { + if geom.coord_count() > limit { + return Err(VERTEX_BUDGET_ERROR.into()); + } out.push(Some(geom), vec![]); } } @@ -4847,6 +4863,31 @@ mod tests { assert!(ran * 2 >= tested, "only {ran} of {tested} layer tools ran on these inputs"); } + /// The budget covers accumulated features and pre-union multipoint discs, including + /// dissolve requests. Exhaustion is an error, never a partial layer or a silent skip. + #[test] + fn buffer_vertex_budget_covers_features_and_multipoint_discs() { + let point = Coord::new(0.0, 0.0); + let vertices = geodesic_buffer(&Geometry::Point(point), 1000.0).unwrap().coord_count(); + let mut ds = Dataset::new("p", GeometryType::Point, vec![]); + for c in [point, Coord::new(0.5, 0.0)] { + ds.push(Some(Geometry::Point(c)), vec![]); + } + let mut m = Map(HashMap::from([("p".to_string(), ds)])); + for dissolve in ["None", "All"] { + let params = serde_json::json!({"input": "p", "distance": 1, "dissolve": dissolve}); + assert!(buffer_with_vertex_limit(&m, ¶ms, 2 * vertices).is_ok()); + assert!(buffer_with_vertex_limit(&m, ¶ms, vertices).unwrap_err().contains("vertex budget")); + } + // All discs overlap: the final union fits, but collecting the discs must still stop. + m.0.get_mut("p").unwrap().features.clear(); + m.0.get_mut("p").unwrap().push(Some(Geometry::MultiPoint(vec![point; 3])), vec![]); + for dissolve in ["None", "All"] { + let params = serde_json::json!({"input": "p", "distance": 1, "dissolve": dissolve}); + assert!(buffer_with_vertex_limit(&m, ¶ms, 2 * vertices).unwrap_err().contains("vertex budget")); + } + } + #[test] fn dissolved_buffers_merge() { let mut pts = Dataset::new("p", GeometryType::Point, vec![]); @@ -4991,7 +5032,7 @@ mod tests { /// 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.012% of the distance (or 1 cm) along a line's sides, and no further inside an arc - /// than 0.013% of the distance in the projection plane (#1221). The old 32-edge circles + /// 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. #[test] fn buffer_outlines_keep_their_distance_between_vertices() { @@ -5000,7 +5041,7 @@ mod tests { polys.iter().flatten().flatten().copied().chain(along_edges(polys, &fractions)).collect() }; // Bound the geometric error independently of how many vertices the buffer emits. - let arc_error = |d: f64| d * 1.3e-4; + let arc_error = |d: f64| d * 3.1e-4; let centre = Coord::new(10.125, 46.125); for c in points(&buffered(Geometry::Point(centre), 1.0)) { @@ -5034,7 +5075,7 @@ mod tests { 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.012% or 1 cm (and the sampler's millimetre); round the - // caps, add the 0.013% arc approximation in the projection plane. + // caps, add the 0.031% arc approximation in the projection plane. let side = (d0 * 1.2e-4).max(0.01) + 1e-3; let bound = if beside { side } else { arc_error(d0) + side }; assert!(d < d0 + side && d > d0 - bound, "{km} km: {c:?} is {d} m from the line"); diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 8dd1067ac..89a83e044 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -214,11 +214,15 @@ projection. Each feature is buffered in an azimuthal equidistant projection on t 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.012 % of the distance (or 1 cm), and within 0.05 % along -a winding line or outline buffered in stretches. Its curves are drawn with straight edges (200 -to a point's circle, about 100 to a semicircular line end). Between vertices, those edges cut up -to 0.013 % of the distance inside the true curve in the projection plane, and up to 0.05 % at +a winding line or outline buffered in stretches. Its curves are drawn with straight edges (128 +to a point's circle, 64 to a semicircular line end). Between vertices, those edges cut up +to 0.031 % of the distance inside the true curve in the projection plane, and up to 0.125 % at polygon corners and line joins, in addition to the projection error. Negative Buffer distances shrink polygons. Buffer and Select Layer By Location limit distances to 5,000 km. +Buffer limits retained outlines to 2,000,000 vertices before dissolving and checks the dissolved +result too. Multipoint discs count before union even when they overlap. A request that exceeds +the budget fails with an error; select fewer features and run smaller batches. This limits retained +geometry, not every temporary allocation inside a complex polygon operation. **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 regardless of ring winding, and multipart measurements are summed. Ring directions are