diff --git a/crates/analysis/src/geodesic.rs b/crates/analysis/src/geodesic.rs index 05e855665..784d58a05 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 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, orient::Direction, + BooleanOps, Buffer, Contains, Coord as GC, CoordsIter, 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}; @@ -65,19 +67,63 @@ const EARTH_R_M: f64 = rustgis_geom::EARTH_RADIUS_M; /// a tenth of the turn at 45°, as at about 5.7° of latitude. const FLAT_ALLOWANCE: f64 = 0.1; +/// 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. +const ARC_STEP: f64 = std::f64::consts::PI / 64.0; + +pub(crate) const VERTEX_BUDGET_ERROR: &str = "Buffer exceeds its vertex budget; split large features or select fewer features"; + +/// 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 { + geodesic_buffer_with_limit(g, metres, usize::MAX).ok().and_then(|b| b.geometry) +} + +/// A buffered outline and its budget charge: at least the retained coordinate count, or the +/// pre-union disc total for a multipoint, whichever is larger. +pub(crate) struct BudgetedBuffer { + pub geometry: Option, + pub charged_vertices: usize, +} + +/// A buffer that must fit the caller's remaining vertex budget. Multipoint discs are counted +/// before collecting them for union, and their charge is returned for the caller to accumulate. +/// 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 { if !metres.is_finite() { - return None; + return Ok(BudgetedBuffer { geometry: None, charged_vertices: 0 }); } + let mut pre_union_vertices = 0; 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); + } + pre_union_vertices = max_vertices - remaining; + 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); + } + let charged_vertices = pre_union_vertices.max(mp.coords_count()); + Ok(BudgetedBuffer { geometry: to_geometry(&mp), charged_vertices }) } /// Buffers of this many features or more are made on several threads. @@ -91,32 +137,51 @@ const STACK_BYTES: usize = 8 << 20; /// The buffers of `geoms` at `metres`, in their order, as [`geodesic_buffer`] makes them: on /// several threads when there are many (one in the browser), since each is made alone. pub fn geodesic_buffers(geoms: &[&Geometry], metres: f64) -> Vec> { + geodesic_buffers_with_limit(geoms, metres, usize::MAX).unwrap_or_default() +} + +/// Ordered parallel buffers sharing an outline budget, including pre-union multipoint discs. +/// `usize::MAX` disables cumulative charging for the public unlimited API. Each worker checks its current +/// allowance before buffering, then atomically charges the result before retaining it. In-flight +/// operations may allocate temporary geometry on each worker (one worker in the browser). +pub(crate) fn geodesic_buffers_with_limit(geoms: &[&Geometry], metres: f64, limit: usize) -> Result>, &'static str> { let threads = if cfg!(target_arch = "wasm32") || geoms.len() < PARALLEL_FEATURES { 1 } else { std::thread::available_parallelism().map_or(1, |n| n.get()).min(16) }; - buffers_on(geoms, metres, threads) + buffers_on(geoms, metres, threads, limit) } -/// [`geodesic_buffers`] on `threads` threads, this one among them. Each takes the next -/// [`BATCH`] of features until none are left, so the work is shared however it is spread, and -/// the results go back in their places. A thread that can't start (a limit on threads, say) -/// leaves its share to the others; one that panics passes its panic on, as a buffer made here -/// would. -fn buffers_on(geoms: &[&Geometry], metres: f64, threads: usize) -> Vec> { - if threads <= 1 { - return geoms.iter().map(|g| geodesic_buffer(g, metres)).collect(); - } - let next = std::sync::atomic::AtomicUsize::new(0); +/// [`geodesic_buffers_with_limit`] on `threads` threads, this one among them. Each takes the next +/// [`BATCH`] of features until none are left, and results go back in their places. The shared +/// remaining budget is charged before a worker retains an outline. A thread that cannot start +/// leaves its work to the others; a worker panic is passed on as in a buffer made here. +fn buffers_on(geoms: &[&Geometry], metres: f64, threads: usize, limit: usize) -> Result>, &'static str> { + use std::sync::atomic::{AtomicUsize, Ordering}; + let next = AtomicUsize::new(0); + let remaining = AtomicUsize::new(limit); let work = || { let mut made = Vec::new(); loop { - let from = next.fetch_add(BATCH, std::sync::atomic::Ordering::Relaxed); + let from = next.fetch_add(BATCH, Ordering::Relaxed); let Some(batch) = geoms.get(from..geoms.len().min(from + BATCH)).filter(|b| !b.is_empty()) else { break }; - made.extend(batch.iter().enumerate().map(|(i, g)| (from + i, geodesic_buffer(g, metres)))); + for (i, g) in batch.iter().enumerate() { + let buffered = geodesic_buffer_with_limit(g, metres, remaining.load(Ordering::Relaxed))?; + if limit != usize::MAX { + let mut left = remaining.load(Ordering::Relaxed); + loop { + let next = left.checked_sub(buffered.charged_vertices).ok_or(VERTEX_BUDGET_ERROR)?; + match remaining.compare_exchange_weak(left, next, Ordering::Relaxed, Ordering::Relaxed) { + Ok(_) => break, + Err(current) => left = current, + } + } + } + made.push((from + i, buffered.geometry)); + } } - made + Ok::<_, &'static str>(made) }; let mut out: Vec> = vec![None; geoms.len()]; let mut place = |made: Vec<(usize, Option)>| { @@ -126,17 +191,22 @@ fn buffers_on(geoms: &[&Geometry], metres: f64, threads: usize) -> Vec = - (1..threads).filter_map(|_| std::thread::Builder::new().stack_size(STACK_BYTES).spawn_scoped(scope, work).ok()).collect(); - place(work()); - for h in others { - // Not a new panic but the worker's, raised again here as a serial run would raise it, - // for `analysis::run` to catch (AGENTS.md's never-crash rule is about new ones). - place(h.join().unwrap_or_else(|e| std::panic::resume_unwind(e))); - } - }); - out + if threads <= 1 { + place(work()?); + } else { + std::thread::scope(|scope| { + let others: Vec<_> = + (1..threads).filter_map(|_| std::thread::Builder::new().stack_size(STACK_BYTES).spawn_scoped(scope, work).ok()).collect(); + place(work()?); + for h in others { + // Not a new panic but the worker's, raised again here as a serial run would raise it, + // for `analysis::run` to catch (AGENTS.md's never-crash rule is about new ones). + place(h.join().unwrap_or_else(|e| std::panic::resume_unwind(e))?); + } + Ok::<_, &'static str>(()) + })?; + } + Ok(out) } /// An azimuthal equidistant plane whose centre's longitude is within 180° of `reference`, a @@ -329,7 +399,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. @@ -338,7 +408,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)) } @@ -362,7 +433,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 @@ -397,7 +468,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() { @@ -458,6 +529,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. @@ -513,11 +610,23 @@ mod tests { let refs: Vec<&Geometry> = points.iter().collect(); let alone: Vec> = refs.iter().map(|g| geodesic_buffer(g, 5_000.0)).collect(); assert_eq!(geodesic_buffers(&refs, 5_000.0), alone); - // However many threads, and fewer features than threads or batches. + let vertices: usize = alone.iter().flatten().map(Geometry::coord_count).sum(); + // The exact shared limit succeeds; one vertex less fails on every thread count. + for threads in [1, 3, 16] { + assert_eq!(buffers_on(&refs, 5_000.0, threads, vertices).unwrap(), alone); + assert_eq!(buffers_on(&refs, 5_000.0, threads, vertices - 1).unwrap_err(), VERTEX_BUDGET_ERROR); + assert_eq!(buffers_on(&refs, 5_000.0, threads, usize::MAX).unwrap(), alone, "{threads} threads"); + assert_eq!(buffers_on(&refs[..3], 5_000.0, threads, usize::MAX).unwrap(), alone[..3], "{threads} threads, 3 features"); + } + assert!(buffers_on(&[], 5_000.0, 3, 0).unwrap().is_empty()); + // Overlapping multipoint unions stay small, but every disc is charged across workers. + let multipoint = Geometry::MultiPoint(vec![Coord::new(0.0, 0.0); 2]); + let refs = vec![&multipoint; 2 * BATCH]; + let disc_vertices = geodesic_buffer(&Geometry::Point(Coord::new(0.0, 0.0)), 5_000.0).unwrap().coord_count(); + let charge = refs.len() * 2 * disc_vertices; for threads in [1, 3, 16] { - assert_eq!(buffers_on(&refs, 5_000.0, threads), alone, "{threads} threads"); - assert_eq!(buffers_on(&refs[..3], 5_000.0, threads), alone[..3], "{threads} threads, 3 features"); + assert!(buffers_on(&refs, 5_000.0, threads, charge).is_ok()); + assert_eq!(buffers_on(&refs, 5_000.0, threads, charge - 1).unwrap_err(), VERTEX_BUDGET_ERROR); } - assert!(buffers_on(&[], 5_000.0, 3).is_empty()); } } diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index f639fc8fb..135eb3d58 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -1,6 +1,8 @@ //! The tools. -use geo::{BooleanOps, Contains, ConvexHull, InteriorPoint, Intersects, MultiPolygon, Orient, Relate, Simplify, orient::Direction}; +use geo::{ + BooleanOps, Contains, ConvexHull, CoordsIter, InteriorPoint, Intersects, MultiPolygon, Orient, Relate, Simplify, orient::Direction, +}; use rustgis_data::raster::Crs; use rustgis_data::{Coord, Dataset, Feature, Field, FieldType, Geometry, GeometryType, Value}; use rustgis_geom::Extent; @@ -9,7 +11,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, geodesic_buffers}; +use crate::geodesic::{VERTEX_BUDGET_ERROR, geodesic_buffer_with_limit, geodesic_buffers_with_limit}; use crate::{Inputs, Output, OutputLayer, ParamKind as K, ParamSpec as P, ToolResult, ToolSpec}; const UNITS: &[&str] = &["Kilometers", "Meters", "Miles", "Feet"]; @@ -1137,7 +1139,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"))?; @@ -1149,7 +1162,7 @@ fn buffer(inp: &dyn Inputs, p: &Json) -> Result { let picked: Vec<&Feature> = rows(inp, p, "input", ds).into_iter().filter_map(|r| ds.features.get(r)).collect(); let geoms: Vec<(&Feature, &Geometry)> = picked.iter().filter_map(|f| Some((*f, f.geometry.as_ref()?))).collect(); let without = picked.len() - geoms.len(); - let buffers = geodesic_buffers(&geoms.iter().map(|(_, g)| *g).collect::>(), metres); + let buffers = geodesic_buffers_with_limit(&geoms.iter().map(|(_, g)| *g).collect::>(), metres, limit)?; // Buffers with nothing in them: points at a negative distance, polygons that shrink away, // features with no finite coordinates. let mut empty = 0; @@ -1167,6 +1180,9 @@ fn buffer(inp: &dyn Inputs, p: &Json) -> Result { if !all.is_empty() { out.fields.clear(); let u = union_all(all.iter().flat_map(|m| m.0.iter())); + if u.coords_count() > limit { + return Err(VERTEX_BUDGET_ERROR.into()); + } if let Some(geom) = from_multipolygon(&u) { out.push(Some(geom), vec![]); } @@ -2183,22 +2199,33 @@ fn simplify(inp: &dyn Inputs, p: &Json) -> Result { /// Select against nearby selectors, indexed in the projection used by the predicate and after /// any distance buffer. Inversion and input row order keep the exhaustive selection's semantics. fn select_by_location(inp: &dyn Inputs, p: &Json) -> Result { + select_by_location_with_vertex_limit(inp, p, MAX_BUFFER_VERTICES) +} + +/// Bounds retained distance-buffer outlines across selectors before building the spatial index. +/// Zero-distance selection uses original geometries and does not consume the buffer budget. +fn select_by_location_with_vertex_limit(inp: &dyn Inputs, p: &Json, limit: usize) -> Result { let ds = input(inp, p, "input")?; let sel_ds = input(inp, p, "select")?; let dist = to_metres(num(p, "distance").unwrap_or(0.0).max(0.0), s(p, "units"))?; if dist > 5_000_000.0 { return Err("Distance is too large (limit 5000 km)".into()); } - let selectors: Vec<(geo::Geometry, Extent)> = rows(inp, p, "select", sel_ds) - .iter() - .filter_map(|r| sel_ds.features.get(*r)?.geometry.as_ref()) - .map(|g| { - let b = if dist > 0.0 { geodesic_buffer(g, dist) } else { None }; - let gg = to_geo(b.as_ref().unwrap_or(g)); - let ext = merc_extent(&gg); - (gg, ext) - }) - .collect(); + let mut remaining = limit; + let mut selectors = Vec::new(); + for r in rows(inp, p, "select", sel_ds) { + let Some(g) = sel_ds.features.get(r).and_then(|f| f.geometry.as_ref()) else { continue }; + let b = if dist > 0.0 { + let buffered = geodesic_buffer_with_limit(g, dist, remaining)?; + remaining = remaining.saturating_sub(buffered.charged_vertices); + buffered.geometry + } else { + None + }; + let gg = to_geo(b.as_ref().unwrap_or(g)); + let ext = merc_extent(&gg); + selectors.push((gg, ext)); + } let index = ExtentIndex::new(selectors.iter().map(|(_, ext)| *ext)); let invert = s(p, "invert") == "Yes"; let hits: Vec = (0..ds.features.len()) @@ -4207,6 +4234,7 @@ fn join_field_tool(inp: &dyn Inputs, p: &Json) -> Result { #[cfg(test)] mod tests { use super::*; + use crate::geodesic::geodesic_buffer; use rustgis_data::Feature; use std::collections::HashMap; @@ -4864,6 +4892,70 @@ 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")); + } + // Crossing thin rectangles add intersection vertices during dissolve even though no + // original vertex lies inside the other rectangle. The inputs fit; the union does not. + let rectangle = |x: f64, y: f64| { + Geometry::Polygon(vec![vec![vec![ + Coord::new(-x, -y), + Coord::new(x, -y), + Coord::new(x, y), + Coord::new(-x, y), + Coord::new(-x, -y), + ]]]) + }; + let mut cross = Dataset::new("cross", GeometryType::Polygon, vec![]); + for g in [rectangle(0.02, 0.001), rectangle(0.001, 0.02)] { + cross.push(Some(g), vec![]); + } + let outline_vertices = cross + .features + .iter() + .filter_map(|f| f.geometry.as_ref()) + .filter_map(|g| geodesic_buffer(g, 0.0)) + .map(|g| g.coord_count()) + .sum(); + let cross_inputs = Map(HashMap::from([("cross".into(), cross)])); + let mut params = serde_json::json!({"input":"cross", "distance":0, "dissolve":"None"}); + assert!(buffer_with_vertex_limit(&cross_inputs, ¶ms, outline_vertices).is_ok()); + params["dissolve"] = serde_json::json!("All"); + let union = layer(buffer_with_vertex_limit(&cross_inputs, ¶ms, usize::MAX).unwrap()); + assert!(union.features[0].geometry.as_ref().unwrap().coord_count() > outline_vertices); + assert!(buffer_with_vertex_limit(&cross_inputs, ¶ms, outline_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")); + } + // Each feature's discs fit alone, but their combined pre-union charge does not. + m.0.get_mut("p").unwrap().features.clear(); + for _ in 0..2 { + m.0.get_mut("p").unwrap().push(Some(Geometry::MultiPoint(vec![point; 2])), vec![]); + } + for dissolve in ["None", "All"] { + let params = serde_json::json!({"input":"p", "distance":1, "dissolve":dissolve}); + assert!(buffer_with_vertex_limit(&m, ¶ms, 4 * vertices).is_ok()); + assert!(buffer_with_vertex_limit(&m, ¶ms, 3 * vertices).unwrap_err().contains("vertex budget")); + } + } + #[test] fn buffers_made_on_threads_keep_their_features_and_say_which_had_no_geometry() { // Enough features for several threads: each output keeps its own attributes and buffer. @@ -5040,21 +5132,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.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() { 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 * 3.1e-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 @@ -5083,9 +5175,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.031% 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"); } } @@ -5096,7 +5188,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. @@ -6269,6 +6361,41 @@ mod tests { assert_eq!(out.features.iter().filter_map(|f| f.geometry.clone()).collect::>(), expected); } + /// Distance selection shares one outline budget, including overlapping multipoint discs. + /// Exhaustion must fail before producing a partial selection, even when inverted. + #[test] + fn location_selection_bounds_accumulated_buffer_outlines() { + let point = Coord::new(0.0, 0.0); + let vertices = geodesic_buffer(&Geometry::Point(point), 1000.0).unwrap().coord_count(); + let mut selectors = Dataset::new("selectors", GeometryType::Point, vec![]); + for c in [point, Coord::new(0.5, 0.0)] { + selectors.push(Some(Geometry::Point(c)), vec![]); + } + let mut features = Dataset::new("features", GeometryType::Point, vec![]); + features.push(Some(Geometry::Point(point)), vec![]); + let mut m = Map(HashMap::from([("selectors".into(), selectors), ("features".into(), features)])); + for invert in ["No", "Yes"] { + let p = serde_json::json!({"input":"features", "select":"selectors", "distance":1, "invert":invert}); + let result = select_by_location_with_vertex_limit(&m, &p, 2 * vertices).unwrap(); + let Output::Selection { rows, .. } = result.output else { panic!("expected selection") }; + assert_eq!(rows, if invert == "No" { vec![0] } else { vec![] }); + assert!(select_by_location_with_vertex_limit(&m, &p, vertices).unwrap_err().contains("vertex budget")); + } + // A zero-distance predicate does not allocate buffer outlines. + let p = serde_json::json!({"input":"features", "select":"selectors", "distance":0}); + assert!(select_by_location_with_vertex_limit(&m, &p, 0).is_ok()); + m.0.get_mut("selectors").unwrap().features.clear(); + m.0.get_mut("selectors").unwrap().push(Some(Geometry::MultiPoint(vec![point; 3])), vec![]); + let p = serde_json::json!({"input":"features", "select":"selectors", "distance":1}); + assert!(select_by_location_with_vertex_limit(&m, &p, 2 * vertices).unwrap_err().contains("vertex budget")); + m.0.get_mut("selectors").unwrap().features.clear(); + for _ in 0..2 { + m.0.get_mut("selectors").unwrap().push(Some(Geometry::MultiPoint(vec![point; 2])), vec![]); + } + assert!(select_by_location_with_vertex_limit(&m, &p, 4 * vertices).is_ok()); + assert!(select_by_location_with_vertex_limit(&m, &p, 3 * vertices).unwrap_err().contains("vertex budget")); + } + #[test] fn indexed_selection_matches_every_selector_including_buffers() { for selection in [None, Some((0..80).filter(|i| i % 3 != 0).collect()), Some(vec![83])] { 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 e09b1885f..01dfe1491 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -214,11 +214,17 @@ 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 (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, cumulatively across features. A request that exceeds +the budget fails with an error; split large features or select fewer features and run smaller batches. This limits retained +geometry, not every temporary allocation inside a complex polygon operation. Select Layer By +Location applies the same shared budget to its distance-buffered selector outlines; selection at +zero distance uses the original geometries without this buffer limit. **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