From d7c685810d1d8800f82e2ee0ba2e453970362b0a Mon Sep 17 00:00:00 2001 From: Dominik Strebel Date: Tue, 3 Nov 2026 13:08:46 +0100 Subject: [PATCH 1/3] perf(analysis): Zonal Statistics works out its zones on several cores (#1252) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Each zone's statistics are its own, its cells summed in the same order on any thread, so the zones run on parallel::try_map_on and the rows are then made in input order: the same values, bit for bit (zonal_stats_on_several_threads_is_zonal_stats_on_one). A thread holds a zone's mask, up to the raster's size, so the threads are also bounded by the scratch budget (parallel::threads_for_scratch). The tool now honours Cancel and reports progress by zone. 20,000 small zones over a 5,000 × 5,000 raster: 0.40 s on one thread, 0.07 s on 16. 2,000 zones of 300 × 300 cells: 3.33 s on one thread, 0.40 s on 11 (the masks' memory allows 11). --- crates/analysis/src/parallel.rs | 6 ++ crates/analysis/src/tools.rs | 111 ++++++++++++++++++++++++-------- 2 files changed, 89 insertions(+), 28 deletions(-) diff --git a/crates/analysis/src/parallel.rs b/crates/analysis/src/parallel.rs index 85408220b..79e74ad8c 100644 --- a/crates/analysis/src/parallel.rs +++ b/crates/analysis/src/parallel.rs @@ -30,6 +30,12 @@ pub(crate) fn threads_for(items: usize, min_items: usize) -> usize { cores().min(items.div_ceil(BATCH)).max(1) } +/// [`threads_for`], each thread holding `scratch_bytes` at once (a feature's raster mask, say): no +/// more than [`SCRATCH_BUDGET`] holds. +pub(crate) fn threads_for_scratch(items: usize, min_items: usize, scratch_bytes: usize) -> usize { + threads_for(items, min_items).min(1 + SCRATCH_BUDGET / scratch_bytes.max(1)) +} + /// How many threads to make `rows` raster rows on, `work` in all, each thread taking /// `scratch_bytes` of scratch: one below [`MIN_ROW_WORK`] and in the browser, no more than there /// are rows, and no more than [`SCRATCH_BUDGET`] holds (a scratch per point, for millions of diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index a387546d9..a6efc8c58 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -4943,6 +4943,11 @@ fn parse_stats(p: &Json) -> Result, String> { } fn zonal_stats(inp: &dyn Inputs, p: &Json) -> Result { + zonal_stats_on(inp, p, None) +} + +/// [`zonal_stats`] on `threads` threads (none: as many as pay off and the masks' memory allows). +fn zonal_stats_on(inp: &dyn Inputs, p: &Json, threads: Option) -> Result { let r_name = { let n = s(p, "raster"); if n.is_empty() { s(p, "layer") } else { n } @@ -5007,38 +5012,47 @@ fn zonal_stats(inp: &dyn Inputs, p: &Json) -> Result { out.fields.push(Field::new(name.clone(), *ty)); } - let target_rows = rows(inp, p, poly_param, poly_ds); + let target_rows: Vec = rows(inp, p, poly_param, poly_ds).into_iter().filter(|&r| r < poly_ds.features.len()).collect(); + // Each feature's statistics are its own, its cells summed in the same order on any thread, so + // the features are worked out on several threads when there are many (#1252). A thread holds + // a feature's mask, up to the raster's size, so their memory bounds the threads too. + let mask_bytes = (r.width as usize).saturating_mul(r.height as usize); + let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_scratch(target_rows.len(), PARALLEL_FEATURES, mask_bytes)); + crate::control::begin(target_rows.len()); + // A feature's statistics by band, whether its geometry was read, its parts that couldn't be + // projected and what went wrong. + let worked = crate::parallel::try_map_on( + &target_rows, + threads, + || (), + |_, _, &r_idx| { + let none = || vec![crate::rasterops::ZoneStats::default(); bands.len()]; + let Some(f) = poly_ds.features.get(r_idx) else { return Ok((none(), false, 0, None)) }; + Ok::<_, String>(match &f.geometry { + Some(Geometry::Polygon(polys)) if !polys.is_empty() => { + // Dataset geometry is already lon/lat; crs records the source system. + match crate::rasterops::transform_feature_polygons_to_lonlat(polys, None) { + Ok(polys_lonlat) => { + let mut failed_parts = Vec::new(); + let st = crate::rasterops::zonal_statistics_with_failures(r, &bands, &polys_lonlat, &mut failed_parts)?; + let detail = (!failed_parts.is_empty()) + .then(|| format!("feature {}: raster projection failed for polygon part(s) {failed_parts:?}", f.id)); + (st, true, failed_parts.len(), detail) + } + Err(e) => (none(), false, 0, Some(format!("feature {}: {e}", f.id))), + } + } + _ => (none(), false, 0, Some(format!("feature {}: missing or empty geometry", f.id))), + }) + }, + )?; let mut skipped_features = 0usize; let mut skipped_parts = 0usize; let mut failure_details = Vec::new(); - for r_idx in target_rows { + for (&r_idx, (stats_by_band, feature_valid, failed_parts, detail)) in target_rows.iter().zip(worked) { let Some(f) = poly_ds.features.get(r_idx) else { continue }; - let (stats_by_band, feature_valid) = match &f.geometry { - Some(Geometry::Polygon(polys)) if !polys.is_empty() => { - // Dataset geometry is already lon/lat; crs records the source system. - match crate::rasterops::transform_feature_polygons_to_lonlat(polys, None) { - Ok(polys_lonlat) => { - let mut failed_parts = Vec::new(); - let st = crate::rasterops::zonal_statistics_with_failures(r, &bands, &polys_lonlat, &mut failed_parts)?; - skipped_parts += failed_parts.len(); - if !failed_parts.is_empty() { - failure_details - .push(format!("feature {}: raster projection failed for polygon part(s) {failed_parts:?}", f.id)); - } - (st, true) - } - Err(e) => { - failure_details.push(format!("feature {}: {e}", f.id)); - (vec![crate::rasterops::ZoneStats::default(); bands.len()], false) - } - } - } - _ => { - failure_details.push(format!("feature {}: missing or empty geometry", f.id)); - (vec![crate::rasterops::ZoneStats::default(); bands.len()], false) - } - }; - + skipped_parts += failed_parts; + failure_details.extend(detail); if !feature_valid { skipped_features += 1; } @@ -9641,6 +9655,47 @@ mod tests { assert!(crate::param_applies("buffer", "units", &with("", ""))); } + /// Zones worked out on several threads give the same statistics, bit for bit, as on one, in + /// input order (#1252): zones without geometry, past the raster's edge and overlapping ones + /// included. + #[test] + fn zonal_stats_on_several_threads_is_zonal_stats_on_one() { + let (w, h) = (200u32, 150u32); + let band = |k: f32| (0..w * h).map(|i| ((i % w) as f32 * 0.37 + (i / w) as f32 * 1.13 + k).sin() * 1e3 + 0.1).collect::>(); + let gt = rustgis_data::raster::GeoTransform { origin_x: 0.0, origin_y: 15.0, pixel_width: 0.1, pixel_height: -0.1 }; + let raster = rustgis_data::raster::Raster::new( + "elev", + w, + h, + vec![band(0.0), band(2.0)], + gt, + rustgis_data::raster::Crs::LonLat, + None, + rustgis_data::raster::ColorHint::Gray, + "Float32", + ) + .unwrap(); + let mut zones = Dataset::new("zones", GeometryType::Polygon, vec![Field::new("id", FieldType::Integer)]); + for i in 0..600i64 { + let (x, y, d) = ((i * 37 % 230) as f64 / 10.0 - 1.0, (i * 53 % 170) as f64 / 10.0 - 1.0, 0.3 + (i % 7) as f64 * 0.4); + zones.push(if i % 61 == 0 { None } else { Some(sq(x, y, d)) }, vec![Value::Int(i)]); + } + let m = + RasterMap { datasets: HashMap::from([("zones".to_string(), zones)]), rasters: HashMap::from([("elev".to_string(), &raster)]) }; + let p = serde_json::json!({"raster": "elev", "polygons": "zones", "bands": [1, 2], "stats": ["count", "sum", "mean", "min", "max", "stddev"]}); + let one = zonal_stats_on(&m, &p, Some(1)).unwrap(); + let attrs = |r: &ToolResult| match &r.output { + Output::Layer(l) => l.features.iter().map(|f| format!("{:?}", f.attrs)).collect::>(), + _ => panic!("a layer"), + }; + assert!(attrs(&one).iter().any(|a| a.contains("Float")), "some zones meet the raster"); + for threads in [2, 7, 16] { + let many = zonal_stats_on(&m, &p, Some(threads)).unwrap(); + assert_eq!(attrs(&many), attrs(&one), "{threads} threads"); + assert_eq!(many.messages, one.messages, "{threads} threads"); + } + } + #[test] fn zonal_stats_tool_cases() { let (w, h) = (10u32, 10u32); From d261f3f31cd816043b77e5dbaf0faead9e66ecbb Mon Sep 17 00:00:00 2001 From: Dominik Strebel Date: Tue, 3 Nov 2026 13:47:03 +0100 Subject: [PATCH 2/3] perf(analysis): Zonal Statistics counts a zone's mask and its trimmed copy against the threads' memory MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit From the review of #1377: polygon_mask_with_failures keeps a zone's bounding-box mask while it makes the trimmed copy, so a thread holds up to two raster-sized masks at its peak. The threads are now bounded by two each (2,000 zones of 300 × 300 cells over a 5,000 × 5,000 raster: 6 threads, 0.61 s, against 3.20 s on one). --- crates/analysis/src/parallel.rs | 4 ++-- crates/analysis/src/tools.rs | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/crates/analysis/src/parallel.rs b/crates/analysis/src/parallel.rs index 79e74ad8c..e9e0eccb1 100644 --- a/crates/analysis/src/parallel.rs +++ b/crates/analysis/src/parallel.rs @@ -30,8 +30,8 @@ pub(crate) fn threads_for(items: usize, min_items: usize) -> usize { cores().min(items.div_ceil(BATCH)).max(1) } -/// [`threads_for`], each thread holding `scratch_bytes` at once (a feature's raster mask, say): no -/// more than [`SCRATCH_BUDGET`] holds. +/// [`threads_for`], each thread holding up to `scratch_bytes` at its peak (a feature's raster +/// masks, say): no more than [`SCRATCH_BUDGET`] holds. pub(crate) fn threads_for_scratch(items: usize, min_items: usize, scratch_bytes: usize) -> usize { threads_for(items, min_items).min(1 + SCRATCH_BUDGET / scratch_bytes.max(1)) } diff --git a/crates/analysis/src/tools.rs b/crates/analysis/src/tools.rs index a6efc8c58..6c413f9a4 100644 --- a/crates/analysis/src/tools.rs +++ b/crates/analysis/src/tools.rs @@ -5015,8 +5015,9 @@ fn zonal_stats_on(inp: &dyn Inputs, p: &Json, threads: Option) -> Result< let target_rows: Vec = rows(inp, p, poly_param, poly_ds).into_iter().filter(|&r| r < poly_ds.features.len()).collect(); // Each feature's statistics are its own, its cells summed in the same order on any thread, so // the features are worked out on several threads when there are many (#1252). A thread holds - // a feature's mask, up to the raster's size, so their memory bounds the threads too. - let mask_bytes = (r.width as usize).saturating_mul(r.height as usize); + // a feature's mask, up to the raster's size, and briefly its trimmed copy too, so their memory + // bounds the threads too. + let mask_bytes = (r.width as usize).saturating_mul(r.height as usize).saturating_mul(2); let threads = threads.unwrap_or_else(|| crate::parallel::threads_for_scratch(target_rows.len(), PARALLEL_FEATURES, mask_bytes)); crate::control::begin(target_rows.len()); // A feature's statistics by band, whether its geometry was read, its parts that couldn't be From 862adc0c2e501c551ef2b6dce95d1c3bd521443d Mon Sep 17 00:00:00 2001 From: Dominik Strebel Date: Tue, 3 Nov 2026 14:02:05 +0100 Subject: [PATCH 3/3] fix(analysis): try_map counts progress and looks at Cancel item by item on several threads too From the review of #1377: on several threads try_map_on counted a batch of 16 items at once and looked at Cancel only before each batch, so with items that take long (a large zone) the progress bar lagged and Cancel waited for up to 16 of them per thread. Both are now done after and before each item, as on one thread; the count is one shared atomic add an item. --- crates/analysis/src/parallel.rs | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/crates/analysis/src/parallel.rs b/crates/analysis/src/parallel.rs index e9e0eccb1..ed16bdbb2 100644 --- a/crates/analysis/src/parallel.rs +++ b/crates/analysis/src/parallel.rs @@ -166,11 +166,13 @@ pub(crate) fn try_map_on>( while !failed.load(Ordering::Relaxed) { let from = next.fetch_add(BATCH, Ordering::Relaxed); let Some(batch) = items.get(from..items.len().min(from.saturating_add(BATCH))).filter(|b| !b.is_empty()) else { break }; - if control::cancelled() { - failed.store(true, Ordering::Relaxed); - return Err((from, Cancelled.into())); - } + // Cancel is looked at, and progress counted, item by item, as on one thread: an item + // can take long (a large zone, say), so a batch of them would hold both back. for (k, item) in batch.iter().enumerate() { + if control::cancelled() { + failed.store(true, Ordering::Relaxed); + return Err((from + k, Cancelled.into())); + } match f(&mut s, from + k, item) { Ok(r) => made.push((from + k, r)), Err(e) => { @@ -178,9 +180,8 @@ pub(crate) fn try_map_on>( return Err((from + k, e)); } } + control::advance(1); } - // Counted a batch at a time: the threads share one count. - control::advance(batch.len()); } Ok(made) };