diff --git a/crates/data/src/raster/mod.rs b/crates/data/src/raster/mod.rs index cb2107d5f..faf7a78aa 100644 --- a/crates/data/src/raster/mod.rs +++ b/crates/data/src/raster/mod.rs @@ -43,8 +43,20 @@ pub fn filled(w: usize, h: usize, fill: T) -> Result> { Ok(v) } +/// A `w × h` band of `values`, its memory reserved first ([`grid`]). More than `w × h` values are +/// an error, not cut off. +pub fn collected(w: usize, h: usize, values: impl IntoIterator) -> Result> { + let mut v = grid(w, h)?; + let mut values = values.into_iter(); + v.extend(values.by_ref().take(w * h)); + if values.next().is_some() { + return Err(DataError::Format(format!("more values than a {w} × {h} raster holds"))); + } + Ok(v) +} + /// Room for `n` samples of a `w × h` band (fewer than all of them, when fewer can come). -pub(crate) fn room(n: usize, w: usize, h: usize) -> Result> { +pub fn room(n: usize, w: usize, h: usize) -> Result> { let mut v = Vec::new(); v.try_reserve_exact(n).map_err(|_| no_memory(w, h))?; Ok(v) @@ -1136,6 +1148,10 @@ pub(crate) mod tests { assert!(room::(usize::MAX / 2, 1, 1).is_err(), "too many bytes to ask for"); assert_eq!(filled(3, 2, f32::NAN).map(|v| v.len()).ok(), Some(6)); assert!(grid::(4, 5).is_ok_and(|g| g.is_empty() && g.capacity() >= 20)); + assert!(collected(huge, huge, std::iter::repeat(0.0f32)).is_err()); + assert!(collected(3, 2, 0..7).is_err(), "more than w × h values"); + assert_eq!(collected(3, 2, 0..6).ok(), Some(vec![0, 1, 2, 3, 4, 5])); + assert_eq!(collected(3, 2, 0..4).ok(), Some(vec![0, 1, 2, 3])); } /// A PNG of `w`×`h` pixels, red in its top half and blue below. diff --git a/crates/engine/src/cmd/cube.rs b/crates/engine/src/cmd/cube.rs index 87c7c4167..2821e5740 100644 --- a/crates/engine/src/cmd/cube.rs +++ b/crates/engine/src/cmd/cube.rs @@ -621,7 +621,7 @@ impl TrendRun { // first piece is worked out, so a run cancelled or refused before then hasn't taken it. if let (0, Some(total)) = (piece, whole) { for band in &mut bands { - band.reserve_exact(total); + reserve(band, total)?; } } for (k, band) in bands.iter_mut().enumerate() { @@ -684,7 +684,8 @@ impl TrendRun { if let Some(e) = too_long.and_then(|f| f(cells)) { return Err(e); } - block = vec![f32::NAN; cells * n]; + // Reserved first: a block the machine can't hold is an error, not a crash (#1291). + block = raster::filled(cells, n, f32::NAN).map_err(|e| e.to_string())?; grid = Some(r.clone()); } Some(g) => off_grid(self.cube.as_ref(), i, g, &r)?, @@ -723,7 +724,7 @@ impl TrendRun { done: &AtomicUsize, ) -> Result, String> { let n = self.steps.len(); - let mut out = vec![f32::NAN; cells * OUTPUTS]; + let mut out = raster::filled(cells, OUTPUTS, f32::NAN).map_err(|e| e.to_string())?; let stopped = AtomicBool::new(false); // Cancel is looked at after about the same work whatever the dates: a cell compares n²⁄2 // pairs, so every 4096 cells for a short series, every cell for a long one. @@ -1270,19 +1271,27 @@ pub(crate) struct Totals { } impl Totals { - pub(crate) fn new(stat: Stat, cells: usize) -> Totals { + /// Totals for `cells` cells, their memory reserved first: totals the machine can't hold are an + /// error, not a crash (#1291). + pub(crate) fn new(stat: Stat, cells: usize) -> Result { + fn each(cells: usize, v: T) -> Result, String> { + let mut out = Vec::new(); + reserve(&mut out, cells)?; + out.resize(cells, v); + Ok(out) + } let acc = match stat { Stat::Count | Stat::Mode | Stat::Variety | Stat::ModeCount | Stat::LongestRun | Stat::LongestRunLength | Stat::Frequency => { Vec::new() } - Stat::Min => vec![f64::INFINITY; cells], - Stat::Max => vec![f64::NEG_INFINITY; cells], - _ => vec![0.0; cells], + Stat::Min => each(cells, f64::INFINITY)?, + Stat::Max => each(cells, f64::NEG_INFINITY)?, + _ => each(cells, 0.0)?, }; - let m2 = if matches!(stat, Stat::Std) { vec![0.0; cells] } else { Vec::new() }; - let last = if stat.in_order() { vec![f32::NAN; cells] } else { Vec::new() }; - let (run, longest) = if stat.runs() { (vec![0; cells], vec![(f32::NAN, 0); cells]) } else { (Vec::new(), Vec::new()) }; - Totals { stat, n: vec![0; cells], acc, m2, last, run, longest, classes: Vec::new(), counts: Vec::new(), refused: None } + let m2 = if matches!(stat, Stat::Std) { each(cells, 0.0)? } else { Vec::new() }; + let last = if stat.in_order() { each(cells, f32::NAN)? } else { Vec::new() }; + let (run, longest) = if stat.runs() { (each(cells, 0)?, each(cells, (f32::NAN, 0))?) } else { (Vec::new(), Vec::new()) }; + Ok(Totals { stat, n: each(cells, 0)?, acc, m2, last, run, longest, classes: Vec::new(), counts: Vec::new(), refused: None }) } /// How many class counts it holds: its cells × the classes it has met. @@ -1577,7 +1586,11 @@ impl ReduceRun { } Some(first) => off_grid(self.cube.as_ref(), i, first, &r)?, } - let t = totals.get_mut(g).ok_or("a date in no group")?.get_or_insert_with(|| Totals::new(self.stat, vals.len())); + let slot = totals.get_mut(g).ok_or("a date in no group")?; + if slot.is_none() { + *slot = Some(Totals::new(self.stat, vals.len())?); + } + let t = slot.as_mut().ok_or("a date in no group")?; let nodata = r.nodata.map(|v| v as f32); for (c, v) in vals.iter().enumerate().take(t.cells()) { if v.is_finite() && Some(*v) != nodata { @@ -1623,14 +1636,14 @@ impl ReduceRun { let (columns, piece_mask) = match (finished.get_mut(g).and_then(Option::take), totals.get_mut(g).and_then(Option::take)) { (Some(out), _) => out, (None, Some(t)) => finish(&t, self.min_share, dates, frequency), - (None, None) => finish(&Totals::new(self.stat, cells), self.min_share, dates, frequency), + (None, None) => finish(&Totals::new(self.stat, cells)?, self.min_share, dates, frequency), }; // Their full length when the cube knows it, so they don't grow by doubling. - stack_columns(band, mask, columns, &piece_mask, whole); + stack_columns(band, mask, columns, &piece_mask, whole)?; if let (0, Some(total)) = (piece, whole) && frequency { - mask.reserve_exact(total); + reserve(mask, total)?; } mask.extend(piece_mask); } @@ -1706,7 +1719,7 @@ fn classes_over(totals: &[Option], finished: &[Option<(Columns, Vec /// Puts a piece's columns under a group's so far, by key: a class the piece didn't meet gets the /// piece's `piece_mask`, and one met first in this piece the pieces' before it (`mask`). A column /// made here takes its full length, `whole` cells, when the cube knows it. -fn stack_columns(outs: &mut Columns, mask: &[f32], columns: Columns, piece_mask: &[f32], whole: Option) { +fn stack_columns(outs: &mut Columns, mask: &[f32], columns: Columns, piece_mask: &[f32], whole: Option) -> Result<(), String> { let order = |a: &Option, b: &Option| match (a, b) { (Some(a), Some(b)) => a.total_cmp(b), _ => a.is_some().cmp(&b.is_some()), @@ -1720,7 +1733,8 @@ fn stack_columns(outs: &mut Columns, mask: &[f32], columns: Columns, piece_mask: } } Err(at) => { - let mut b = Vec::with_capacity(whole.unwrap_or(0).max(mask.len() + column.len())); + let mut b = Vec::new(); + reserve(&mut b, whole.unwrap_or(0).max(mask.len() + column.len()))?; b.extend_from_slice(mask); b.extend(column); outs.insert(at, (key, b)); @@ -1732,6 +1746,13 @@ fn stack_columns(outs: &mut Columns, mask: &[f32], columns: Columns, piece_mask: b.extend_from_slice(piece_mask); } } + Ok(()) +} + +/// Room for `cells` more values in `v`, reserved first: outputs the machine can't hold are an +/// error, not a crash (#1291). +fn reserve(v: &mut Vec, cells: usize) -> Result<(), String> { + v.try_reserve_exact(cells).map_err(|_| format!("not enough memory for the values of {cells} cells")) } /// Prepares `cube.reduce`: the cube, the statistic, and the dates of the period (and months) @@ -1934,16 +1955,16 @@ mod tests { let nan = f32::NAN; let mut outs = Vec::new(); let mut mask = Vec::new(); - super::stack_columns(&mut outs, &mask, vec![(Some(3.0), vec![1.0, nan])], &[0.0, nan], None); + super::stack_columns(&mut outs, &mask, vec![(Some(3.0), vec![1.0, nan])], &[0.0, nan], None).unwrap(); mask.extend([0.0, nan]); - super::stack_columns(&mut outs, &mask, vec![(Some(5.0), vec![4.0])], &[0.0], Some(3)); + super::stack_columns(&mut outs, &mask, vec![(Some(5.0), vec![4.0])], &[0.0], Some(3)).unwrap(); let got: Vec<(Option, String)> = outs.iter().map(|(k, b)| (*k, format!("{b:?}"))).collect(); assert_eq!(got, vec![(Some(3.0), "[1.0, NaN, 0.0]".to_string()), (Some(5.0), "[0.0, NaN, 4.0]".to_string())]); // One output, keyed None, is extended piece after piece. let mut one = Vec::new(); - super::stack_columns(&mut one, &[], vec![(None, vec![1.0])], &[], Some(3)); + super::stack_columns(&mut one, &[], vec![(None, vec![1.0])], &[], Some(3)).unwrap(); assert!(one.first().is_some_and(|(_, b)| b.capacity() >= 3), "made at its full length"); - super::stack_columns(&mut one, &[], vec![(None, vec![2.0, 3.0])], &[], Some(3)); + super::stack_columns(&mut one, &[], vec![(None, vec![2.0, 3.0])], &[], Some(3)).unwrap(); assert_eq!(one, vec![(None, vec![1.0, 2.0, 3.0])]); } @@ -1952,7 +1973,7 @@ mod tests { #[test] fn groups_being_read_share_the_class_cap() { let group = |classes: &[f64]| { - let mut t = super::Totals::new(super::Stat::Mode, 10); + let mut t = super::Totals::new(super::Stat::Mode, 10).unwrap(); for (c, v) in classes.iter().enumerate() { t.add(c, *v); } diff --git a/crates/engine/src/cmd/imagery.rs b/crates/engine/src/cmd/imagery.rs index 013962622..125162d41 100644 --- a/crates/engine/src/cmd/imagery.rs +++ b/crates/engine/src/cmd/imagery.rs @@ -7,7 +7,7 @@ use std::sync::Arc; use rustgis_analysis::georef::{self, Gcp, Kind}; use rustgis_analysis::rasterops::{self, NODATA}; use rustgis_analysis::warp; -use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster}; +use rustgis_data::raster::{self, ColorHint, Crs, GeoTransform, Raster}; use rustgis_data::{Geometry, GeometryType}; use rustgis_geom::{Coord, Extent}; use rustgis_map::RasterLayer; @@ -174,6 +174,22 @@ fn replace_layer(s: &mut Session, layer: &str, mut new: rustgis_map::RasterLayer Ok(json!({"layer": id})) } +/// A copy of `r`'s bands, whole, their memory reserved first: a raster the machine can't hold +/// twice is an error, not a crash (#1291). A band of the wrong length is copied as it is, for +/// `Raster::new` to refuse. +fn copied(r: &Raster) -> Result>, String> { + let (w, h) = (r.width as usize, r.height as usize); + r.bands + .iter() + .map(|b| { + let mut copy = raster::room(b.data.len(), w, h)?; + copy.extend_from_slice(&b.data); + Ok(copy) + }) + .collect::>() + .map_err(|e| e.to_string()) +} + fn rebuild(r: &Raster, transform: GeoTransform, crs: Crs, bands: Vec>) -> Result { let mut out = Raster::new(r.name.clone(), r.width, r.height, bands, transform, crs, r.nodata, r.color.clone(), r.data_type.clone()) .map_err(|e| e.to_string())?; @@ -238,7 +254,7 @@ fn georeference(s: &mut Session, p: &Value) -> CmdResult { } GeoTransform { origin_x: bx, origin_y: by, pixel_width: ax, pixel_height: ay } }; - let mut out = rebuild(&r, gt, Crs::LonLat, r.bands.iter().map(|b| b.data.clone()).collect())?; + let mut out = rebuild(&r, gt, Crs::LonLat, copied(&r)?)?; // The same pixels, so the same values rounded when the file was read. out.rounded = r.rounded.clone(); let e = out.extent(); @@ -545,24 +561,21 @@ fn change(s: &mut Session, p: &Value) -> CmdResult { return Err(format!("unknown method '{method}'")); } let band = uint(p, "band").unwrap_or(1).saturating_sub(1) as usize; - let av: Vec = a.bands.get(band).ok_or("no such band")?.data.clone(); + let av = &a.bands.get(band).ok_or("no such band")?.data; let bv = rasterops::resample_to(&bfr, band, &a)?; - let out: Vec = av - .iter() - .zip(&bv) - .map(|(x, y)| { - if !a.is_valid(*x) || *y == NODATA { - return NODATA; - } - let v = match method { - "ratio" if *x != 0.0 => y / x, - "percent" if *x != 0.0 => (y - x) / x * 100.0, - "difference" => y - x, - _ => f32::NAN, - }; - if v.is_finite() { v } else { NODATA } - }) - .collect(); + let values = av.iter().zip(&bv).map(|(x, y)| { + if !a.is_valid(*x) || *y == NODATA { + return NODATA; + } + let v = match method { + "ratio" if *x != 0.0 => y / x, + "percent" if *x != 0.0 => (y - x) / x * 100.0, + "difference" => y - x, + _ => f32::NAN, + }; + if v.is_finite() { v } else { NODATA } + }); + let out = raster::collected(a.width as usize, a.height as usize, values).map_err(|e| e.to_string())?; add(s, rasterops::derive(&a, &name(p, format!("{}_change", a.name)), vec![out])?, false) } @@ -798,13 +811,15 @@ fn process(s: &mut Session, p: &Value) -> CmdResult { let k = (uint(p, "factor").unwrap_or(2).clamp(2, 64) as u32).min(r.width.min(r.height)).max(1); let (w, h) = ((r.width / k).max(1), (r.height / k).max(1)); // Each output pixel is the mean of the valid pixels it covers. - let bands: Vec> = r + let bands = r .bands .iter() .map(|b| { - (0..h) - .flat_map(|row| (0..w).map(move |col| (row, col))) - .map(|(row, col)| { + let cells = (0..h).flat_map(|row| (0..w).map(move |col| (row, col))); + raster::collected( + w as usize, + h as usize, + cells.map(|(row, col)| { let (mut sum, mut n) = (0.0f64, 0u32); for dy in 0..k { for dx in 0..k { @@ -819,10 +834,11 @@ fn process(s: &mut Session, p: &Value) -> CmdResult { } } if n > 0 { (sum / n as f64) as f32 } else { NODATA } - }) - .collect() + }), + ) }) - .collect(); + .collect::, _>>() + .map_err(|e| e.to_string())?; let gt = GeoTransform { pixel_width: r.transform.pixel_width * k as f64, pixel_height: r.transform.pixel_height * k as f64, @@ -940,7 +956,10 @@ fn raster_mosaic(s: &mut Session, p: &Value) -> CmdResult { let (w, h) = (w as u32, h as u32); let nodata = first.nodata.or(Some(NODATA as f64)); let fill = nodata.map(|v| v as f32).unwrap_or(NODATA); - let mut bands: Vec> = vec![vec![fill; w as usize * h as usize]; first.band_count()]; + let mut bands = (0..first.band_count()) + .map(|_| raster::filled(w as usize, h as usize, fill)) + .collect::, _>>() + .map_err(|e| e.to_string())?; for (r, (pc, pr)) in rs.iter().zip(&places) { let (dc, dr) = ((pc - c0) as usize, (pr - r0) as usize); for (out, b) in bands.iter_mut().zip(&r.bands) { @@ -1046,20 +1065,22 @@ fn clip_window( let (w, h) = (cols.end.saturating_sub(cols.start), rows.end.saturating_sub(rows.start)); let nodata = if mask.is_some() { r.nodata.or(Some(NODATA as f64)) } else { r.nodata }; let fill = nodata.map(|v| v as f32).unwrap_or(NODATA); - let bands: Vec> = r + let bands = r .bands .iter() .map(|b| { - rows.clone() - .flat_map(|row| cols.clone().map(move |col| (row, col))) - .enumerate() - .map(|(i, (row, col))| match mask.map(|m| m.get(i).copied().unwrap_or(false)) { + let cells = rows.clone().flat_map(|row| cols.clone().map(move |col| (row, col))).enumerate(); + raster::collected( + w as usize, + h as usize, + cells.map(|(i, (row, col))| match mask.map(|m| m.get(i).copied().unwrap_or(false)) { Some(false) => fill, _ => b.data.get((row * r.width + col) as usize).copied().unwrap_or(fill), - }) - .collect() + }), + ) }) - .collect(); + .collect::, _>>() + .map_err(|e| e.to_string())?; let o = r.transform.to_world(f64::from(cols.start), f64::from(rows.start)); let gt = GeoTransform { origin_x: o.x, origin_y: o.y, ..r.transform }; let mut out = Raster::new(name, w, h, bands, gt, r.crs, nodata, r.color.clone(), r.data_type.clone()).map_err(|e| e.to_string())?; @@ -1125,7 +1146,7 @@ fn edit_pixels(s: &mut Session, p: &Value) -> CmdResult { - 1, ), }; - let mut bands: Vec> = r.bands.iter().map(|b| b.data.clone()).collect(); + let mut bands = copied(&r)?; let mut n = 0usize; // Only the pixels under the polygon's box are visited. let (cols, rows) = pixel_window(&r, &e).ok_or("the polygon is outside the raster")?; @@ -1340,7 +1361,7 @@ fn mosaic(s: &mut Session, p: &Value) -> CmdResult { return Err(format!("the mosaic would be {cols} × {rows} cells; set a larger cell size in the workspace options")); } let nb = imgs.iter().map(|r| r.band_count()).min().unwrap_or(1).min(if w.kind == "reality" { 1 } else { 4 }); - let mut bands = vec![vec![NODATA; cols * rows]; nb]; + let mut bands = (0..nb).map(|_| raster::filled(cols, rows, NODATA)).collect::, _>>().map_err(|e| e.to_string())?; for row in 0..rows { for col in 0..cols { let c = Coord::new(e.min_x + (col as f64 + 0.5) * cell, e.max_y - (row as f64 + 0.5) * cell); diff --git a/crates/engine/src/cmd/zarr.rs b/crates/engine/src/cmd/zarr.rs index a1dc508bf..d8add1968 100644 --- a/crates/engine/src/cmd/zarr.rs +++ b/crates/engine/src/cmd/zarr.rs @@ -935,7 +935,7 @@ impl ReduceRun { use std::sync::atomic::Ordering; // The first slice's grid (its size, place and coordinate system), all the run keeps of it. let mut first: Option<(u32, u32, raster::GeoTransform, Crs)> = None; - let mut totals = super::cube::Totals::new(self.stat, 0); + let mut totals = super::cube::Totals::new(self.stat, 0)?; // Steps never written (no chunk at all), read as having no value. let mut unwritten = 0usize; for &i in &self.steps { @@ -965,7 +965,7 @@ impl ReduceRun { if len > MAX_REDUCE_CELLS { return Err(format!("{len} cells, more than a summary takes ({MAX_REDUCE_CELLS})")); } - totals = super::cube::Totals::new(self.stat, len); + totals = super::cube::Totals::new(self.stat, len)?; first = Some((r.width, r.height, r.transform, r.crs)); } // One grid for every step: the same size and coordinate system, its place within a