diff --git a/crates/engine/src/cmd/cube.rs b/crates/engine/src/cmd/cube.rs index 7895df4bf..50dab03d7 100644 --- a/crates/engine/src/cmd/cube.rs +++ b/crates/engine/src/cmd/cube.rs @@ -27,7 +27,7 @@ pub fn commands() -> Vec { r#"{"layers"?: [raster layers, a date each] | "layer"?: a Zarr layer (or "path" | "url", "variable", "index", "crs" as for zarr.reduce) or a raster of several bands, a date each (default: the active layer) | "collection": a STAC collection (desktop; read in place over an area, in pieces when large; with url, asset, bbox, maxCloud, cellSize and resampling below, and "start", "end"), "url"?: the collection's STAC API (as for stac.addCube), "asset"?: the collection's band, such as NDVI, "bbox"?: [w, s, e, n] | "view" (the collection's area), "maxCloud"?: a collection's cloud cover cap, "cellSize"?: a collection's cell size in the files' units (the overview level nearest it is read), "resampling"?: "nearest" (default: that overview's pixels as stored) | "average" (each cell the mean of the full images' pixels in it, which reads them), "band"?: 1-based (a stack's band), "dates"?: [dates] (a stack's or the bands' dates, in their order; default: read from the layers' names or files, or the bands' names), "start"?, "end"?: dates (a whole number from 1000 to 9999 is a year; an end given as a year or a month takes all of it), "months"?: [1–12] (one season across years), "test"?: "hamedRao" (default) | "yueWang" | "seasonal" | "original", "seasons"?: parts of a year for the seasonal test (default 12), "alpha"?: 0.05 (above 0, below 0.5), "minValid"?: dates a cell needs (default: half of them, at least 4; the seasonal test also needs 4 in seasons of two or more), "per"?: "year" (default) | "decade", "name"?: text, "ran"?: (internal: a trend the app has worked out)} → four raster layers: Sen's slope, the test's Z and p-value, and the slope where p < alpha"#, ), C::new("cube.reduce", "Summarise a Data Cube over Time", "Analyze › Image Analysis", reduce_cmd).params( - r#"{"layers"?: [raster layers, a date each] | "layer"?: a Zarr layer or a raster of several bands, a date each (default: the active layer) | "path"?: a Zarr store's folder | "url"?: a Zarr store's address, or a STAC collection's API (as for stac.addCube) | "collection": a STAC collection (desktop; read in place over an area, in pieces when large; with url, asset, bbox, maxCloud, cellSize and resampling below, and "start", "end"), "variable"?, "index"?, "crs"?: a Zarr variable's (as for zarr.reduce), "asset"?: the collection's band, such as NDVI, "bbox"?: [w, s, e, n] | "view" (the collection's area), "maxCloud"?: a collection's cloud cover cap, "cellSize"?: a collection's cell size in the files' units (the overview level nearest it is read), "resampling"?: "nearest" (default) | "average" (as for cube.trend), "band"?: 1-based (a stack's band), "dates"?: [dates] (a stack's or the bands' dates, in their order; default: read from the layers' names or files, or the bands' names), "stat": "mean"|"min"|"max"|"sum"|"std" (the population standard deviation, n)|"count"|"countAbove", "above"?: number (countAbove's threshold), "by"?: "all" (default: one layer) | "year" (a layer for each calendar year) | "month" (a layer for each month of the year, across years), "start"?, "end"?: dates (as for cube.trend), "months"?: [1–12] (a season across years), "minValid"?: 0–1 (a cell with a smaller share of its group's dates is left empty; default 0), "name"?: text} → raster layers of the statistic, one for each group"#, + r#"{"layers"?: [raster layers, a date each] | "layer"?: a Zarr layer or a raster of several bands, a date each (default: the active layer) | "path"?: a Zarr store's folder | "url"?: a Zarr store's address, or a STAC collection's API (as for stac.addCube) | "collection": a STAC collection (desktop; read in place over an area, in pieces when large; with url, asset, bbox, maxCloud, cellSize and resampling below, and "start", "end"), "variable"?, "index"?, "crs"?: a Zarr variable's (as for zarr.reduce), "asset"?: the collection's band, such as NDVI, "bbox"?: [w, s, e, n] | "view" (the collection's area), "maxCloud"?: a collection's cloud cover cap, "cellSize"?: a collection's cell size in the files' units (the overview level nearest it is read), "resampling"?: "nearest" (default) | "average" (as for cube.trend), "band"?: 1-based (a stack's band), "dates"?: [dates] (a stack's or the bands' dates, in their order; default: read from the layers' names or files, or the bands' names), "stat": "mean"|"min"|"max"|"sum"|"std" (the population standard deviation, n)|"count"|"countAbove"|"mode" or "majority" (a categorical cube's most frequent class; the least of those tied)|"modeCount" (how many dates had the mode)|"variety" (how many classes)|"changes" (how many times the class changed from one valued date to the next)|"longestRun" (the class of the longest run of dates in one class; the least of those tied)|"longestRunLength" (how many dates that run lasted)|"frequency" (a layer for each class met: its count of dates)|"countIn" (how many dates in the classes given), "above"?: number (countAbove's threshold), "classes"?: [whole numbers 0–255, at most 256] (countIn's classes), "by"?: "all" (default: one layer) | "year" (a layer for each calendar year) | "month" (a layer for each month of the year, across years), "start"?, "end"?: dates (as for cube.trend), "months"?: [1–12] (a season across years), "minValid"?: 0–1 (a cell with a smaller share of its group's dates is left empty; default 0), "name"?: text} → raster layers of the statistic, one for each group (frequency: one for each class of each group)"#, )] } @@ -652,7 +652,7 @@ impl TrendRun { done.fetch_add(1, Ordering::Relaxed); Ok(()) }; - read_slices(self.cube.as_ref(), &self.steps, piece, cancel, &mut place)?; + read_slices(self.cube.as_ref(), &self.steps, piece, cancel, false, &mut place)?; let grid = grid.ok_or("no date to read")?; Ok((grid, block, cells)) } @@ -743,10 +743,13 @@ fn read_slices( steps: &[usize], piece: usize, cancel: &AtomicBool, + in_order: bool, place: &mut dyn FnMut(usize, Option<(Arc, usize)>) -> Result<(), String>, ) -> Result<(), String> { - // One at a time in the browser, which has no threads to read on. - let readers = if cfg!(target_arch = "wasm32") { 1 } else { cube.readers().clamp(1, 16).min(steps.len().max(1)) }; + // One at a time in the browser, which has no threads to read on, and where the slices must come + // in their order (`in_order`): several readers place them as they are read, and one slow read + // lets the others run ahead of it. + let readers = if cfg!(target_arch = "wasm32") || in_order { 1 } else { cube.readers().clamp(1, 16).min(steps.len().max(1)) }; if readers == 1 { for (k, &i) in steps.iter().enumerate() { if cancel.load(Ordering::Relaxed) { @@ -964,6 +967,23 @@ pub enum Stat { Count, /// How many of the values are above the threshold. CountAbove(f64), + /// A categorical cube's most frequent class (the least of those tied). + Mode, + /// How many classes a cell had. + Variety, + /// How many times a cell's class changed from one valued date to the next, in date order. + Changes, + /// How many of the values are one of these classes, 0 to 255 (a bit each). + CountIn([u64; 4]), + /// How many dates had the most frequent class (the mode's count). + ModeCount, + /// The class of the longest run of dates with one class, in date order (the least of those + /// tied). + LongestRun, + /// How many dates the longest run lasted. + LongestRunLength, + /// How many dates had each class: an output for each class met. + Frequency, } impl Stat { @@ -976,12 +996,52 @@ impl Stat { Stat::Std => "std", Stat::Count => "count", Stat::CountAbove(_) => "countAbove", + Stat::Mode => "mode", + Stat::Variety => "variety", + Stat::Changes => "changes", + Stat::CountIn(_) => "countIn", + Stat::ModeCount => "modeCount", + Stat::LongestRun => "longestRun", + Stat::LongestRunLength => "longestRunLength", + Stat::Frequency => "frequency", } } - /// The statistic asked for: `stat`, and countAbove's threshold, `above`. + /// Whether it takes the values as classes (whole numbers), as of a land-cover map. + pub(crate) fn categorical(self) -> bool { + matches!( + self, + Stat::Mode + | Stat::Variety + | Stat::Changes + | Stat::CountIn(_) + | Stat::ModeCount + | Stat::LongestRun + | Stat::LongestRunLength + | Stat::Frequency + ) + } + + /// Whether it compares each date with the one before, so its dates must come in date order. + pub(crate) fn in_order(self) -> bool { + matches!(self, Stat::Changes | Stat::LongestRun | Stat::LongestRunLength) + } + + /// Whether it follows each cell's runs of dates in one class. + fn runs(self) -> bool { + matches!(self, Stat::LongestRun | Stat::LongestRunLength) + } + + /// Whether it counts each class met in each cell. + fn counts_classes(self) -> bool { + matches!(self, Stat::Mode | Stat::Variety | Stat::ModeCount | Stat::Frequency) + } + + /// The statistic asked for: `stat`, countAbove's threshold (`above`) and countIn's classes + /// (`classes`). pub(crate) fn parse(p: &Value) -> Result { - Ok(match ps(p, "stat").map(str::trim).ok_or("stat is required (mean, min, max, sum, std, count or countAbove)")? { + const STATS: &str = "mean, min, max, sum, std, count, countAbove, mode (or majority), modeCount, variety, changes, longestRun, longestRunLength, frequency or countIn"; + Ok(match ps(p, "stat").map(str::trim).ok_or_else(|| format!("stat is required ({STATS})"))? { x if x.eq_ignore_ascii_case("mean") => Stat::Mean, x if x.eq_ignore_ascii_case("min") || x.eq_ignore_ascii_case("minimum") => Stat::Min, x if x.eq_ignore_ascii_case("max") || x.eq_ignore_ascii_case("maximum") => Stat::Max, @@ -991,11 +1051,39 @@ impl Stat { x if x.eq_ignore_ascii_case("countAbove") => { Stat::CountAbove(super::f(p, "above").ok_or("countAbove takes \"above\", the threshold")?) } - x => return Err(format!("unknown stat '{x}' (mean, min, max, sum, std, count or countAbove)")), + x if x.eq_ignore_ascii_case("mode") || x.eq_ignore_ascii_case("majority") => Stat::Mode, + x if x.eq_ignore_ascii_case("variety") => Stat::Variety, + x if x.eq_ignore_ascii_case("changes") => Stat::Changes, + x if x.eq_ignore_ascii_case("countIn") => Stat::CountIn(classes(p)?), + x if x.eq_ignore_ascii_case("modeCount") => Stat::ModeCount, + x if x.eq_ignore_ascii_case("longestRun") => Stat::LongestRun, + x if x.eq_ignore_ascii_case("longestRunLength") => Stat::LongestRunLength, + x if x.eq_ignore_ascii_case("frequency") => Stat::Frequency, + x => return Err(format!("unknown stat '{x}' ({STATS})")), }) } } +/// countIn's classes (`classes`): whole numbers from 0 to 255, as a land-cover map's codes are. +fn classes(p: &Value) -> Result<[u64; 4], String> { + const NEED: &str = "countIn takes \"classes\", the classes to count: [whole numbers from 0 to 255]"; + let list = p.get("classes").and_then(Value::as_array).filter(|l| !l.is_empty()).ok_or(NEED)?; + let mut set = [0u64; 4]; + if list.len() > 256 { + return Err(NEED.into()); + } + for v in list { + let c = v.as_f64().filter(|c| c.fract() == 0.0 && (0.0..=255.0).contains(c)).ok_or(NEED)? as usize; + if let Some(word) = set.get_mut(c / 64) { + *word |= 1 << (c % 64); + } + } + Ok(set) +} + +/// The most classes a categorical summary tells apart in a run of cells: more is no land-cover map. +pub const MAX_CLASSES: usize = 1000; + /// The share of its dates a cell needs for a summary's value (`minValid`, 0 to 1; 0 by default). pub(crate) fn min_share(p: &Value) -> Result { match p.get("minValid").filter(|v| !v.is_null()) { @@ -1005,25 +1093,55 @@ pub(crate) fn min_share(p: &Value) -> Result { } /// A summary's running totals over a run of cells, only those its statistic needs: each cell's -/// count of values, one running value (the mean, the least, the most, the sum or the count above) -/// and, for std, the squared deviations (Welford's method). 4 to 20 bytes a cell. +/// count of values, one running value (the mean, the least, the most, the sum, a count, or the +/// changes so far) and, for std, the squared deviations (Welford's method); for changes, each +/// cell's last class; for the longest run (its class or length), its last class, the run it is +/// in and the longest so far; for mode, modeCount, variety and frequency, each class's count in each cell. 4 to 20 bytes +/// a cell, and 2 more for each class those four have met. pub(crate) struct Totals { stat: Stat, n: Vec, acc: Vec, m2: Vec, + /// Changes and the longest run: each cell's class at its last valued date. + last: Vec, + /// The longest run: each cell's run of dates in its last class, and its longest run's class + /// and length so far. + run: Vec, + longest: Vec<(f32, u16)>, + /// Mode, modeCount, variety and frequency: the classes met, in order, and each one's count in + /// each cell. + classes: Vec, + counts: Vec>, + /// Why the values can't be summarised as classes, at the first that can't. + refused: Option, } impl Totals { pub(crate) fn new(stat: Stat, cells: usize) -> Totals { let acc = match stat { - Stat::Count => Vec::new(), + 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], }; let m2 = if matches!(stat, Stat::Std) { vec![0.0; cells] } else { Vec::new() }; - Totals { stat, n: vec![0; cells], acc, m2 } + 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 } + } + + /// How many class counts it holds: its cells × the classes it has met. + pub(crate) fn held(&self) -> u64 { + (self.classes.len() as u64).saturating_mul(self.n.len() as u64) + } + + /// Why the values can't be summarised as classes: one that isn't a whole number, or more + /// classes than [`MAX_CLASSES`] or than the run's cells × classes cap ([`MAX_VALUES`]). + pub(crate) fn refused(&self) -> Option<&str> { + self.refused.as_deref() } /// How many cells it totals. @@ -1031,10 +1149,24 @@ impl Totals { self.n.len() } - /// Cell `c`'s next value, a finite one. + /// Cell `c`'s next value, a finite one; a categorical statistic's in date order. pub(crate) fn add(&mut self, c: usize, v: f64) { + if self.refused.is_some() || c >= self.n.len() { + return; + } + if self.stat.categorical() && v.fract() != 0.0 { + self.refused = Some(format!("{v} isn't a class: {} takes whole numbers, as a land-cover map's", self.stat.code())); + return; + } + if self.stat.counts_classes() && !self.count_class(c, v as f32) { + return; + } let Some(k) = self.n.get_mut(c) else { return }; *k += 1; + if self.stat.runs() { + self.extend_run(c, v as f32); + return; + } let k = f64::from(*k); let (Some(a), m) = (self.acc.get_mut(c), self.m2.get_mut(c)) else { return }; match self.stat { @@ -1049,15 +1181,73 @@ impl Totals { Stat::Max => *a = a.max(v), Stat::Sum => *a += v, Stat::CountAbove(t) if v > t => *a += 1.0, - Stat::CountAbove(_) | Stat::Count => {} + Stat::CountIn(set) if in_classes(&set, v) => *a += 1.0, + Stat::Changes => { + if let Some(last) = self.last.get_mut(c) { + if last.is_finite() && *last != v as f32 { + *a += 1.0; + } + *last = v as f32; + } + } + Stat::CountAbove(_) + | Stat::CountIn(_) + | Stat::Count + | Stat::Mode + | Stat::Variety + | Stat::ModeCount + | Stat::LongestRun + | Stat::LongestRunLength + | Stat::Frequency => {} } } + /// Cell `c`'s run of dates in one class, with class `v` next: one date longer in the same class, + /// a new run of one in another; the longest so far kept, the least class of those tied. + fn extend_run(&mut self, c: usize, v: f32) { + let (Some(last), Some(run), Some((class, length))) = (self.last.get_mut(c), self.run.get_mut(c), self.longest.get_mut(c)) else { + return; + }; + *run = if *last == v { run.saturating_add(1) } else { 1 }; + *last = v; + if *run > *length || (*run == *length && v < *class) { + (*class, *length) = (v, *run); + } + } + + /// Counts class `v` in cell `c`, a class met for the first time taking a count for every + /// cell; false, and said, past the classes' caps. + fn count_class(&mut self, c: usize, v: f32) -> bool { + let at = match self.classes.binary_search_by(|x| x.total_cmp(&v)) { + Ok(at) => at, + Err(at) => { + let cells = self.n.len(); + let held = (self.classes.len() as u64 + 1).saturating_mul(cells as u64); + let cap = class_cap(self.stat); + if self.classes.len() >= MAX_CLASSES || held > cap { + self.refused = Some(format!( + "more than {} classes in {cells} cells: {} takes a categorical raster (a land-cover map), up to {MAX_CLASSES} classes and {cap} cells × classes; take a smaller area", + self.classes.len(), + self.stat.code() + )); + return false; + } + self.classes.insert(at, v); + self.counts.insert(at, vec![0; cells]); + at + } + }; + if let Some(n) = self.counts.get_mut(at).and_then(|counts| counts.get_mut(c)) { + *n = n.saturating_add(1); + } + true + } + /// Each cell's statistic over `dates` dates: empty with fewer values than `min_share` of them, /// except that with no share asked for a count is 0 where a cell never had a value. pub(crate) fn finish(&self, min_share: f64, dates: usize) -> Vec { let need = (min_share.clamp(0.0, 1.0) * dates as f64 - 1e-9).ceil().max(1.0) as u32; - let count = matches!(self.stat, Stat::Count | Stat::CountAbove(_)); + let count = self.counts_zero(); self.n .iter() .enumerate() @@ -1069,11 +1259,70 @@ impl Totals { (match self.stat { Stat::Std => (self.m2.get(c).copied().unwrap_or(f64::NAN) / f64::from((*k).max(1))).sqrt(), Stat::Count => f64::from(*k), + Stat::Variety => self.counts.iter().filter(|n| n.get(c).is_some_and(|n| *n > 0)).count() as f64, + Stat::Mode => self.mode(c).map_or(f64::NAN, |(_, class)| f64::from(class)), + Stat::ModeCount => self.mode(c).map_or(0.0, |(n, _)| f64::from(n)), + Stat::LongestRun => self.longest.get(c).map_or(f64::NAN, |(class, _)| f64::from(*class)), + Stat::LongestRunLength => self.longest.get(c).map_or(0.0, |(_, length)| f64::from(*length)), _ => a, }) as f32 }) .collect() } + + /// Whether a cell with no value counts 0 (with no share asked for) rather than none: a count's. + fn counts_zero(&self) -> bool { + matches!( + self.stat, + Stat::Count + | Stat::CountAbove(_) + | Stat::CountIn(_) + | Stat::Changes + | Stat::Variety + | Stat::ModeCount + | Stat::LongestRunLength + | Stat::Frequency + ) + } + + /// Cell `c`'s most frequent class and its count: the least of the classes tied, as they are in + /// order. + fn mode(&self, c: usize) -> Option<(u16, f32)> { + let mut best: Option<(u16, f32)> = None; + for (class, counts) in self.classes.iter().zip(&self.counts) { + let n = counts.get(c).copied().unwrap_or(0); + if n > 0 && best.is_none_or(|(most, _)| n > most) { + best = Some((n, *class)); + } + } + best + } + + /// The outputs, as columns of cells: one for the statistic (keyed `None`), or, for frequency, + /// one for each class met, in order, each cell's count of dates in it. + pub(crate) fn columns(&self, min_share: f64, dates: usize) -> Vec<(Option, Vec)> { + if self.stat != Stat::Frequency { + return vec![(None, self.finish(min_share, dates))]; + } + let mask = self.mask(min_share, dates); + self.classes + .iter() + .zip(&self.counts) + .map(|(class, counts)| (Some(*class), mask.iter().zip(counts).map(|(m, n)| m + f32::from(*n)).collect())) + .collect() + } + + /// Each cell's value in a class's column where it never had that class: 0, or none where it + /// has fewer values than `min_share` of its dates (a count's rule). + pub(crate) fn mask(&self, min_share: f64, dates: usize) -> Vec { + let need = (min_share.clamp(0.0, 1.0) * dates as f64 - 1e-9).ceil().max(1.0) as u32; + self.n.iter().map(|k| if *k < need && min_share != 0.0 { f32::NAN } else { 0.0 }).collect() + } +} + +/// Whether class `v` is one of countIn's. +fn in_classes(set: &[u64; 4], v: f64) -> bool { + (0.0..=255.0).contains(&v) && set.get(v as usize / 64).is_some_and(|w| w >> (v as usize % 64) & 1 == 1) } /// The most values a summary's outputs hold (cells × outputs): 32 million on the desktop, 8 @@ -1081,9 +1330,13 @@ impl Totals { /// are read when it doesn't. Its running totals are a piece's at a time and only for the groups /// whose dates are being read: one year's in a summary by year read one date at a time (those of /// the dates in flight with several readers), but every month's in one by month. They take 4 bytes -/// a cell for count, 12 for mean, min, max, sum and countAbove, and 20 for std, so a summary by -/// month at the cap holds up to five times its outputs' 4 bytes a value: some 640 MB for std on the -/// desktop. A STAC cube is read in pieces whose cells × dates stay within [`MAX_VALUES`], so a long +/// a cell for count, 12 for mean, min, max, sum, countAbove and countIn, 16 for changes and the +/// longest run (its class or length) and 20 for std, so a summary by month at the cap holds up to five times its outputs' +/// 4 bytes a value: some 640 MB for std on the desktop. Mode, modeCount, variety and frequency take +/// 4 bytes a cell and 2 more for each class met, the cells × classes of the groups being read at +/// most [`MAX_VALUES`] in all (128 MB of counts; for frequency [`MAX_REDUCED`], its outputs' cap), +/// checked as each slice is placed. Frequency's outputs are a column for each class of each group, all within +/// this cap. A STAC cube is read in pieces whose cells × dates stay within [`MAX_VALUES`], so a long /// series keeps each piece's totals small. pub const MAX_REDUCED: u64 = MAX_VALUES / 2; @@ -1122,8 +1375,9 @@ impl ReduceRun { } /// Reads the slices piece by piece ([`TimeCube::pieces`]), each folded into its group's - /// running totals, and puts the pieces' outputs together row by row: a raster for each group. - /// `done` counts the slices read; `cancel` stops it between slices. + /// running totals, and puts the pieces' outputs together row by row: a raster for each group + /// (for frequency, one for each class of each group). `done` counts the slices read; `cancel` + /// stops it between slices. pub fn run(&self, cancel: &AtomicBool, done: &AtomicUsize) -> Result, String> { let outputs = self.groups.len(); let too_many = |cells: usize| { @@ -1138,7 +1392,11 @@ impl ReduceRun { return Err(e); } let mut stitch = Stitch::default(); - let mut bands: Vec> = vec![Vec::new(); outputs]; + // Each group's columns so far, by class for frequency; and for frequency each group's + // value where a class is missing, for a class met first in a later piece. + let mut bands: Vec = vec![Vec::new(); outputs]; + let mut masks: Vec> = vec![Vec::new(); outputs]; + let frequency = self.stat == Stat::Frequency; // The cells of the pieces read so far, for a cube that didn't know its size. let mut read_cells = 0usize; for piece in 0..self.cube.pieces().max(1) { @@ -1151,7 +1409,7 @@ impl ReduceRun { // order of arrival can move an output's last bit at most, and only then. let mut totals: Vec> = (0..outputs).map(|_| None).collect(); let mut seen = vec![0usize; outputs]; - let mut finished: Vec>> = vec![None; outputs]; + let mut finished: Vec)>> = vec![None; outputs]; let mut place = |k: usize, slice: Option<(Arc, usize)>| -> Result<(), String> { let g = self.group_of.get(k).copied().filter(|g| *g < outputs).ok_or("a date in no group")?; // A date with nothing written only lowers the counts. @@ -1174,6 +1432,13 @@ impl ReduceRun { t.add(c, f64::from(*v)); } } + if let Some(e) = t.refused() { + return Err(e.to_string()); + } + // The groups being read at once (every month's, by month) share the cap. + if let Some(e) = classes_over(&totals, &finished, class_cap(self.stat)) { + return Err(e); + } } let n = seen.get_mut(g).ok_or("a date in no group")?; *n += 1; @@ -1181,20 +1446,18 @@ impl ReduceRun { if *n >= dates && let (Some(t), Some(out)) = (totals.get_mut(g).and_then(Option::take), finished.get_mut(g)) { - *out = Some(t.finish(self.min_share, dates)); + *out = Some(finish(&t, self.min_share, dates, frequency)); + if let Some(e) = classes_over(&totals, &finished, class_cap(self.stat)) { + return Err(e); + } } done.fetch_add(1, Ordering::Relaxed); Ok(()) }; - read_slices(self.cube.as_ref(), &self.steps, piece, cancel, &mut place)?; + // Changes and the longest run compare each date with the one before: read in date order. + read_slices(self.cube.as_ref(), &self.steps, piece, cancel, self.stat.in_order(), &mut place)?; let grid = grid.ok_or("none of the dates in the period has data written")?; stitch.check(&grid)?; - // Their full length when the cube knows it, so they don't grow by doubling. - if let (0, Some(total)) = (piece, whole) { - for band in &mut bands { - band.reserve_exact(total); - } - } let cells = (grid.width as usize).saturating_mul(grid.height as usize); read_cells = read_cells.saturating_add(cells); if whole.is_none() @@ -1202,15 +1465,30 @@ impl ReduceRun { { return Err(e); } - for (g, band) in bands.iter_mut().enumerate() { + for (g, (band, mask)) in bands.iter_mut().zip(&mut masks).enumerate() { let dates = self.groups.get(g).map_or(0, |(_, d)| *d); // A group with no data written in this piece: no values, as its totals would give. - let out = match (finished.get_mut(g).and_then(Option::take), totals.get_mut(g).and_then(Option::take)) { + 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)) => t.finish(self.min_share, dates), - (None, None) => Totals::new(self.stat, cells).finish(self.min_share, dates), + (None, Some(t)) => finish(&t, self.min_share, dates, frequency), + (None, None) => finish(&Totals::new(self.stat, cells), self.min_share, dates, frequency), }; - band.extend(out); + // Their full length when the cube knows it, so they don't grow by doubling. + stack_columns(band, mask, columns, &piece_mask, whole); + if let (0, Some(total)) = (piece, whole) + && frequency + { + mask.reserve_exact(total); + } + mask.extend(piece_mask); + } + // Frequency's columns are known as they are met: their values within the outputs' cap. + let columns: usize = bands.iter().map(Vec::len).sum(); + let values = (read_cells as u64).saturating_mul(columns as u64); + if values > MAX_REDUCED { + return Err(format!( + "{read_cells} cells × {columns} classes and groups is {values} values, more than {MAX_REDUCED}: take a smaller area, a coarser cell size or fewer groups" + )); } stitch.add(grid); } @@ -1218,28 +1496,92 @@ impl ReduceRun { if let Some(e) = self.cube.nothing_read() { return Err(e); } + if bands.iter().all(Vec::is_empty) { + return Err("none of the dates in the period has a class".into()); + } done.store(self.total(), Ordering::Relaxed); + let grid = &grid; self.groups .iter() .zip(bands) - .map(|((name, _), band)| { - Raster::new( - name.clone(), - grid.width, - height, - vec![band], - grid.transform, - grid.crs, - None, - raster::ColorHint::Gray, - "Float32", - ) - .map_err(|e| e.to_string()) + .flat_map(|((name, _), columns)| { + columns.into_iter().map(move |(class, band)| { + let name = match class { + Some(class) => format!("{name} class {class}"), + None => name.clone(), + }; + Raster::new(name, grid.width, height, vec![band], grid.transform, grid.crs, None, raster::ColorHint::Gray, "Float32") + .map_err(|e| e.to_string()) + }) }) .collect() } } +/// A group's outputs, as columns of cells: one keyed `None`, or frequency's one for each class. +type Columns = Vec<(Option, Vec)>; + +/// A group's columns from its totals, and for frequency the value of a class it didn't meet. +fn finish(t: &Totals, min_share: f64, dates: usize, frequency: bool) -> (Columns, Vec) { + let mask = if frequency { t.mask(min_share, dates) } else { Vec::new() }; + (t.columns(min_share, dates), mask) +} + +/// The most cells × classes the groups being read count at once: [`MAX_VALUES`], or for frequency, +/// whose classes are each an output column, the outputs' cap [`MAX_REDUCED`], so it is refused +/// before its counts and columns are made rather than after. +fn class_cap(stat: Stat) -> u64 { + if stat == Stat::Frequency { MAX_REDUCED } else { MAX_VALUES } +} + +/// Why the class counts of the groups being read (`totals`) are too many: more than `cap` cells × +/// classes in all, including columns retained by completed groups. +fn classes_over(totals: &[Option], finished: &[Option<(Columns, Vec)>], cap: u64) -> Option { + let completed = finished + .iter() + .flatten() + .flat_map(|(columns, _)| columns.iter()) + .map(|(_, values)| values.len() as u64) + .fold(0u64, u64::saturating_add); + let held = totals.iter().flatten().map(Totals::held).fold(completed, u64::saturating_add); + (held > cap).then(|| { + format!( + "{held} cells × classes across the groups being read, more than {cap}: take a smaller area, a coarser cell size, or a group by year rather than by month" + ) + }) +} + +/// 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) { + let order = |a: &Option, b: &Option| match (a, b) { + (Some(a), Some(b)) => a.total_cmp(b), + _ => a.is_some().cmp(&b.is_some()), + }; + let before = outs.first().map_or(mask.len(), |(_, b)| b.len()); + for (key, column) in columns { + match outs.binary_search_by(|(k, _)| order(k, &key)) { + Ok(at) => { + if let Some((_, b)) = outs.get_mut(at) { + b.extend(column); + } + } + Err(at) => { + let mut b = Vec::with_capacity(whole.unwrap_or(0).max(mask.len() + column.len())); + b.extend_from_slice(mask); + b.extend(column); + outs.insert(at, (key, b)); + } + } + } + for (_, b) in outs.iter_mut() { + if b.len() == before { + b.extend_from_slice(piece_mask); + } + } +} + /// Prepares `cube.reduce`: the cube, the statistic, and the dates of the period (and months) /// asked for, in their groups. pub fn reduce_start(s: &Session, p: &Value) -> Result { @@ -1253,8 +1595,13 @@ pub fn reduce_start(s: &Session, p: &Value) -> Result { Some(x) => return Err(format!("by is all, year or month, not {x}")), }; let cube = cube_from(s, p, SUMMARY)?; - let steps = steps(cube.dates(), cube.period(), p)?; + let mut steps = steps(cube.dates(), cube.period(), p)?; let date = |i: usize| cube.dates().get(i).copied().unwrap_or(f64::NAN); + // Changes and runs are counted from each date to the next: the dates in order, which a Zarr + // axis out of order doesn't give. + if stat.in_order() { + steps.sort_by(|x, y| date(*x).total_cmp(&date(*y))); + } // A date that can't be read would go to no year or month: said, not put in a group of its own. if let Some(&i) = steps.iter().find(|&&i| !date(i).is_finite()) { return Err(format!("the cube's date {} can't be read: give the dates (\"dates\")", i + 1)); @@ -1346,4 +1693,85 @@ mod tests { assert!(super::same_grid(&a, &hair)); assert!(!super::same_grid(&a, &half)); } + + /// Several readers place a cube's slices as they are read; in order, one reads them, in the + /// steps' order, however long each read takes. + #[cfg(not(target_arch = "wasm32"))] + #[test] + fn slices_asked_in_order_come_in_order() { + use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster}; + use std::sync::Arc; + use std::sync::atomic::AtomicBool; + struct Slow(Vec); + impl super::TimeCube for Slow { + fn dates(&self) -> &[f64] { + &self.0 + } + fn slice(&self, i: usize) -> Result, usize)>, String> { + // The early dates the slowest, so readers that place as they read finish them last. + std::thread::sleep(std::time::Duration::from_millis(((12 - i.min(12)) * 3) as u64)); + let gt = GeoTransform { origin_x: 0.0, origin_y: 0.0, pixel_width: 1.0, pixel_height: -1.0 }; + let r = Raster::new("s", 1, 1, vec![vec![i as f32]], gt, Crs::LonLat, None, ColorHint::Gray, "Float32") + .map_err(|e| e.to_string())?; + Ok(Some((Arc::new(r), 0))) + } + fn label(&self) -> String { + "slow".into() + } + fn readers(&self) -> usize { + 8 + } + } + let cube = Slow((0..12).map(f64::from).collect()); + let steps: Vec = (0..12).collect(); + let mut came = Vec::new(); + let read = super::read_slices(&cube, &steps, 0, &AtomicBool::new(false), true, &mut |k, _| { + came.push(k); + Ok(()) + }); + assert!(read.is_ok(), "{read:?}"); + assert_eq!(came, steps); + } + + /// Frequency's columns over pieces: a class met first in a later piece gets the earlier + /// pieces' mask (0, or none where a cell had too few values), and a class a piece didn't meet + /// gets that piece's. + #[test] + fn a_class_met_in_a_later_piece_is_counted_0_before() { + 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); + mask.extend([0.0, nan]); + super::stack_columns(&mut outs, &mask, vec![(Some(5.0), vec![4.0])], &[0.0], Some(3)); + 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)); + 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)); + assert_eq!(one, vec![(None, vec![1.0, 2.0, 3.0])]); + } + + /// The class counts of the groups being read share one cap: two groups each within it can be + /// over it together. + #[test] + fn groups_being_read_share_the_class_cap() { + let group = |classes: &[f64]| { + let mut t = super::Totals::new(super::Stat::Mode, 10); + for (c, v) in classes.iter().enumerate() { + t.add(c, *v); + } + t + }; + // 10 cells × 3 classes and 10 × 2: 30 and 20 counts, 50 in all. + let totals = vec![Some(group(&[1.0, 2.0, 3.0])), None, Some(group(&[1.0, 2.0]))]; + assert!(super::classes_over(&totals, &[], 50).is_none()); + let finished = vec![Some((vec![(Some(1.0), vec![0.0; 20])], Vec::new()))]; + assert!(super::classes_over(&totals, &finished, 70).is_none()); + assert!(super::classes_over(&totals, &finished, 69).is_some(), "finished groups still consume the shared cap"); + let e = super::classes_over(&totals, &[], 49).unwrap_or_default(); + assert!(e.contains("50 cells × classes across the groups being read, more than 49"), "{e}"); + } } diff --git a/crates/engine/src/cmd/zarr.rs b/crates/engine/src/cmd/zarr.rs index e0a4a4943..e451ad201 100644 --- a/crates/engine/src/cmd/zarr.rs +++ b/crates/engine/src/cmd/zarr.rs @@ -975,6 +975,13 @@ impl ReduceRun { #[cfg(not(target_arch = "wasm32"))] pub fn reduce_start(s: &Session, p: &Value) -> Result { let stat = Stat::parse(p)?; + // A variable of classes is summarised by cube.reduce, which reads its dates in order. + if stat.categorical() { + return Err(format!( + "zarr.reduce takes mean, min, max, sum, std, count or countAbove: {} is cube.reduce's (give it this layer or store)", + stat.code() + )); + } let min_valid = super::cube::min_share(p)?; let ZarrCube { zt, period } = zarr_cube(s, p)?; let steps = super::cube::steps(&zt.dates, period, p)?; diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index a3cea6946..dbdc028f3 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -1316,6 +1316,11 @@ fn hostile_params_never_panic() { // A summary over time: a statistic given, so its threshold, groups and share are reached. json!({"stat": "countAbove", "above": "x", "by": 7, "minValid": -1e308, "months": [13, "x", null], "name": {"a": 1}}), json!({"stat": "std", "by": "month", "minValid": f64::NAN, "dates": ["x", 7, null]}), + // A categorical summary: countIn's classes, and a mode by groups with a share out of range. + json!({"stat": "countIn", "classes": [1e308, -1, "x", null, 3.5, {"a": 1}]}), + json!({"stat": "countIn", "classes": "all", "by": "year"}), + json!({"stat": "mode", "by": "month", "minValid": 2, "months": []}), + json!({"stat": "frequency", "by": "year", "minValid": -1, "name": []}), // MySQL addresses that can't be read, so nothing connects. json!({"database": "mysql://u:p#@[::1/x?ssl-mode=nope", "password": 7, "layer": ":", "crs": "EPSG:-1"}), json!({"database": "MARIADB://u@h:99999999/%zz?ssl-mode=VERIFY_CA&ssl-ca=/nonexistent", "sql": "SELECT 1"}), @@ -9327,6 +9332,169 @@ fn a_short_cube_is_summarised_but_not_trended() { assert!(e.contains("summary's stack needs 2 layers"), "{e}"); } +/// A stack of land-cover maps summarised as classes: each cell's most frequent class (the least +/// of those tied) and how many dates had it, how many classes it had, how many times its class +/// changed and the class and length of its longest run, in date order (the layers added out of it), how many +/// dates it was in some classes, and a layer for each class with its count of dates. A cell left +/// without a class at a date (no data) is skipped. What isn't a class, too many classes, and a categorical +/// statistic asked of zarr.reduce are refused. +#[cfg(not(target_arch = "wasm32"))] +#[test] +fn a_categorical_stack_summarised_over_time() { + use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster, geotiff}; + let gt = GeoTransform { origin_x: -44.0, origin_y: -7.0, pixel_width: 0.01, pixel_height: -0.01 }; + let add = |s: &mut Session, name: &str, cells: Vec| { + let r = Raster::new(name, cells.len() as u32, 1, vec![cells], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap(); + let bytes = geotiff::write(&r, 0).unwrap(); + crate::cmd::raster::add_raster(s, name, rustgis_map::RasterLayer::from_bytes(r, format!("{name}.tif"), bytes)).unwrap(); + }; + let nan = f32::NAN; + // Four cells over six years: 3 3 4 3 4 4; always 15; 21 21 33 with no data between; and + // 5 7 7 7 5 5, whose mode (5, the least of two tied) isn't its longest run's class (7). + let years: [(u32, [f32; 4]); 6] = [ + (2000, [3.0, 15.0, nan, 5.0]), + (2001, [3.0, 15.0, 21.0, 7.0]), + (2002, [4.0, 15.0, nan, 7.0]), + (2003, [3.0, 15.0, 21.0, 7.0]), + (2004, [4.0, 15.0, 33.0, 5.0]), + (2005, [4.0, 15.0, nan, 5.0]), + ]; + let mut s = Session::new(Project::default()); + let mut names = Vec::new(); + // Added out of date order: changes go by the dates, not the layers' order. + for k in [3, 0, 5, 1, 4, 2] { + let (y, cells) = years[k]; + let name = format!("cobertura_{y}"); + add(&mut s, &name, cells.to_vec()); + names.push(name); + } + let run = |s: &mut Session, p: Value| { + let r = s.run("cube.reduce", &p).unwrap(); + let id = r["layers"][0].as_u64().unwrap() as u32; + let l = s.map().layer(id).unwrap(); + (l.name.clone(), l.raster.as_ref().unwrap().data.as_ref().unwrap().bands[0].data.clone()) + }; + let (name, mode) = run(&mut s, json!({"layers": names, "stat": "mode"})); + assert_eq!(mode, vec![3.0, 15.0, 21.0, 5.0]); + assert_eq!(name, "cobertura mode 2000-01-01–2005-01-01"); + assert_eq!(run(&mut s, json!({"layers": names, "stat": "majority"})).1, mode); + assert_eq!(run(&mut s, json!({"layers": names, "stat": "modeCount"})).1, vec![3.0, 6.0, 2.0, 3.0]); + assert_eq!(run(&mut s, json!({"layers": names, "stat": "variety"})).1, vec![2.0, 1.0, 2.0, 2.0]); + assert_eq!(run(&mut s, json!({"layers": names, "stat": "changes"})).1, vec![3.0, 0.0, 1.0, 2.0]); + // Runs of 2 in 3 and in 4 (the least), 6 in 15, 21 through its gap, 3 in 7. + assert_eq!(run(&mut s, json!({"layers": names, "stat": "longestRun"})).1, vec![3.0, 15.0, 21.0, 7.0]); + assert_eq!(run(&mut s, json!({"layers": names, "stat": "longestRunLength"})).1, vec![2.0, 6.0, 2.0, 3.0]); + assert_eq!(run(&mut s, json!({"layers": names, "stat": "countIn", "classes": [4, 33]})).1, vec![3.0, 0.0, 1.0, 0.0]); + // A share of the dates: the third cell has a class at half of them, so it's left out at 0.6. + let (_, half) = run(&mut s, json!({"layers": names, "stat": "mode", "minValid": 0.6})); + assert_eq!((half[0], half[1], half[3]), (3.0, 15.0, 5.0)); + assert!(half[2].is_nan(), "{half:?}"); + // Frequency: a layer for each class met, each cell's count of dates in it. + let layers = |s: &Session, r: &Value| -> Vec<(String, Vec)> { + r["layers"] + .as_array() + .unwrap() + .iter() + .map(|id| { + let l = s.map().layer(id.as_u64().unwrap() as u32).unwrap(); + (l.name.clone(), l.raster.as_ref().unwrap().data.as_ref().unwrap().bands[0].data.clone()) + }) + .collect() + }; + let r = s.run("cube.reduce", &json!({"layers": names, "stat": "frequency"})).unwrap(); + let got = layers(&s, &r); + let want: Vec<(String, Vec)> = [ + (3, [3.0, 0.0, 0.0, 0.0]), + (4, [3.0, 0.0, 0.0, 0.0]), + (5, [0.0, 0.0, 0.0, 3.0]), + (7, [0.0, 0.0, 0.0, 3.0]), + (15, [0.0, 6.0, 0.0, 0.0]), + (21, [0.0, 0.0, 2.0, 0.0]), + (33, [0.0, 0.0, 1.0, 0.0]), + ] + .iter() + .map(|(c, n)| (format!("cobertura frequency 2000-01-01–2005-01-01 class {c}"), n.to_vec())) + .collect(); + assert_eq!(got, want); + // With a share, a cell with too few values is left out of every class's layer. + let r = s.run("cube.reduce", &json!({"layers": names, "stat": "frequency", "minValid": 0.6, "name": "freq"})).unwrap(); + let got = layers(&s, &r); + assert_eq!(got.len(), 7); + assert!(got.iter().all(|(_, n)| n[2].is_nan()), "{got:?}"); + assert_eq!(got[0].0, "freq class 3"); + assert_eq!((got[0].1[0], got[0].1[1], got[0].1[3]), (3.0, 0.0, 0.0)); + // A raster whose bands are the years out of order, their dates given in the bands' order: + // changes still go by the dates. + let order = [3, 0, 5, 1, 4, 2]; + let bands: Vec> = order.iter().map(|&k| years[k].1.to_vec()).collect(); + let dates: Vec = order.iter().map(|&k| format!("{}-01-01", years[k].0)).collect(); + let r = Raster::new("bandas", 4, 1, bands, gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap(); + let bytes = geotiff::write(&r, 0).unwrap(); + crate::cmd::raster::add_raster(&mut s, "bandas", rustgis_map::RasterLayer::from_bytes(r, "bandas.tif".to_string(), bytes)).unwrap(); + assert_eq!(run(&mut s, json!({"layer": "bandas", "dates": dates, "stat": "changes"})).1, vec![3.0, 0.0, 1.0, 2.0]); + assert_eq!(run(&mut s, json!({"layer": "bandas", "dates": dates, "stat": "longestRun"})).1, vec![3.0, 15.0, 21.0, 7.0]); + assert_eq!(run(&mut s, json!({"layer": "bandas", "dates": dates, "stat": "longestRunLength"})).1, vec![2.0, 6.0, 2.0, 3.0]); + // By year, a change is one within the year: none, with one date a year. + let r = s.run("cube.reduce", &json!({"layers": names, "stat": "changes", "by": "year"})).unwrap(); + assert_eq!(r["layers"].as_array().unwrap().len(), 6); + // Every key a categorical summary takes is read, none said to be ignored. + let mark = s.notice_mark(); + run(&mut s, json!({"layers": names, "stat": "countIn", "classes": [3], "minValid": 0.5, "name": "floresta"})); + let ignored: Vec = s.notices_since(mark).iter().map(|n| n.text.clone()).filter(|t| t.contains("ignores")).collect(); + assert!(ignored.is_empty(), "{ignored:?}"); + // Refused, nothing added. + let before = s.map().layers.len(); + for (bad, said) in [ + (json!({"layers": names, "stat": "countIn"}), "countIn takes \"classes\""), + (json!({"layers": names, "stat": "countIn", "classes": []}), "countIn takes \"classes\""), + (json!({"layers": names, "stat": "countIn", "classes": [3.5]}), "whole numbers from 0 to 255"), + (json!({"layers": names, "stat": "countIn", "classes": [256]}), "whole numbers from 0 to 255"), + (json!({"layers": names, "stat": "countIn", "classes": "3"}), "countIn takes \"classes\""), + (json!({"layers": names, "stat": "countIn", "classes": vec![3; 257]}), "countIn takes \"classes\""), + ] { + let e = s.run("cube.reduce", &bad).unwrap_err(); + assert!(e.contains(said), "{bad}: {e}"); + } + // Values that aren't classes (an NDVI): said, with the value. + let ndvi = ["ndvi_2020", "ndvi_2021"]; + for (k, name) in ndvi.iter().enumerate() { + add(&mut s, name, vec![0.25 + k as f32 * 0.5]); + } + let before_ndvi = s.map().layers.len(); + for stat in ["mode", "variety", "changes", "countIn"] { + let e = s.run("cube.reduce", &json!({"layers": ndvi, "stat": stat, "classes": [0]})).unwrap_err(); + assert!(e.contains("0.25 isn't a class") && e.contains(stat), "{stat}: {e}"); + } + assert_eq!(s.map().layers.len(), before_ndvi); + // More classes than a land-cover map has: refused, not a count for each. + let many = ["muitas_2020", "muitas_2021"]; + for name in many { + add(&mut s, name, (0..1001).map(|c| c as f32).collect()); + } + let after_many = s.map().layers.len(); + let e = s.run("cube.reduce", &json!({"layers": many, "stat": "variety"})).unwrap_err(); + assert!(e.contains("more than 1000 classes") && e.contains("variety"), "{e}"); + assert_eq!(s.map().layers.len(), after_many); + assert_eq!(before_ndvi - before, ndvi.len()); + // A Zarr axis out of date order (days 31, 0, 59) whose classes are 2, 1, 2 along it: in date + // order they are 1, 2, 2, one change, not the axis's two, and a longest run in 2, where the + // axis's runs are all of one date (and the least class, 1, would win). + let (dir, path) = zarr_store_at("cube-reduce-classes-order", &[31, 0, 59]); + for (t, class) in [(0, 2.0f32), (1, 1.0), (2, 2.0)] { + let chunk: Vec = std::iter::repeat_n(class, 12).flat_map(|x| x.to_le_bytes()).collect(); + std::fs::write(std::path::Path::new(&path).join(format!("t/{t}.0.0")), chunk).unwrap(); + } + assert_eq!(run(&mut s, json!({"path": path, "variable": "t", "stat": "changes"})).1, vec![1.0; 12]); + assert_eq!(run(&mut s, json!({"path": path, "variable": "t", "stat": "longestRun"})).1, vec![2.0; 12]); + assert_eq!(run(&mut s, json!({"path": path, "variable": "t", "stat": "longestRunLength"})).1, vec![2.0; 12]); + let _ = std::fs::remove_dir_all(dir); + // zarr.reduce doesn't take them: it points to cube.reduce. + let (dir, path) = zarr_store_over_time("cube-reduce-classes"); + let e = s.run("zarr.reduce", &json!({"path": path, "variable": "t", "stat": "mode"})).unwrap_err(); + assert!(e.contains("cube.reduce"), "{e}"); + let _ = std::fs::remove_dir_all(dir); +} + /// A Zarr axis in decreasing order: a trend refuses it, a summary takes it and names its period from /// the earliest date to the latest. #[test] diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 052823219..5a2270730 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -472,14 +472,42 @@ variable, a STAC collection's band): **mean**, **min**, **max**, **sum**, **std* standard deviation), **count** of values, or **countAbove** a threshold (`"above"`), as Summarise over time does for a Zarr variable (see [Data](data.md)). +- **Classes:** a series of categorical rasters, such as yearly land-cover maps (MapBiomas, ESA + WorldCover), is summarised as classes with eight more statistics: + - **mode**, each cell's most frequent class, the least of those tied (`"majority"` too); + - **modeCount**, at how many dates it had that class; + - **variety**, how many classes it had; + - **changes**, how many times its class changed from one date with a class to the next, in date + order. A date with no data is skipped, so A, no data, B is one change. Within each group, so by + year a change is one within the year; + - **longestRun**, the class it kept for the most dates in a row, in date order, the least of those + tied. A date with no data neither breaks a run nor adds to it, so A, no data, A is a run of two; + - **longestRunLength**, how many dates that run lasted; + - **frequency**, a layer for each class met (`… frequency 1985–2023 class 3`), each cell's count + of dates in that class. By year or month, a layer for each class of each group; + - **countIn**, at how many dates it was one of the classes given (`"classes": [3, 4, 5]`, whole + numbers from 0 to 255), such as the years a cell was forest. + + They take the values as classes: a value that isn't a whole number (an NDVI) is refused with the + value, as is a run with more than 1,000 classes, or more than 64 million cells × classes in the + groups being read at once (16 million in the web app; for frequency, whose classes are each a + layer, the outputs' 32 million values). These are the data cube summary's only: Summarise over time for a + Zarr variable refuses them and points here. Mode, modeCount, variety and frequency keep a count + of each class in each cell while they read, 2 bytes a class, for each group being read: by + month, every month's at once, so a large area by month may need a smaller area or a group by + year. **changes**, + **longestRun** and **longestRunLength** read a STAC collection one date at a time, in order, so + they are slower there than the other statistics, which read several dates at once. + - **Groups:** `"by": "all"` (the default) makes one layer of the whole period, `"year"` one for each calendar year (an annual mean NDVI from 16-day composites), `"month"` one for each month of the year across the years (a climatology). `"start"`, `"end"` and `"months"` pick the dates as for a trend. - **Missing values:** no-data and missing values are skipped cell by cell. A cell keeps no value where it had fewer values than `"minValid"` (a share of its group's dates, 0 by default: 0.5 - wants half of a year's composites); with the default, **count** and **countAbove** give 0 where - a cell never had a value. + wants half of a year's composites); with the default, **count**, **countAbove**, **countIn**, + **modeCount**, **variety**, **changes**, **longestRunLength** and **frequency** give 0 where a cell + never had a value, and **mode** and **longestRun** give none. - **Size:** the slices are folded into running totals as they are read, so a long period takes no more memory than a short one; a STAC collection's area is read in pieces, as for a trend. The outputs hold up to 32 million values (cells times groups; 8 million in the web app). A summary by @@ -491,6 +519,7 @@ over time does for a Zarr variable (see [Data](data.md)). ```sh rustgis-cli run --cmd 'cube.reduce={"layers":["ndvi_2015","ndvi_2016","ndvi_2017","ndvi_2018"],"stat":"max"}' +rustgis-cli run --data cobertura_1985.tif --data cobertura_1986.tif --data cobertura_1987.tif --cmd 'cube.reduce={"layers":["cobertura_1985","cobertura_1986","cobertura_1987"],"stat":"changes"}' rustgis-cli run --cmd 'cube.reduce={"url":"https://data.inpe.br/bdc/stac/v1","collection":"mod13q1-6.1","asset":"NDVI","bbox":[-44.8,-7.5,-44.3,-7.0],"start":"2001","end":"2025","cellSize":1000,"stat":"mean","by":"year","minValid":0.5,"name":"NDVI"}' --save urucui.rgis ```