diff --git a/crates/data/src/raster/geotiff.rs b/crates/data/src/raster/geotiff.rs index 54d77b4f8..66d9569f1 100644 --- a/crates/data/src/raster/geotiff.rs +++ b/crates/data/src/raster/geotiff.rs @@ -878,6 +878,13 @@ pub(super) fn read_limited(r: Reader<'_>, name: &str, max: usize, open: Option, name: &str, max: usize, open: Option, name: &str, max: usize, open: Option(2, 2).unwrap(); + let e = img.encoder(); + e.write_tag(Tag::ModelPixelScaleTag, &[0.5f64, 0.5, 0.0][..]).unwrap(); + e.write_tag(Tag::ModelTiepointTag, &[0.0f64, 0.0, 0.0, 10.0, 50.0, 0.0][..]).unwrap(); + e.write_tag(Tag::GeoKeyDirectoryTag, &[1u16, 1, 0, 2, 1024, 0, 1, 2, 2048, 0, 1, 4326][..]).unwrap(); + let xml = r#"NDVI 2024-06-01"#; + e.write_tag(Tag::Unknown(42112), xml).unwrap(); + img.write_data(&[0.1f32, 0.2, 0.3, 0.4]).unwrap(); + } + let r = read(std::io::Cursor::new(out.into_inner()), "ndvi").unwrap(); + assert_eq!(r.bands[0].name, "NDVI 2024-06-01"); + // Without the tag, the band keeps its number. + let plain = read(std::io::Cursor::new(grey8(2, 2, &[1, 2, 3, 4], None)), "plain").unwrap(); + assert_ne!(plain.bands[0].name, "NDVI 2024-06-01"); + } + #[test] fn a_files_own_projection_is_read_and_written_back() { // The NLCD's Albers, as its GeoKeys give it. diff --git a/crates/data/src/raster/pam.rs b/crates/data/src/raster/pam.rs index a75daff98..6e6d32fdc 100644 --- a/crates/data/src/raster/pam.rs +++ b/crates/data/src/raster/pam.rs @@ -118,6 +118,42 @@ pub fn parse(text: &str) -> Option { Some(pam) } +/// The band descriptions a GeoTIFF holds itself, in GDAL's metadata tag (`GDAL_METADATA`, 42112), +/// where GDAL and rasterio write them: `None` when it isn't such a document or describes no band. +/// Its `sample` counts from 0. +/// +/// ```xml +/// +/// NDVI 2024-06-01 +/// +/// ``` +pub fn parse_gdal_metadata(text: &str) -> Option { + let opts = ParsingOptions { allow_dtd: false, nodes_limit: crate::io::XML_NODES, ..ParsingOptions::default() }; + let doc = crate::xml::parse(text, opts).ok()?; + let root = doc.root_element(); + if root.tag_name().name() != "GDALMetadata" { + return None; + } + let mut pam = Pam::default(); + // The bands described so far, so a long document isn't searched again for each item. + let mut seen = std::collections::HashSet::new(); + for item in root.children().filter(|n| n.has_tag_name("Item") && n.attribute("role") == Some("description")) { + let band = item + .attribute("sample") + .and_then(|v| v.trim().parse::().ok()) + .and_then(|s| s.checked_add(1)) + .filter(|n| *n <= MAX_BANDS); + let description = item.text().map(str::trim).filter(|d| !d.is_empty()).map(|d| d.chars().take(256).collect::()); + // A band's first description, if it is given twice. + if let (Some(band), Some(description)) = (band, description) + && seen.insert(band) + { + pam.bands.push(PamBand { band, description: Some(description), ..PamBand::default() }); + } + } + (!pam.bands.is_empty()).then_some(pam) +} + #[cfg(test)] mod tests { use super::*; @@ -164,6 +200,29 @@ mod tests { assert!(parse(&deep).is_none()); } + #[test] + fn reads_band_descriptions_from_gdal_metadata() { + let xml = r#" + Area + NDVI 2024-07-01 + NDVI 2024-06-01 + again + 0.0001 + no band + past the last band + past the last band + +"#; + let pam = parse_gdal_metadata(xml).unwrap(); + assert_eq!(pam.bands.len(), 2, "{pam:?}"); + assert_eq!(pam.band(1).unwrap().description.as_deref(), Some("NDVI 2024-06-01")); + assert_eq!(pam.band(2).unwrap().description.as_deref(), Some("NDVI 2024-07-01")); + assert_eq!(pam.nodata(), None); + for text in ["", "not xml", "", "", "1"] { + assert!(parse_gdal_metadata(text).is_none(), "{text}"); + } + } + #[test] fn sidecar_names() { let p = sidecar_paths(Path::new("/d/scene.tif")); diff --git a/crates/engine/src/cmd/cube.rs b/crates/engine/src/cmd/cube.rs index dd8b65697..11527ec30 100644 --- a/crates/engine/src/cmd/cube.rs +++ b/crates/engine/src/cmd/cube.rs @@ -4,8 +4,9 @@ //! A [`TimeCube`] gives its dates and the slice at each one. The sources: //! - a Zarr variable with a time dimension (#554), read a slice at a time; //! - a stack of raster layers, each a date (from the `dates` given, or found in the layer's name -//! or file); -//! - a raster layer of several bands, each a date (from `dates`, or the bands' names). +//! or file), those on another grid than the first's resampled onto it as they are read; +//! - a raster layer of several bands, each a date (from `dates`, or the bands' names: the +//! descriptions GDAL writes in the file or its `.aux.xml`). //! //! `cube.trend` reads the slices into one block, cell by date, and works out each cell's trend //! with [`rustgis_analysis::trend`]: Sen's slope and the modified Mann–Kendall test. @@ -15,6 +16,7 @@ use std::sync::Arc; use std::sync::atomic::{AtomicBool, AtomicUsize, Ordering}; +use rustgis_analysis::rasterops; use rustgis_analysis::trend::{self, Test}; use rustgis_data::raster::{self, Raster}; use serde_json::{Value, json}; @@ -24,10 +26,10 @@ use crate::{CmdResult, CommandSpec as C, Session}; pub fn commands() -> Vec { vec![C::new("cube.trend", "Trend Analysis", "Analyze › Image Analysis", trend_cmd).params( - r#"{"cube"?: an entry of cube.catalog (its STAC API, collection, band, cell size and cloud cap, or its Zarr store and variable) | "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"#, + r#"{"cube"?: an entry of cube.catalog (its STAC API, collection, band, cell size and cloud cap, or its Zarr store and variable) | "layers"?: [raster layers, a date each; one on another grid is resampled (nearest) onto the first's] | "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#"{"cube"?: an entry of cube.catalog (as for cube.trend) | "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)"#, + r#"{"cube"?: an entry of cube.catalog (as for cube.trend) | "layers"?: [raster layers, a date each; one on another grid is resampled (nearest) onto the first's] | "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)"#, )] } @@ -100,6 +102,10 @@ pub struct Stack { items: Vec<(f64, Arc, usize)>, dates: Vec, label: String, + /// The grid the slices are put on, for a stack of layers on more than one: the first layer's. + /// Each slice on another is resampled onto it (nearest neighbour, as Raster Functions › + /// Resample does) as it is read, so only one is held at a time. + grid: Option>, } impl Stack { @@ -107,7 +113,7 @@ impl Stack { fn new(mut items: Vec<(f64, Arc, usize)>, label: String) -> Stack { items.sort_by(|a, b| a.0.total_cmp(&b.0)); let dates = items.iter().map(|i| i.0).collect(); - Stack { items, dates, label } + Stack { items, dates, label, grid: None } } } @@ -116,11 +122,23 @@ impl TimeCube for Stack { &self.dates } fn slice(&self, i: usize) -> Result, usize)>, String> { - self.items.get(i).map(|(_, r, b)| Some((r.clone(), *b))).ok_or_else(|| format!("no date {i}")) + let (_, r, b) = self.items.get(i).ok_or_else(|| format!("no date {i}"))?; + match &self.grid { + Some(g) if !one_grid(g, r) => { + let values = rasterops::resample_to(r, *b, g)?; + Ok(Some((Arc::new(rasterops::derive(g, &r.name, vec![values])?), 0))) + } + _ => Ok(Some((r.clone(), *b))), + } } fn label(&self) -> String { self.label.clone() } + /// The grid the slices are put on, known before any is read (or resampled). + fn size(&self) -> Option<(usize, usize)> { + let g = self.grid.as_ref().or_else(|| self.items.first().map(|(_, r, _)| r))?; + Some((g.width as usize, g.height as usize)) + } } /// A Zarr variable along its time dimension. @@ -269,6 +287,9 @@ pub fn cube_from(s: &Session, p: &Value, needs: Needs) -> Result> = None; + let mut grids = false; for (k, v) in list.iter().enumerate() { let id = s.layer_id(&json!({ "layer": v }))?; let l = s.map().layer(id).ok_or("no such layer")?; @@ -276,23 +297,29 @@ pub fn cube_from(s: &Session, p: &Value, needs: Needs) -> Result= r.bands.len() { return Err(format!("{} has {} bands", l.name, r.bands.len())); } - // The file's own name, not its folders' (`/data/2020/ndvi_a.tif` would date them all). - let file = l - .raster - .as_ref() - .and_then(|r| r.path.as_deref()) - .and_then(|f| std::path::Path::new(f).file_name()) - .map(|f| f.to_string_lossy().into_owned()) - .unwrap_or_default(); let date = match &dates { Some(d) => d.get(k).copied().ok_or("dates has a date for each layer, in their order")?, - None => date_in_text(&l.name) - .or_else(|| date_in_text(&file)) + None => date_of(&l.name, l.raster.as_ref().and_then(|r| r.path.as_deref())) .ok_or_else(|| format!("no date in the name of {}: give \"dates\"", l.name))?, }; if label.is_empty() { label = l.name.clone(); } + // A layer on another grid is resampled through longitude/latitude, so both need + // coordinates (a scan not yet georeferenced has none). + match &first { + None => first = Some(r.clone()), + Some(g) if !one_grid(g, &r) => { + g.require_coordinates()?; + r.require_coordinates()?; + // One that misses the first's area entirely would be a date with no value. + if !g.extent().intersects(&r.extent()) { + return Err(format!("{} lies outside the first layer's area, whose grid the stack is put on", l.name)); + } + grids = true; + } + Some(_) => {} + } // One layer a date: two on a day would count as two ordered values in the test. if items.iter().any(|(d, ..): &(f64, _, usize)| *d == date) { return Err(format!( @@ -303,7 +330,9 @@ pub fn cube_from(s: &Session, p: &Value, needs: Needs) -> Result Result) -> Option { + date_in_text(name).or_else(|| { + let file = std::path::Path::new(path?).file_name()?; + date_in_text(&file.to_string_lossy()) + }) +} + +/// The date a stack takes raster layer `layer` at when no `dates` are given, for the Workflows +/// pane's list of layers. +pub fn layer_date(s: &Session, layer: &str) -> Option { + let l = s.map().layer(s.layer_id(&json!({ "layer": layer })).ok()?)?; + date_of(&l.name, l.raster.as_ref().and_then(|r| r.path.as_deref())) +} + fn band_index(p: &Value) -> Result { match p.get("band").filter(|v| !v.is_null()) { None => Ok(0), @@ -549,11 +594,16 @@ impl TrendRun { format!("{cells} cells × {n} dates is too long a run: take a smaller area, a coarser cell size or fewer dates{how}") }) }; - // A cube that knows its size is refused before it reads anything. + // A cube that knows its size is refused before it reads anything; one read whole (a stack, + // whose slices may be resampled) by its values too, which a cube read in pieces checks + // piece by piece. let whole = self.cube.size().map(|(w, h)| w.saturating_mul(h)); if let Some(e) = whole.and_then(too_long) { return Err(e); } + if let Some(e) = whole.filter(|_| self.cube.pieces() <= 1).and_then(|cells| too_many_values(cells, n)) { + return Err(e); + } let need = self.min_valid.unwrap_or((n / 2).max(trend::MIN_VALUES)).max(trend::MIN_VALUES); let mut stitch = Stitch::default(); // The four outputs, a band each, filled piece by piece (never held interleaved whole). @@ -628,11 +678,8 @@ impl TrendRun { match &grid { None => { cells = vals.len(); - let values = (cells as u64).saturating_mul(n as u64); - if values > MAX_VALUES { - return Err(format!( - "{cells} cells × {n} dates is {values} values, more than {MAX_VALUES}: take a smaller area or fewer dates" - )); + if let Some(e) = too_many_values(cells, n) { + return Err(e); } if let Some(e) = too_long.and_then(|f| f(cells)) { return Err(e); @@ -737,6 +784,13 @@ impl TrendRun { } } +/// Why a block of `cells` × `n` dates is too large to read at once ([`MAX_VALUES`]), if it is. +fn too_many_values(cells: usize, n: usize) -> Option { + let values = (cells as u64).saturating_mul(n as u64); + (values > MAX_VALUES) + .then(|| format!("{cells} cells × {n} dates is {values} values, more than {MAX_VALUES}: take a smaller area or fewer dates")) +} + /// Reads piece `piece` of each of `steps`' slices and hands it to `place` with its place among the /// steps: in their order when read one at a time, as they come when the cube's /// [`TimeCube::readers`] read several at once. `cancel` stops it between slices. @@ -802,9 +856,15 @@ fn read_slices( Ok(()) } +/// Whether rasters `a` and `b` are on one grid: the same size and coordinate system, and +/// [`same_grid`]. +fn one_grid(a: &Raster, b: &Raster) -> bool { + (a.width, a.height, a.crs) == (b.width, b.height, b.crs) && same_grid(&a.transform, &b.transform) +} + /// Refuses slice `r`, at date `i`, when it isn't on the grid of the run's first slice `g`. fn off_grid(cube: &dyn TimeCube, i: usize, g: &Raster, r: &Raster) -> Result<(), String> { - if (g.width, g.height, g.crs) == (r.width, r.height, r.crs) && same_grid(&g.transform, &r.transform) { + if one_grid(g, r) { return Ok(()); } Err(format!( @@ -1771,6 +1831,11 @@ mod tests { assert_eq!(super::date_in_text("NDVI 2015"), day("2015-01-01")); assert_eq!(super::date_in_text("ndvi_v2015"), None); assert_eq!(super::date_in_text("MOD13Q1.A2024153"), day("2024-06-01")); + // A layer renamed keeps its file's date; a folder's doesn't date it. + assert_eq!(super::date_of("Sorriso NDVI", Some("/data/MOD13Q1.A2024153.tif")), day("2024-06-01")); + assert_eq!(super::date_of("ndvi_2015-06-01", Some("/data/MOD13Q1.A2024153.tif")), day("2015-06-01")); + assert_eq!(super::date_of("ndvi", Some("/data/2020/ndvi_a.tif")), None); + assert_eq!(super::date_of("ndvi", None), None); } /// A long run's thread that panics gives an error naming the command and the cause, and its @@ -1799,6 +1864,24 @@ mod tests { assert!(!super::same_grid(&a, &half)); } + /// A stack knows its size before any slice is read: the grid it is put on when its layers' grids + /// differ, whatever the earliest date's layer is; its first date's otherwise. + #[test] + fn a_stack_knows_its_size_before_reading() { + use super::TimeCube; + use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster}; + use std::sync::Arc; + let gt = GeoTransform { origin_x: -56.0, origin_y: -12.0, pixel_width: 0.01, pixel_height: -0.01 }; + let raster = |w: u32, h: u32| { + Arc::new(Raster::new("r", w, h, vec![vec![0.0; (w * h) as usize]], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap()) + }; + let (grid, early) = (raster(3, 2), raster(6, 4)); + let mut stack = super::Stack::new(vec![(2.0, grid.clone(), 0), (1.0, early, 0)], "r".into()); + assert_eq!(stack.size(), Some((6, 4)), "one grid: the first date's"); + stack.grid = Some(grid); + assert_eq!(stack.size(), Some((3, 2)), "the grid the slices are put on"); + } + /// 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"))] diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index e3a15e5da..a8e5600e4 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -10215,12 +10215,7 @@ fn a_trend_over_a_data_cube() { let r = s.run("cube.trend", &json!({"path": gap, "variable": "t", "test": "original"})).unwrap(); assert!(value(&s, &r["layers"][0], 0) > 1000.0, "{}", value(&s, &r["layers"][0], 0)); let _ = std::fs::remove_dir_all(gap_dir); - // Refused: too few layers, no dates, dates of another length, unknown options, another grid. - let other = Raster::new("ndvi_2030-01-01", 3, 3, vec![vec![0.0; 9]], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap(); - let bytes = geotiff::write(&other, 0).unwrap(); - crate::cmd::raster::add_raster(&mut s, "ndvi_2030-01-01", rustgis_map::RasterLayer::from_bytes(other, "o.tif", bytes)).unwrap(); - let mut mixed = names.clone(); - mixed.push("ndvi_2030-01-01".into()); + // Refused: too few layers, no dates, dates of another length, unknown options. for bad in [ json!({"layers": &names[..3]}), json!({"layers": names, "dates": ["2001", "2002"]}), @@ -10228,7 +10223,6 @@ fn a_trend_over_a_data_cube() { json!({"layers": names, "alpha": 0.7}), json!({"layers": names, "per": "week"}), json!({"layers": names, "start": "2030"}), - json!({"layers": mixed}), json!({"layer": names[0]}), ] { assert!(s.run("cube.trend", &bad).is_err(), "{bad}"); @@ -10291,6 +10285,238 @@ fn a_trend_over_a_data_cube() { } } +/// A multi-band Float32 GeoTIFF in longitude/latitude, its bands' values interleaved, as GDAL +/// writes it: `descriptions` in its GDAL_METADATA tag when given, and -9999 for no-data. +#[cfg(not(target_arch = "wasm32"))] +fn multiband_tiff(w: u32, h: u32, bands: &[Vec], origin: (f64, f64), px: f64, descriptions: Option<&[String]>) -> Vec { + let n = bands.len() as u16; + let mut out = vec![b'I', b'I', 42, 0, 0, 0, 0, 0]; + for p in 0..(w * h) as usize { + for b in bands { + out.extend(b[p].to_le_bytes()); + } + } + let bytes = out.len() as u32 - 8; + let shorts = |v: &[u16]| v.iter().flat_map(|x| x.to_le_bytes()).collect::>(); + let longs = |v: &[u32]| v.iter().flat_map(|x| x.to_le_bytes()).collect::>(); + let doubles = |v: &[f64]| v.iter().flat_map(|x| x.to_le_bytes()).collect::>(); + let ascii = |s: &str| s.bytes().chain([0]).collect::>(); + let mut e: Vec<(u16, u16, u32, Vec)> = vec![ + (256, 4, 1, longs(&[w])), + (257, 4, 1, longs(&[h])), + (258, 3, u32::from(n), shorts(&vec![32; n as usize])), + (259, 3, 1, shorts(&[1])), + (262, 3, 1, shorts(&[1])), + (273, 4, 1, longs(&[8])), + (277, 3, 1, shorts(&[n])), + (278, 4, 1, longs(&[h])), + (279, 4, 1, longs(&[bytes])), + (284, 3, 1, shorts(&[1])), + (339, 3, u32::from(n), shorts(&vec![3; n as usize])), + (33550, 12, 3, doubles(&[px, px, 0.0])), + (33922, 12, 6, doubles(&[0.0, 0.0, 0.0, origin.0, origin.1, 0.0])), + (34735, 3, 12, shorts(&[1, 1, 0, 2, 1024, 0, 1, 2, 2048, 0, 1, 4326])), + ]; + if let Some(d) = descriptions { + let items: String = + d.iter().enumerate().map(|(k, d)| format!(r#"{d}"#)).collect(); + let xml = ascii(&format!("{items}")); + e.push((42112, 2, xml.len() as u32, xml)); + } + let nodata = ascii("-9999"); + e.push((42113, 2, nodata.len() as u32, nodata)); + let ifd_at = out.len() + out.len() % 2; + out.resize(ifd_at, 0); + let values_at = ifd_at + 2 + 12 * e.len() + 4; + let mut values: Vec = Vec::new(); + let mut dir = (e.len() as u16).to_le_bytes().to_vec(); + for (tag, ty, count, bytes) in &mut e { + dir.extend(tag.to_le_bytes()); + dir.extend(ty.to_le_bytes()); + dir.extend(count.to_le_bytes()); + if bytes.len() <= 4 { + bytes.resize(4, 0); + dir.extend(bytes.iter()); + } else { + dir.extend(((values_at + values.len()) as u32).to_le_bytes()); + values.extend(bytes.iter()); + // Each value on a word boundary. + if values.len() % 2 == 1 { + values.push(0); + } + } + } + out[4..8].copy_from_slice(&(ifd_at as u32).to_le_bytes()); + out.extend(dir); + out.extend([0; 4]); + out.extend(values); + out +} + +/// A trend over GeoTIFF files (#579, step 3): a stack of files dated by their names (ISO and MODIS +/// dates), two of them on other grids, resampled onto the first's; and a multi-band file, a date a +/// band, dated by the band descriptions GDAL writes in it or by the dates given. The cells with a +/// planted trend have a significant slope; the noisy ones have a slope but no significant one. +#[cfg(not(target_arch = "wasm32"))] +#[test] +fn a_trend_over_stacked_and_multiband_geotiffs() { + use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster, geotiff}; + let mut s = Session::new(Project::default()); + let dir = std::env::temp_dir().join(format!("rustgis-cube-files-{}", std::process::id())); + std::fs::create_dir_all(&dir).unwrap(); + let (x0, y0, px) = (-56.0, -12.0, 0.01); + // Twelve years on a grid of 4 × 2 cells: 0 and 5 rise 0.02 a year, 3 and 7 fall 0.01 a year, + // 1, 2 and 6 are noise (the same noise turned over or around, never significant), and 4 has + // three values only, the file's no-data value on the other dates. + let noise = [0.03, -0.02, 0.05, -0.04, 0.01, 0.0, -0.05, 0.04, -0.01, 0.02, -0.03, 0.03]; + let cells = |k: usize| -> Vec { + let rise = 0.3 + 0.02 * k as f32; + let fall = 0.8 - 0.01 * k as f32; + let wobble = (0.001 * noise[k]) as f32; + vec![ + rise + wobble, + 0.5 + noise[k] as f32, + 0.5 - noise[k] as f32, + fall - wobble, + if k < 3 { 0.4 } else { -9999.0 }, + rise - wobble, + 0.5 + noise[11 - k] as f32, + fall + wobble, + ] + }; + let years: Vec = (2013..2025).collect(); + let mut names = Vec::new(); + for (k, y) in years.iter().enumerate() { + // Named by a date, or as MODIS names its files (1 July is day 182, 183 in a leap year). + let name = if k % 2 == 0 { format!("ndvi_{y}-07-01") } else { format!("MOD13Q1.A{y}{}", if y % 4 == 0 { 183 } else { 182 }) }; + let v = cells(k); + let (w, h, gt, data) = match k { + // Cells half as wide: each value four times. + 5 => { + let fine = (0..32).map(|i| v[(i / 16) * 4 + (i % 8) / 2]).collect(); + (8, 4, GeoTransform { origin_x: x0, origin_y: y0, pixel_width: px / 2.0, pixel_height: -px / 2.0 }, fine) + } + // A larger area, the grid's cells two columns and a row in from its corner. + 8 => { + let inside = |i: usize| (2..6).contains(&(i % 7)) && (1..3).contains(&(i / 7)); + let wide = (0..28).map(|i| if inside(i) { v[(i / 7 - 1) * 4 + i % 7 - 2] } else { 9.0 }).collect(); + (7, 4, GeoTransform { origin_x: x0 - 2.0 * px, origin_y: y0 + px, pixel_width: px, pixel_height: -px }, wide) + } + _ => (4, 2, GeoTransform { origin_x: x0, origin_y: y0, pixel_width: px, pixel_height: -px }, v), + }; + let r = Raster::new(&name, w, h, vec![data], gt, Crs::LonLat, Some(-9999.0), ColorHint::Gray, "Float32").unwrap(); + let path = dir.join(format!("{name}.tif")); + std::fs::write(&path, geotiff::write(&r, 0).unwrap()).unwrap(); + s.run("layer.addData", &json!({"path": path.display().to_string()})).unwrap(); + names.push(name); + } + let band = |s: &Session, id: &Value| { + s.map().layer(id.as_u64().unwrap() as u32).unwrap().raster.as_ref().unwrap().data.clone().unwrap().bands[0].data.clone() + }; + // The planted trends are found; the noise, and the cell with three values, are left empty. + let found = |s: &Session, r: &Value| { + let (slope, p, significant) = (band(s, &r["layers"][0]), band(s, &r["layers"][2]), band(s, &r["layers"][3])); + for (cell, want) in [(0, 0.02), (5, 0.02), (3, -0.01), (7, -0.01)] { + assert!((significant[cell] - want).abs() < 1e-3, "cell {cell}: {} (p {})", significant[cell], p[cell]); + } + for cell in [1, 2, 6] { + assert!(slope[cell].is_finite() && p[cell] > 0.05 && significant[cell].is_nan(), "cell {cell}: p {}", p[cell]); + } + assert!(slope[4].is_nan() && significant[4].is_nan()); + }; + let r = s.run("cube.trend", &json!({"layers": names})).unwrap(); + assert_eq!(r["dates"], 12, "{r}"); + found(&s, &r); + let name = |s: &Session, id: &Value| s.map().layer(id.as_u64().unwrap() as u32).unwrap().name.clone(); + assert_eq!(name(&s, &r["layers"][0]), "ndvi 2013-07-01–2024-07-01 Sen slope per year"); + // Every date counts in every cell, the two resampled ones included, and cell 4 has its three; + // the larger file's cells outside the grid are left out. + let count = s.run("cube.reduce", &json!({"layers": names, "stat": "count"})).unwrap(); + assert_eq!(band(&s, &count["layers"][0]), [12.0, 12.0, 12.0, 12.0, 3.0, 12.0, 12.0, 12.0]); + let max = s.run("cube.reduce", &json!({"layers": names, "stat": "max"})).unwrap(); + assert!(band(&s, &max["layers"][0]).iter().all(|v| *v < 1.0), "{:?}", band(&s, &max["layers"][0])); + // The first layer listed gives the grid: listed first, the finer file gives its 8 × 4 cells. + let mut fine_first = names.clone(); + fine_first.swap(0, 5); + let r = s.run("cube.trend", &json!({"layers": fine_first, "name": "fine"})).unwrap(); + let fine = s.map().layer(r["layers"][0].as_u64().unwrap() as u32).unwrap().raster.as_ref().unwrap().data.clone().unwrap(); + assert_eq!((fine.width, fine.height), (8, 4)); + // A scan not yet georeferenced, on another grid, can't be resampled: refused before the run. + let placed = Raster::new( + "scan", + 3, + 3, + vec![vec![0.5; 9]], + GeoTransform { origin_x: x0, origin_y: y0, pixel_width: px, pixel_height: -px }, + Crs::LonLat, + None, + ColorHint::Gray, + "Float32", + ) + .unwrap(); + let bytes = geotiff::write(&placed, 0).unwrap(); + let mut scan = placed; + scan.georeferenced = false; + crate::cmd::raster::add_raster(&mut s, "scan_2030-07-01", rustgis_map::RasterLayer::from_bytes(scan, "scan.tif", bytes)).unwrap(); + let mut with_scan = names.clone(); + with_scan.push("scan_2030-07-01".into()); + let e = s.run("cube.trend", &json!({"layers": with_scan})).unwrap_err(); + assert!(e.contains("no coordinates"), "{e}"); + // A layer away from the first's area would give its date no value: refused. + let gt_away = GeoTransform { origin_x: x0 + 10.0, origin_y: y0, pixel_width: px, pixel_height: -px }; + let away = Raster::new("away_2030-07-01", 3, 3, vec![vec![0.5; 9]], gt_away, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap(); + let bytes = geotiff::write(&away, 0).unwrap(); + crate::cmd::raster::add_raster(&mut s, "away_2030-07-01", rustgis_map::RasterLayer::from_bytes(away, "away.tif", bytes)).unwrap(); + let mut with_away = names.clone(); + with_away.push("away_2030-07-01".into()); + let e = s.run("cube.trend", &json!({"layers": with_away})).unwrap_err(); + assert!(e.contains("outside the first layer's area"), "{e}"); + // A layer in another coordinate system (Web Mercator, 500 m cells) over the grid is resampled + // onto it too, through longitude/latitude: every cell gets its date. + let merc = |lon: f64, lat: f64| { + (6_378_137.0 * lon.to_radians(), 6_378_137.0 * (std::f64::consts::FRAC_PI_4 + lat.to_radians() / 2.0).tan().ln()) + }; + let (mx, my) = merc(x0 - 0.05, y0 + 0.03); + let gt_merc = GeoTransform { origin_x: mx, origin_y: my, pixel_width: 500.0, pixel_height: -500.0 }; + let m = + Raster::new("merc_2030-07-01", 40, 40, vec![vec![0.5; 1600]], gt_merc, Crs::WebMercator, None, ColorHint::Gray, "Float32").unwrap(); + let bytes = geotiff::write(&m, 0).unwrap(); + crate::cmd::raster::add_raster(&mut s, "merc_2030-07-01", rustgis_map::RasterLayer::from_bytes(m, "merc.tif", bytes)).unwrap(); + let mut with_merc = names.clone(); + with_merc.push("merc_2030-07-01".into()); + let count = s.run("cube.reduce", &json!({"layers": with_merc, "stat": "count", "name": "merc"})).unwrap(); + assert_eq!(band(&s, &count["layers"][0]), [13.0, 13.0, 13.0, 13.0, 4.0, 13.0, 13.0, 13.0]); + // A multi-band file, a date a band, its bands named by the descriptions in it. + let bands: Vec> = (0..12).map(cells).collect(); + let described: Vec = years.iter().map(|y| format!("NDVI {y}-07-01")).collect(); + let path = dir.join("ndvi_bands.tif"); + std::fs::write(&path, multiband_tiff(4, 2, &bands, (x0, y0), px, Some(&described))).unwrap(); + s.run("layer.addData", &json!({"path": path.display().to_string()})).unwrap(); + let file = s.map().layer_by_name("ndvi_bands").unwrap().raster.as_ref().unwrap().data.clone().unwrap(); + assert_eq!((file.bands.len(), file.bands[0].name.as_str()), (12, "NDVI 2013-07-01")); + let r = s.run("cube.trend", &json!({"layer": "ndvi_bands"})).unwrap(); + assert_eq!(r["dates"], 12, "{r}"); + found(&s, &r); + // An `.aux.xml` next to the file names the bands it describes, the file's own tag the others. + let path = dir.join("ndvi_side.tif"); + std::fs::write(&path, multiband_tiff(4, 2, &bands, (x0, y0), px, Some(&described))).unwrap(); + let sidecar = r#"NDVI 2012-07-01"#; + std::fs::write(dir.join("ndvi_side.tif.aux.xml"), sidecar).unwrap(); + s.run("layer.addData", &json!({"path": path.display().to_string()})).unwrap(); + let side = s.map().layer_by_name("ndvi_side").unwrap().raster.as_ref().unwrap().data.clone().unwrap(); + assert_eq!((side.bands[0].name.as_str(), side.bands[1].name.as_str()), ("NDVI 2012-07-01", "NDVI 2014-07-01")); + // Without descriptions, the dates are given. + let path = dir.join("ndvi_plain.tif"); + std::fs::write(&path, multiband_tiff(4, 2, &bands, (x0, y0), px, None)).unwrap(); + s.run("layer.addData", &json!({"path": path.display().to_string()})).unwrap(); + let e = s.run("cube.trend", &json!({"layer": "ndvi_plain"})).unwrap_err(); + assert!(e.contains("give \"dates\""), "{e}"); + let typed: Vec = years.iter().map(|y| format!("{y}-07-01")).collect(); + let r = s.run("cube.trend", &json!({"layer": "ndvi_plain", "dates": typed, "name": "typed"})).unwrap(); + found(&s, &r); + let _ = std::fs::remove_dir_all(dir); +} + /// A data cube summarised over time (#579): each statistic over a stack of layers, a year or a /// month of the year at a time, a cell needing a share of its group's dates, a raster's bands, and /// a Zarr variable as `zarr.reduce` sums it up; each output kept in the project. diff --git a/crates/ui-egui/src/analysispane.rs b/crates/ui-egui/src/analysispane.rs index 1155b103c..b7ce537b6 100644 --- a/crates/ui-egui/src/analysispane.rs +++ b/crates/ui-egui/src/analysispane.rs @@ -72,6 +72,9 @@ pub struct TrendForm { pub end: String, pub months: String, pub min_valid: String, + /// The dates typed, a layer's (or band's) each in their order; empty to read them from the + /// names. + pub dates: String, #[cfg(not(target_arch = "wasm32"))] pub running: Option, } @@ -929,10 +932,15 @@ fn trend(app: &mut RustGisApp, ui: &mut Ui, rasters: &[String]) { if f.test.is_empty() { (f.test, f.seasons, f.alpha) = ("hamedRao".into(), 12, "0.05".into()); } + let was_stack = f.stack; ui.horizontal(|ui| { ui.radio_value(&mut f.stack, true, tl!("Layers, a date each")); ui.radio_value(&mut f.stack, false, tl!("One layer: a Zarr variable or a raster's bands")); }); + // Dates typed for one kind of series don't carry over to the other. + if f.stack != was_stack { + f.dates.clear(); + } if f.stack { f.ticked.retain(|n| rasters.contains(n)); ui.horizontal(|ui| { @@ -946,7 +954,7 @@ fn trend(app: &mut RustGisApp, ui: &mut Ui, rasters: &[String]) { egui::ScrollArea::vertical().id_salt("trend_layers").max_height(160.0).show(ui, |ui| { for n in rasters { let mut on = f.ticked.contains(n); - let date = rustgis_engine::cmd::cube::date_in_text(n) + let date = rustgis_engine::cmd::cube::layer_date(&app.session, n) .map_or_else(|| tl!("no date in the name").to_string(), |d| rustgis_engine::data::time::format(d, "YYYY-MM-DD")); if ui.checkbox(&mut on, format!("{n} ({date})")).changed() { if on { @@ -962,6 +970,9 @@ fn trend(app: &mut RustGisApp, ui: &mut Ui, rasters: &[String]) { } let f = &mut app.workflow.trend; egui::Grid::new("trend_grid").num_columns(2).spacing(vec2(8.0, 6.0)).show(ui, |ui| { + let input_label = ui.label(tl!("Dates")); + ui.add(egui::TextEdit::singleline(&mut f.dates).hint_text(tl!("from the names")).desired_width(240.0)).labelled_by(input_label.id); + ui.end_row(); let input_label = ui.label(tl!("Test")); let tests = [ ("hamedRao", "Modified Mann–Kendall (Hamed & Rao)"), @@ -1018,7 +1029,10 @@ fn trend(app: &mut RustGisApp, ui: &mut Ui, rasters: &[String]) { } let mut p = json!({"test": f.test, "per": if f.decade { "decade" } else { "year" }}); if f.stack { - p["layers"] = json!(f.ticked); + // In the list's order, not the order they were ticked: the dates typed go with them in + // that order, and the first gives the grid. + let ticked: Vec<&String> = rasters.iter().filter(|n| f.ticked.contains(n)).collect(); + p["layers"] = json!(ticked); } else { p["layer"] = json!(app.workflow.layer); } @@ -1046,6 +1060,18 @@ fn trend(app: &mut RustGisApp, ui: &mut Ui, rasters: &[String]) { p[k] = json!(v.trim()); } } + // Dates typed, a layer's or band's each in their order, apart by commas or semicolons (a date + // can hold a space before its time): checked by the command, which says which it can't read. + let parts: Vec<&str> = f.dates.split([',', ';']).map(str::trim).collect(); + // An empty one before the last would move the dates after it onto other layers; a comma at + // the end is let be. + if parts.iter().rev().skip(1).any(|d| d.is_empty()) { + bad = Some(tl!("A date is missing between two commas").to_string()); + } + let dates: Vec<&str> = parts.into_iter().filter(|d| !d.is_empty()).collect(); + if !dates.is_empty() { + p["dates"] = json!(dates); + } let typed: Vec<&str> = f.months.split([',', ' ']).map(str::trim).filter(|m| !m.is_empty()).collect(); let months: Vec = typed.iter().filter_map(|m| m.parse::().ok()).map(|m| json!(m)).collect(); if months.len() != typed.len() { diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 8ae08fb11..3364ebf2b 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -457,12 +457,18 @@ between every two dates, which holds up against outliers such as a cloud left in and a **Mann–Kendall** test of whether the trend is significant. - **The series** (a data cube) is any of: - - raster layers, a date each, dated by their names (`ndvi_2024-06-01`, MODIS's `A2024153`, - `2024-06`, a year) or by the `dates` given: `{"layers": ["ndvi_2015", "ndvi_2016", …]}`; + - raster layers, a date each, dated by their names or their files' (`ndvi_2024-06-01`, MODIS's + `A2024153`, `2024-06`, a year) or by the `dates` given: `{"layers": ["ndvi_2015", "ndvi_2016", + …]}`. A layer on another grid than the first one listed is resampled onto it (nearest + neighbour, as Raster Functions › Resample does) as it is read, so a scene with another cell size + or coordinate system takes part. The first layer sets the cell size: a finer layer gives one of + its pixels to each cell, not their mean. Cells a layer doesn't cover are empty that date. The + Workflows pane's Dates field takes the dates separated by commas; - a Zarr variable along its time dimension, by its layer or its store, as for Summarise over time (see [Data](data.md)): `{"path": "air.zarr", "variable": "t2m"}`; - a raster of several bands, a date each, as Earth Engine and openEO export, dated by the bands' - names or `dates`; + names (the descriptions GDAL writes in the file, or in an `.aux.xml` next to it: `NDVI + 2024-06-01`) or `dates`: `{"layer": "ndvi_2015_2024", "dates": ["2015", "2016", …]}`; - a STAC collection's band (desktop), read in place from its cloud-optimised GeoTIFFs over an area, as Add as Time Series finds them (see [Web data](web-data.md)): `{"url": "https://data.inpe.br/bdc/stac/v1", "collection": "mod13q1-6.1", "asset": "NDVI", "bbox": [-46, @@ -486,8 +492,6 @@ and a **Mann–Kendall** test of whether the trend is significant. "bdc-mod13q1-ndvi", "bbox": [-46, -10.95, -40.35, -2.7], "start": "2001", "end": "2025"}` fills in the collection's API, band, cell size and cloud cap, or a Zarr entry's store and variable. - - The rasters have to be on one grid; put them on one first with Raster Functions › Resample. - **The test:** NDVI and temperature series are autocorrelated, which makes the plain test report too many significant cells, so the default is the **modified Mann–Kendall** test of **Hamed & Rao (1998)**: the variance of S is corrected by the autocorrelation of the detrended series' diff --git a/docs/user-guide/data.md b/docs/user-guide/data.md index 11965849e..55d8b481e 100644 --- a/docs/user-guide/data.md +++ b/docs/user-guide/data.md @@ -694,7 +694,9 @@ A TIFF or PNG without coordinates opens as a scan to place with **Georeference** file next to it (as a Shapefile does); otherwise, and always for uploads, it must be in longitude/latitude degrees. On the desktop, a GDAL `.aux.xml` file next to a raster (`dem.tif.aux.xml` or `dem.aux.xml`) names its bands from their descriptions, gives the no-data value when the file -has none, and supplies the band statistics when only an overview of the raster is read. +has none, and supplies the band statistics when only an overview of the raster is read. The band +descriptions GDAL and rasterio write inside a GeoTIFF (its `GDAL_METADATA` tag) name its bands +too, a sidecar's taking precedence. Exporting a raster opened from an overview writes that overview; the results of adding and of exporting it (`layer.addData`, `share.exportRaster`) give the full size as `overviewOf`, so scripts and the MCP server see it too.