feat(engine): raster commands the machine can't hold are errors, not crashes (#1291 step 1, part 3)

crates/data/src/raster/mod.rs

@@ -43,8 +43,20 @@ pub fn filled<T: Clone>(w: usize, h: usize, fill: T) -> Result<Vec<T>> {
     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<T>(w: usize, h: usize, values: impl IntoIterator<Item = T>) -> Result<Vec<T>> {
+    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<T>(n: usize, w: usize, h: usize) -> Result<Vec<T>> {
+pub fn room<T>(n: usize, w: usize, h: usize) -> Result<Vec<T>> {
     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::<f64>(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::<u16>(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.

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<Vec<f32>, 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<Totals, String> {
+        fn each<T: Clone>(cells: usize, v: T) -> Result<Vec<T>, 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<Totals>], finished: &[Option<(Columns, Vec<f32>
 /// 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<usize>) {
+fn stack_columns(outs: &mut Columns, mask: &[f32], columns: Columns, piece_mask: &[f32], whole: Option<usize>) -> Result<(), String> {
     let order = |a: &Option<f32>, b: &Option<f32>| 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<T>(v: &mut Vec<T>, 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<f32>, 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);
             }

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<Vec<Vec<f32>>, 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::<Result<_, rustgis_data::DataError>>()
+        .map_err(|e| e.to_string())
+}
+
 fn rebuild(r: &Raster, transform: GeoTransform, crs: Crs, bands: Vec<Vec<f32>>) -> Result<Raster, String> {
     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<f32> = 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<f32> = 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<Vec<f32>> = 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::<Result<Vec<_>, _>>()
+                .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<f32>> = 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::<Result<Vec<_>, _>>()
+        .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<Vec<f32>> = 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::<Result<Vec<_>, _>>()
+        .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<Vec<f32>> = 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::<Result<Vec<_>, _>>().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);

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