From 84158383426cb72b5ad25d65a56af1a35a6551f3 Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Sat, 24 Oct 2026 15:10:21 -0300 Subject: [PATCH 1/7] feat(analysis): project a raster onto a chosen grid (raster.project, #821 part 1) raster.project puts a raster on a grid of another coordinate system, cell size, origin and extent (a polygon layer, a box or the view), or exactly on another raster's grid, so rasters from different sources line up cell for cell. Cells are sampled at their centres (nearest, bilinear, cubic), or take every pixel whose centre falls in them (average, mode), which is what a coarser grid of elevation or land-cover classes should hold. It builds on warp.rs: the footprint, the row mapping and the sampling are shared with Georeference, which is unchanged. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/warp.rs | 690 ++++++++++++++++++++++++++----- crates/assistant/src/tests.rs | 8 +- crates/engine/src/cmd/imagery.rs | 111 +++++ crates/engine/src/cmd/plugins.rs | 5 +- crates/engine/src/tests.rs | 84 +++- docs/user-guide/analysis.md | 37 +- 6 files changed, 817 insertions(+), 118 deletions(-) diff --git a/crates/analysis/src/warp.rs b/crates/analysis/src/warp.rs index 76db77772..2fffd59f3 100644 --- a/crates/analysis/src/warp.rs +++ b/crates/analysis/src/warp.rs @@ -1,4 +1,5 @@ -//! Warping a raster onto a north-up grid through a control-point fit ([`crate::georef::Fit`]). +//! Warping a raster onto a north-up grid through a control-point fit ([`crate::georef::Fit`]), and +//! projecting one onto a grid of another coordinate system, cell size and extent ([`project`]). //! //! `GeoTransform` has no rotation terms, so an image that a fit turns, shears or bends can only be //! placed by resampling it. Every cell of the output grid is taken back through the fit (after @@ -206,47 +207,15 @@ pub fn warp(r: &Raster, fit: &Fit, o: &Options) -> Result { let (sx, sy) = r.overview.map_or((1.0, 1.0), |(fw, fh)| (f64::from(fw) / w, f64::from(fh) / h)); let place = |c: f64, row: f64| fit.to_map(c * sx, row * sy).and_then(|(x, y)| re.forward(x, y)); - // The footprint: the image's edge and a lattice inside it, placed (a pole inside the image - // is no extreme of its edge). - let (mut min_x, mut min_y, mut max_x, mut max_y) = (f64::MAX, f64::MAX, f64::MIN, f64::MIN); - let mut placed = 0; - let mut add = |c: f64, row: f64| { - if let Some((x, y)) = place(c, row) { - (min_x, min_y, max_x, max_y) = (min_x.min(x), min_y.min(y), max_x.max(x), max_y.max(y)); - placed += 1; - } - }; - for i in 0..=EDGE_POINTS { - let t = i as f64 / EDGE_POINTS as f64; - for (c, row) in [(t * w, 0.0), (t * w, h), (0.0, t * h), (w, t * h)] { - add(c, row); - } - } - for i in 0..=INTERIOR_POINTS { - for j in 0..=INTERIOR_POINTS { - add(i as f64 / INTERIOR_POINTS as f64 * w, j as f64 / INTERIOR_POINTS as f64 * h); - } - } - if placed < 4 || ![min_x, min_y, max_x, max_y].iter().all(|v| v.is_finite()) { - return Err("the image's edge can't be placed in the output coordinate system".into()); - } + let [min_x, min_y, max_x, max_y] = footprint(w, h, &place).ok_or("the image's edge can't be placed in the output coordinate system")?; if o.crs.is_geographic() && max_x - min_x > 180.0 { return Err("the warped image spans more than 180° of longitude (it crosses the antimeridian, or the fit sends part of it far away): give a projected crs".into()); } if max_x <= min_x || max_y <= min_y { return Err("the fit puts the image on a line, with no area".into()); } - - let nb = r.bands.len(); - let alpha_src = match r.color { - ColorHint::Rgb { alpha: true } if nb >= 4 => Some(3), - ColorHint::GrayAlpha if nb >= 2 => Some(1), - _ => None, - }; - let rgb = matches!(r.color, ColorHint::Rgb { .. }) && nb >= 3; - let gray = matches!(r.color, ColorHint::Gray) && nb == 1 && r.data_type == "UInt8"; - let add_alpha = (rgb || gray) && alpha_src.is_none() && r.nodata.is_none(); - let out_bands = nb + usize::from(add_alpha); + let canvas = Canvas::of(r, true); + let out_bands = canvas.out_bands; // The cell: asked for, or the image pixel's size on the ground (the root of its area, where // the fit turns or skews it), coarser if the grid would hold too many samples. @@ -263,15 +232,7 @@ pub fn warp(r: &Raster, fit: &Fit, o: &Options) -> Result { c } None => { - let (cx, cy) = (w / 2.0, h / 2.0); - let (a, b, d) = (place(cx, cy), place(cx + 1.0, cy), place(cx, cy + 1.0)); - let (Some(a), Some(b), Some(d)) = (a, b, d) else { - return Err("the source pixel's size can't be worked out: give a cellSize".into()); - }; - let mut cell = ((b.0 - a.0) * (d.1 - a.1) - (d.0 - a.0) * (b.1 - a.1)).abs().sqrt(); - if !(cell.is_finite() && cell > 0.0) { - return Err("the source pixel's size can't be worked out: give a cellSize".into()); - } + let mut cell = pixel_size(w, h, &place).ok_or("the source pixel's size can't be worked out: give a cellSize")?; // Each step at least 1% coarser, so a grid rounded up past the limit still shrinks. for _ in 0..64 { let n = cells(cell); @@ -287,20 +248,7 @@ pub fn warp(r: &Raster, fit: &Fit, o: &Options) -> Result { } }; let (ow, oh) = cells(cell); - let (ow, oh) = (ow as usize, oh as usize); - - let method = if matches!(r.color, ColorHint::Palette(_)) { Resampling::Nearest } else { o.resampling }; - let sampler = Sampler { r, method, alpha: alpha_src, round: !r.data_type.starts_with("Float") }; - // Cells with no value: no-data where the raster has some. Colours with alpha are zero there. - let clear = alpha_src.is_some() || add_alpha; - let outside = if clear { 0.0 } else { r.nodata.map_or(f32::NAN, |d| d as f32) }; - // Alpha is read as 0–255 whatever the pixel type. - let opaque = 255.0; - - let mut bands: Vec> = vec![vec![outside; ow * oh]; out_bands]; - let mut coords: Vec> = vec![None; ow]; - let mut nodes: Vec<(usize, Option<(f64, f64)>)> = Vec::new(); - let mut px = vec![f32::NAN; nb]; + let g = Grid { crs: o.crs, min_x, max_y, cell_x: cell, cell_y: cell, width: ow as usize, height: oh as usize }; // The image position of an output coordinate. let source = |x: f64, y: f64| -> Option<(f64, f64)> { let (px, py) = re.back(x, y)?; @@ -308,77 +256,474 @@ pub fn warp(r: &Raster, fit: &Fit, o: &Options) -> Result { let (c, row) = (c / sx, row / sy); (c.is_finite() && row.is_finite()).then_some((c, row)) }; - for oy in 0..oh { - let y = max_y - (oy as f64 + 0.5) * cell; - let x_at = |ox: usize| min_x + (ox as f64 + 0.5) * cell; - nodes.clear(); - let mut ox = 0; - loop { - nodes.push((ox, source(x_at(ox), y))); - if ox + 1 >= ow { - break; - } - ox = (ox + STEP).min(ow - 1); - } - for pair in nodes.windows(2) { - let [(a, pa), (b, pb)] = pair else { continue }; - let linear = match (pa, pb) { - (Some(pa), Some(pb)) if b - a >= 2 => { - let mid = a + (b - a) / 2; - let t = (mid - a) as f64 / (b - a) as f64; - let guess = (pa.0 + (pb.0 - pa.0) * t, pa.1 + (pb.1 - pa.1) * t); - source(x_at(mid), y).is_some_and(|m| (m.0 - guess.0).abs().max((m.1 - guess.1).abs()) <= TOLERANCE) - } - (Some(_), Some(_)) => true, - _ => false, - }; - for i in *a..*b { - let t = (i - a) as f64 / (b - a) as f64; - let c = match (linear, pa, pb) { - (true, Some(pa), Some(pb)) => Some((pa.0 + (pb.0 - pa.0) * t, pa.1 + (pb.1 - pa.1) * t)), - _ => source(x_at(i), y), - }; - if let Some(slot) = coords.get_mut(i) { - *slot = c; - } + sampled(r, &g, canvas, o.resampling, &source, o.name.clone()) +} + +/// The bounds `[min x, min y, max x, max y]` of a `w` × `h` image whose pixel positions `place` +/// puts in output coordinates: its edge and a lattice inside it (a pole inside the image is no +/// extreme of its edge). `None` when too little of it can be placed. +fn footprint(w: f64, h: f64, place: &dyn Fn(f64, f64) -> Option<(f64, f64)>) -> Option<[f64; 4]> { + let (mut min_x, mut min_y, mut max_x, mut max_y) = (f64::MAX, f64::MAX, f64::MIN, f64::MIN); + let mut placed = 0; + let mut add = |c: f64, row: f64| { + if let Some((x, y)) = place(c, row) { + (min_x, min_y, max_x, max_y) = (min_x.min(x), min_y.min(y), max_x.max(x), max_y.max(y)); + placed += 1; + } + }; + for i in 0..=EDGE_POINTS { + let t = i as f64 / EDGE_POINTS as f64; + for (c, row) in [(t * w, 0.0), (t * w, h), (0.0, t * h), (w, t * h)] { + add(c, row); + } + } + for i in 0..=INTERIOR_POINTS { + for j in 0..=INTERIOR_POINTS { + add(i as f64 / INTERIOR_POINTS as f64 * w, j as f64 / INTERIOR_POINTS as f64 * h); + } + } + (placed >= 4 && [min_x, min_y, max_x, max_y].iter().all(|v| v.is_finite())).then_some([min_x, min_y, max_x, max_y]) +} + +/// About the size of a pixel of a `w` × `h` image in output coordinates, at its centre: the root +/// of its area, where `place` turns or skews it. +fn pixel_size(w: f64, h: f64, place: &dyn Fn(f64, f64) -> Option<(f64, f64)>) -> Option { + let (cx, cy) = (w / 2.0, h / 2.0); + let (a, b, d) = (place(cx, cy)?, place(cx + 1.0, cy)?, place(cx, cy + 1.0)?); + let cell = ((b.0 - a.0) * (d.1 - a.1) - (d.0 - a.0) * (b.1 - a.1)).abs().sqrt(); + (cell.is_finite() && cell > 0.0).then_some(cell) +} + +/// A north-up grid: its coordinate system, the top left corner of its top left cell, its cell +/// size each way and how many cells it has. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Grid { + pub crs: Crs, + pub min_x: f64, + pub max_y: f64, + pub cell_x: f64, + pub cell_y: f64, + pub width: usize, + pub height: usize, +} + +impl Grid { + /// The grid `r` is on. + pub fn of(r: &Raster) -> Grid { + let t = r.transform; + let (w, h) = (f64::from(r.width), f64::from(r.height)); + Grid { + crs: r.crs, + min_x: t.origin_x.min(t.origin_x + w * t.pixel_width), + max_y: t.origin_y.max(t.origin_y + h * t.pixel_height), + cell_x: t.pixel_width.abs(), + cell_y: t.pixel_height.abs(), + width: r.width as usize, + height: r.height as usize, + } + } + + /// Cells of `cell` (x, y) in `crs` covering `bounds` (`[min x, min y, max x, max y]`), their + /// corners on `origin` plus whole cells, so grids with the same origin and cell line up. + pub fn snapped(crs: Crs, bounds: [f64; 4], cell: (f64, f64), origin: (f64, f64)) -> Result { + let [x0, y0, x1, y1] = bounds; + let (cx, cy) = cell; + if !(cx.is_finite() && cy.is_finite() && cx > 0.0 && cy > 0.0) { + return Err("cellSize must be a positive number".into()); + } + if !(origin.0.is_finite() && origin.1.is_finite()) { + return Err("origin is [x, y] in the output's units".into()); + } + if !(bounds.iter().all(|v| v.is_finite()) && x0 < x1 && y0 < y1) { + return Err("the area has no extent in the output coordinate system".into()); + } + let min_x = origin.0 + ((x0 - origin.0) / cx).floor() * cx; + let max_x = origin.0 + ((x1 - origin.0) / cx).ceil() * cx; + let min_y = origin.1 + ((y0 - origin.1) / cy).floor() * cy; + let max_y = origin.1 + ((y1 - origin.1) / cy).ceil() * cy; + let (w, h) = (((max_x - min_x) / cx).round().max(1.0), ((max_y - min_y) / cy).round().max(1.0)); + if w * h > MAX_SAMPLES as f64 { + return Err(format!( + "a cell size of {cx} makes a grid of {} × {} cells, more than {MAX_SAMPLES}: give a larger cellSize or a smaller extent", + w as u64, h as u64 + )); + } + Ok(Grid { crs, min_x, max_y, cell_x: cx, cell_y: cy, width: w as usize, height: h as usize }) + } + + /// The top left corner's transform. + fn transform(&self) -> GeoTransform { + GeoTransform { origin_x: self.min_x, origin_y: self.max_y, pixel_width: self.cell_x, pixel_height: -self.cell_y } + } +} + +/// The bounds in `crs` of `r`'s footprint, and about the size of its pixel there. +pub fn footprint_in(r: &Raster, crs: Crs) -> Result<([f64; 4], f64), String> { + r.require_coordinates()?; + let re = Reproject::new(r.crs, crs)?; + let t = r.transform; + let place = |c: f64, row: f64| { + let p = t.to_world(c, row); + re.forward(p.x, p.y) + }; + let (w, h) = (f64::from(r.width), f64::from(r.height)); + let bounds = footprint(w, h, &place).ok_or_else(|| format!("{} can't be placed in {}", r.name, crs.label()))?; + let cell = pixel_size(w, h, &place).ok_or_else(|| format!("the size of {}'s pixel can't be worked out: give a cellSize", r.name))?; + Ok((bounds, cell)) +} + +/// The bounds in `crs` of a longitude/latitude box `[west, south, east, north]`, its edges +/// followed (a projection bends them). +pub fn bounds_in(lonlat: [f64; 4], crs: Crs) -> Result<[f64; 4], String> { + let [west, south, east, north] = lonlat; + let re = Reproject::new(Crs::LonLat, crs)?; + let place = |u: f64, v: f64| re.forward(west + u * (east - west), north - v * (north - south)); + footprint(1.0, 1.0, &place).ok_or_else(|| format!("the area can't be placed in {}", crs.label())) +} + +/// How [`project`] fills a cell: sampled at its centre, or from every source pixel whose centre +/// falls in it. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Method { + /// Sampled at the cell's centre. + Sample(Resampling), + /// The mean of the source pixels in the cell. + Average, + /// The most frequent value of the source pixels in the cell, the least of those tied. + Mode, +} + +impl Method { + pub fn parse(s: &str) -> Option { + match s.trim().to_ascii_lowercase().as_str() { + "average" | "mean" => Some(Method::Average), + "mode" | "majority" => Some(Method::Mode), + other => Resampling::parse(other).map(Method::Sample), + } + } + + pub fn id(self) -> &'static str { + match self { + Method::Sample(r) => r.id(), + Method::Average => "average", + Method::Mode => "mode", + } + } +} + +/// Projects `r` onto the grid `g`, through longitude/latitude when their coordinate systems +/// differ: each cell sampled at its centre, or (average, mode) taken from every pixel of `r` +/// whose centre falls in it. A cell no pixel centre falls in (a grid finer than the raster) takes +/// the pixel under its centre. Cells off the raster are its no-data (else NaN, or 0 in a +/// whole-number type), or see-through for a colour image with neither alpha nor no-data, which +/// gains an alpha band for it. +pub fn project(r: &Raster, g: &Grid, method: Method, name: String) -> Result { + r.require_coordinates()?; + let canvas = Canvas::of(r, false); + let out_bands = canvas.out_bands; + // The samples held, and the buffers a row needs. + if g.width as f64 * (g.height as f64 * out_bands as f64 + COLUMN_COST) > MAX_SAMPLES as f64 { + return Err(format!( + "{} × {} cells of {} band(s) are more than {MAX_SAMPLES} samples: give a larger cellSize or a smaller extent", + g.width, g.height, out_bands + )); + } + let re = Reproject::new(r.crs, g.crs)?; + let t = r.transform; + let source = |x: f64, y: f64| -> Option<(f64, f64)> { + let (px, py) = re.back(x, y)?; + let (c, row) = t.to_pixel(Coord::new(px, py)); + (c.is_finite() && row.is_finite()).then_some((c, row)) + }; + match method { + Method::Sample(m) => sampled(r, g, canvas, m, &source, name), + Method::Average | Method::Mode => aggregated(r, g, canvas, method == Method::Mode, &source, name), + } +} + +/// How an output raster holds a source's bands. +struct Canvas { + /// The source's bands. + nb: usize, + /// Its alpha band, if it has one. + alpha_src: Option, + /// Whether an alpha band is added: a colour image (or, for a picture, a grey byte one) with + /// neither alpha nor no-data, so the cells outside the image are see-through. + add_alpha: bool, + out_bands: usize, + /// Whether a cell with no value is see-through (zero) rather than no-data. + clear: bool, + /// What a cell with no value holds. + outside: f32, + color: ColorHint, +} + +impl Canvas { + /// How `r` is held: as a `picture` (a scan being georeferenced), a grey byte image gains + /// alpha as a colour one does; as data, it keeps its one band, and a cell with no value is + /// its no-data, else NaN in a floating-point type and 0 in a whole-number one, as GDAL leaves + /// it. + fn of(r: &Raster, picture: bool) -> Canvas { + let nb = r.bands.len(); + let alpha_src = match r.color { + ColorHint::Rgb { alpha: true } if nb >= 4 => Some(3), + ColorHint::GrayAlpha if nb >= 2 => Some(1), + _ => None, + }; + let rgb = matches!(r.color, ColorHint::Rgb { .. }) && nb >= 3; + let gray = picture && matches!(r.color, ColorHint::Gray) && nb == 1 && r.data_type == "UInt8"; + let add_alpha = (rgb || gray) && alpha_src.is_none() && r.nodata.is_none(); + // Cells with no value: no-data where the raster has some. Colours with alpha are zero there. + let clear = alpha_src.is_some() || add_alpha; + let none = if picture || r.data_type.starts_with("Float") { f32::NAN } else { 0.0 }; + let outside = if clear { 0.0 } else { r.nodata.map_or(none, |d| d as f32) }; + let color = match (add_alpha, rgb) { + (true, true) => ColorHint::Rgb { alpha: true }, + (true, false) => ColorHint::GrayAlpha, + _ => r.color.clone(), + }; + Canvas { nb, alpha_src, add_alpha, out_bands: nb + usize::from(add_alpha), clear, outside, color } + } + + /// Cell `at` of `bands` set to `values` (no value where NaN), and opaque when an alpha band + /// was added. + fn set(&self, bands: &mut [Vec], at: usize, values: &[f32]) { + for (band, v) in bands.iter_mut().zip(values) { + if let Some(o) = band.get_mut(at) { + // A cell with no value is no-data (or see-through) like one outside. + *o = if v.is_nan() { if self.clear { 0.0 } else { self.outside } } else { *v }; } } - if let Some((last, p)) = nodes.last() - && let Some(c) = coords.get_mut(*last) + // Alpha is read as 0–255 whatever the pixel type. + if self.add_alpha + && let Some(o) = bands.get_mut(self.nb).and_then(|b| b.get_mut(at)) { - *c = *p; + *o = 255.0; + } + } + + /// The output raster on `g`, its bands named as the source's. + fn finish(&self, r: &Raster, g: &Grid, bands: Vec>, name: String) -> Result { + let mut out = Raster::new( + name, + g.width as u32, + g.height as u32, + bands, + g.transform(), + g.crs, + r.nodata, + self.color.clone(), + r.data_type.clone(), + ) + .map_err(|e| e.to_string())?; + for (i, band) in out.bands.iter_mut().enumerate() { + band.name = r.bands.get(i).map_or_else(|| "Alpha".to_string(), |b| b.name.clone()); } + Ok(out) + } +} + +/// The image positions of the `out.len()` points `x_at(i)` along the row `y` of output +/// coordinates: exact every [`STEP`] points, and linear between where the midpoint of a stretch +/// is within [`TOLERANCE`] of exact (otherwise exact throughout the stretch). +fn row_positions( + x_at: &dyn Fn(usize) -> f64, + y: f64, + source: &dyn Fn(f64, f64) -> Option<(f64, f64)>, + nodes: &mut Vec<(usize, Option<(f64, f64)>)>, + out: &mut [Option<(f64, f64)>], +) { + let n = out.len(); + if n == 0 { + return; + } + nodes.clear(); + let mut i = 0; + loop { + nodes.push((i, source(x_at(i), y))); + if i + 1 >= n { + break; + } + i = (i + STEP).min(n - 1); + } + for pair in nodes.windows(2) { + let [(a, pa), (b, pb)] = pair else { continue }; + let linear = match (pa, pb) { + (Some(pa), Some(pb)) if b - a >= 2 => { + let mid = a + (b - a) / 2; + let t = (mid - a) as f64 / (b - a) as f64; + let guess = (pa.0 + (pb.0 - pa.0) * t, pa.1 + (pb.1 - pa.1) * t); + source(x_at(mid), y).is_some_and(|m| (m.0 - guess.0).abs().max((m.1 - guess.1).abs()) <= TOLERANCE) + } + (Some(_), Some(_)) => true, + _ => false, + }; + for k in *a..*b { + let t = (k - a) as f64 / (b - a) as f64; + let c = match (linear, pa, pb) { + (true, Some(pa), Some(pb)) => Some((pa.0 + (pb.0 - pa.0) * t, pa.1 + (pb.1 - pa.1) * t)), + _ => source(x_at(k), y), + }; + if let Some(slot) = out.get_mut(k) { + *slot = c; + } + } + } + if let Some((last, p)) = nodes.last() + && let Some(c) = out.get_mut(*last) + { + *c = *p; + } +} + +/// `r` on the grid `g`, each cell sampled at its centre's image position (`source`). +fn sampled( + r: &Raster, + g: &Grid, + canvas: Canvas, + method: Resampling, + source: &dyn Fn(f64, f64) -> Option<(f64, f64)>, + name: String, +) -> Result { + let (w, h) = (f64::from(r.width), f64::from(r.height)); + let (ow, oh) = (g.width, g.height); + let method = if matches!(r.color, ColorHint::Palette(_)) { Resampling::Nearest } else { method }; + let sampler = Sampler { r, method, alpha: canvas.alpha_src, round: !r.data_type.starts_with("Float") }; + let mut bands: Vec> = vec![vec![canvas.outside; ow * oh]; canvas.out_bands]; + let mut coords: Vec> = vec![None; ow]; + let mut nodes = Vec::new(); + let mut px = vec![f32::NAN; canvas.nb]; + let x_at = |ox: usize| g.min_x + (ox as f64 + 0.5) * g.cell_x; + for oy in 0..oh { + let y = g.max_y - (oy as f64 + 0.5) * g.cell_y; + row_positions(&x_at, y, source, &mut nodes, &mut coords); for (ox, c) in coords.iter().enumerate() { let Some((c, row)) = *c else { continue }; if !(c >= 0.0 && row >= 0.0 && c < w && row < h) { continue; } sampler.sample(c, row, &mut px); - let at = oy * ow + ox; - for (band, v) in bands.iter_mut().zip(&px) { - if let Some(o) = band.get_mut(at) { - // A cell with no value is no-data (or see-through) like one outside. - *o = if v.is_nan() { if clear { 0.0 } else { outside } } else { *v }; + canvas.set(&mut bands, oy * ow + ox, &px); + } + } + canvas.finish(r, g, bands, name) +} + +/// `r` on the grid `g`, each cell the mean (or, with `mode`, the most frequent value, the least +/// of those tied) of the pixels whose centres fall in it. A cell's corners are placed on the +/// image (`source`) and the pixel centres inside the quadrilateral they make are taken; one with +/// none takes the pixel under its centre. A palette image takes the mode. +fn aggregated( + r: &Raster, + g: &Grid, + canvas: Canvas, + mode: bool, + source: &dyn Fn(f64, f64) -> Option<(f64, f64)>, + name: String, +) -> Result { + let (w, h) = (r.width as usize, r.height as usize); + let (ow, oh) = (g.width, g.height); + let mode = mode || matches!(r.color, ColorHint::Palette(_)); + let round = !r.data_type.starts_with("Float"); + let nearest = Sampler { r, method: Resampling::Nearest, alpha: None, round }; + let mut bands: Vec> = vec![vec![canvas.outside; ow * oh]; canvas.out_bands]; + let mut nodes = Vec::new(); + // Cell corners: a row's top edge, then its bottom edge (the next row's top). + let x_at = |i: usize| g.min_x + i as f64 * g.cell_x; + let mut top: Vec> = vec![None; ow + 1]; + let mut bottom: Vec> = vec![None; ow + 1]; + row_positions(&x_at, g.max_y, source, &mut nodes, &mut top); + let mut values: Vec> = vec![Vec::new(); canvas.nb]; + let mut cell = vec![f32::NAN; canvas.nb]; + for oy in 0..oh { + row_positions(&x_at, g.max_y - (oy + 1) as f64 * g.cell_y, source, &mut nodes, &mut bottom); + for ox in 0..ow { + let corners = (top.get(ox), top.get(ox + 1), bottom.get(ox + 1), bottom.get(ox)); + let (Some(Some(tl)), Some(Some(tr)), Some(Some(br)), Some(Some(bl))) = corners else { continue }; + let quad = [*tl, *tr, *br, *bl]; + let (c0, c1) = quad.iter().fold((f64::MAX, f64::MIN), |(lo, hi), p| (lo.min(p.0), hi.max(p.0))); + let (r0, r1) = quad.iter().fold((f64::MAX, f64::MIN), |(lo, hi), p| (lo.min(p.1), hi.max(p.1))); + if c1 <= 0.0 || r1 <= 0.0 || c0 >= w as f64 || r0 >= h as f64 { + continue; + } + // The pixels whose centres (at half integers) may be inside. + let first = |lo: f64| (lo - 0.5).ceil().max(0.0) as usize; + let past = |hi: f64, n: usize| ((hi - 0.5).floor() + 1.0).clamp(0.0, n as f64) as usize; + let (i0, i1, j0, j1) = (first(c0), past(c1, w), first(r0), past(r1, h)); + values.iter_mut().for_each(Vec::clear); + let mut inside = false; + for j in j0..j1 { + for i in i0..i1 { + if !holds(&quad, (i as f64 + 0.5, j as f64 + 0.5)) { + continue; + } + inside = true; + let idx = j * w + i; + for (b, v) in values.iter_mut().enumerate() { + if let Some(x) = r.bands.get(b).and_then(|band| band.data.get(idx)).filter(|x| r.is_valid(**x)) { + v.push(*x); + } + } } } - if add_alpha && let Some(o) = bands.get_mut(nb).and_then(|b| b.get_mut(at)) { - *o = opaque; + let at = oy * ow + ox; + if inside { + for (o, v) in cell.iter_mut().zip(values.iter_mut()) { + *o = if mode { mode_of(v) } else { mean_of(v, round) }; + } + } else { + // A grid finer than the raster: the pixel under the cell's centre. + let (cx, cy) = (quad.iter().map(|p| p.0).sum::() / 4.0, quad.iter().map(|p| p.1).sum::() / 4.0); + if !(cx >= 0.0 && cy >= 0.0 && cx < w as f64 && cy < h as f64) { + continue; + } + nearest.sample(cx, cy, &mut cell); } + canvas.set(&mut bands, at, &cell); } + std::mem::swap(&mut top, &mut bottom); + } + canvas.finish(r, g, bands, name) +} + +/// Whether `p` is inside the quadrilateral `q`: an output cell's top left, top right, bottom +/// right and bottom left corners as image positions. Its top and left edges count and its right +/// and bottom ones don't, so a pixel centre on the edge between two cells is in one of them. +fn holds(q: &[(f64, f64); 4], p: (f64, f64)) -> bool { + let [a, b, c, d] = *q; + // Twice the signed area: which side of each edge is inside. + let turn = (a.0 * b.1 - b.0 * a.1) + (b.0 * c.1 - c.0 * b.1) + (c.0 * d.1 - d.0 * c.1) + (d.0 * a.1 - a.0 * d.1); + if turn == 0.0 || !turn.is_finite() { + return false; } + let s = turn.signum(); + let side = |from: (f64, f64), to: (f64, f64)| s * ((to.0 - from.0) * (p.1 - from.1) - (to.1 - from.1) * (p.0 - from.0)); + side(a, b) >= 0.0 && side(b, c) > 0.0 && side(c, d) > 0.0 && side(d, a) >= 0.0 +} - let color = match (add_alpha, rgb) { - (true, true) => ColorHint::Rgb { alpha: true }, - (true, false) => ColorHint::GrayAlpha, - _ => r.color.clone(), - }; - let transform = GeoTransform { origin_x: min_x, origin_y: max_y, pixel_width: cell, pixel_height: -cell }; - let mut out = Raster::new(o.name.clone(), ow as u32, oh as u32, bands, transform, o.crs, r.nodata, color, r.data_type.clone()) - .map_err(|e| e.to_string())?; - for (i, band) in out.bands.iter_mut().enumerate() { - band.name = r.bands.get(i).map_or_else(|| "Alpha".to_string(), |b| b.name.clone()); +/// The mean of `v`, rounded for a whole-number pixel type; NaN when it is empty. +fn mean_of(v: &[f32], round: bool) -> f32 { + if v.is_empty() { + return f32::NAN; } - Ok(out) + let m = (v.iter().map(|x| f64::from(*x)).sum::() / v.len() as f64) as f32; + if round { m.round() } else { m } +} + +/// The most frequent value of `v`, the least of those tied; NaN when it is empty. Sorts `v`. +fn mode_of(v: &mut [f32]) -> f32 { + v.sort_unstable_by(f32::total_cmp); + let (mut best, mut best_n) = (f32::NAN, 0usize); + let mut k = 0; + while let Some(x) = v.get(k).copied() { + let n = v.get(k..).map_or(0, |rest| rest.iter().take_while(|y| **y == x).count()).max(1); + // Sorted ascending: a later value only wins with more. + if n > best_n { + (best, best_n) = (x, n); + } + k += n; + } + best } #[cfg(test)] @@ -562,4 +907,123 @@ mod tests { let out = warp(&r, &fit, &opts(None, Resampling::Bilinear)).unwrap(); assert_eq!(out.value(3, c, rw), Some(255.0)); } + + /// A `w` × `h` raster in UTM 13N with 1 m pixels, its top left corner at (0, `h`), no-data -1. + fn utm(w: u32, h: u32, values: Vec, data_type: &str) -> Raster { + let gt = GeoTransform { origin_x: 0.0, origin_y: f64::from(h), pixel_width: 1.0, pixel_height: -1.0 }; + let crs = Crs::from_epsg(32613).unwrap(); + Raster::new("r", w, h, vec![values], gt, crs, Some(-1.0), ColorHint::Gray, data_type).unwrap() + } + + fn values(r: &Raster) -> Vec> { + (0..i64::from(r.height)).flat_map(|row| (0..i64::from(r.width)).map(move |c| (c, row))).map(|(c, row)| r.value(0, c, row)).collect() + } + + /// Average and mode take every pixel whose centre falls in a cell, no-data skipped; mode's + /// ties go to the least value. + #[test] + fn average_and_mode_take_every_pixel_in_a_cell() { + // 4 × 4 pixels, value = column + 10 · row; pixel (1, 0) is no-data. + let mut v: Vec = (0..16).map(|i| (i % 4 + 10 * (i / 4)) as f32).collect(); + v[1] = -1.0; + let r = utm(4, 4, v, "Float32"); + let g = Grid::snapped(r.crs, [0.0, 0.0, 4.0, 4.0], (2.0, 2.0), (0.0, 0.0)).unwrap(); + assert_eq!((g.width, g.height), (2, 2)); + let avg = project(&r, &g, Method::Average, "a".into()).unwrap(); + // Top left: 0, 10 and 11 (1 is no-data); top right: 2, 3, 12, 13. + assert_eq!(values(&avg), vec![Some(7.0), Some(7.5), Some(25.5), Some(27.5)]); + // Classes: a cell of 3, 5, 5, 3 ties and takes 3; one of 7, 7, 7, 2 takes 7. + let classes = vec![3.0, 5.0, 7.0, 7.0, 5.0, 3.0, 7.0, 2.0, 1.0, 1.0, 4.0, 4.0, 1.0, 9.0, 4.0, 4.0]; + let r = utm(4, 4, classes, "UInt8"); + let mode = project(&r, &g, Method::Mode, "m".into()).unwrap(); + assert_eq!(values(&mode), vec![Some(3.0), Some(7.0), Some(1.0), Some(4.0)]); + // A whole-number type's mean is rounded: (3 + 5 + 5 + 3) / 4 = 4, (7 + 7 + 7 + 2) / 4 = 5.75. + let avg = project(&r, &g, Method::Average, "a".into()).unwrap(); + assert_eq!(values(&avg)[..2], [Some(4.0), Some(6.0)]); + } + + /// A pixel centre on the edge between cells is in exactly one of them, and a grid finer than + /// the raster takes the pixel under each cell's centre. + #[test] + fn each_pixel_is_in_one_cell_and_a_finer_grid_takes_the_nearest() { + // Value = column + 10 · row. + let r = utm(4, 4, (0..16).map(|i| (i % 4 + 10 * (i / 4)) as f32).collect(), "Float32"); + // Cells of 2 whose edges run through pixel centres: a centre on a cell's top or left edge + // is in it, one on its right or bottom edge in the next, so each is taken once. + let g = Grid::snapped(r.crs, [0.0, 0.0, 4.0, 4.0], (2.0, 2.0), (0.5, 0.5)).unwrap(); + assert_eq!((g.width, g.height, g.min_x, g.max_y), (3, 3, -1.5, 4.5)); + let out = project(&r, &g, Method::Average, "a".into()).unwrap(); + let expected = [None, Some(0.5), Some(2.5), None, Some(15.5), Some(17.5), None, Some(30.5), Some(32.5)]; + assert_eq!(values(&out), expected); + // Cells of 0.5: each takes the pixel it is in. + let g = Grid::snapped(r.crs, [0.0, 0.0, 4.0, 4.0], (0.5, 0.5), (0.0, 0.0)).unwrap(); + let out = project(&r, &g, Method::Average, "a".into()).unwrap(); + assert_eq!((out.width, out.height), (8, 8)); + assert_eq!((out.value(0, 0, 0), out.value(0, 3, 1), out.value(0, 7, 7)), (Some(0.0), Some(1.0), Some(33.0))); + } + + /// Projected data keeps its bands: a grey byte raster with no no-data value gains no alpha + /// band (a scan being georeferenced does), and its cells off the raster are 0, as GDAL leaves + /// them; a floating-point one's are NaN, which is no value. + #[test] + fn projected_data_keeps_its_bands() { + let gt = GeoTransform { origin_x: 0.0, origin_y: 4.0, pixel_width: 1.0, pixel_height: -1.0 }; + let crs = Crs::from_epsg(32613).unwrap(); + let g = Grid::snapped(crs, [-2.0, 0.0, 4.0, 4.0], (1.0, 1.0), (0.0, 0.0)).unwrap(); + let classes = Raster::new("c", 4, 4, vec![vec![3.0; 16]], gt, crs, None, ColorHint::Gray, "UInt8").unwrap(); + let out = project(&classes, &g, Method::Mode, "m".into()).unwrap(); + assert_eq!((out.band_count(), out.color.clone(), out.nodata), (1, ColorHint::Gray, None)); + assert_eq!((out.value(0, 0, 0), out.value(0, 2, 0)), (Some(0.0), Some(3.0))); + let floats = Raster::new("f", 4, 4, vec![vec![3.5; 16]], gt, crs, None, ColorHint::Gray, "Float32").unwrap(); + let out = project(&floats, &g, Method::Sample(Resampling::Bilinear), "b".into()).unwrap(); + assert_eq!((out.band_count(), out.value(0, 0, 0), out.value(0, 3, 1)), (1, None, Some(3.5))); + } + + /// Grids with the same cell and origin line up wherever their areas are; a grid too large to + /// hold is refused before anything is made. + #[test] + fn snapped_grids_line_up_and_are_bounded() { + let crs = Crs::from_epsg(31983).unwrap(); + let a = Grid::snapped(crs, [5_123_456.7, 9_000_001.0, 5_200_000.0, 9_100_000.0], (100.0, 100.0), (0.0, 0.0)).unwrap(); + let b = Grid::snapped(crs, [5_150_050.0, 9_050_049.0, 5_300_000.0, 9_200_000.0], (100.0, 100.0), (0.0, 0.0)).unwrap(); + for g in [a, b] { + assert_eq!((g.min_x % 100.0, g.max_y % 100.0), (0.0, 0.0)); + } + assert!(a.min_x <= 5_123_456.7 && a.min_x + a.width as f64 * 100.0 >= 5_200_000.0); + let e = Grid::snapped(crs, [0.0, 0.0, 1e6, 1e6], (0.01, 0.01), (0.0, 0.0)).unwrap_err(); + assert!(e.contains("give a larger cellSize"), "{e}"); + for bad in [(0.0, 1.0), (-1.0, 1.0), (f64::NAN, 1.0), (1.0, f64::INFINITY)] { + assert!(Grid::snapped(crs, [0.0, 0.0, 10.0, 10.0], bad, (0.0, 0.0)).is_err(), "{bad:?}"); + } + assert!(Grid::snapped(crs, [10.0, 0.0, 0.0, 10.0], (1.0, 1.0), (0.0, 0.0)).is_err()); + assert!(Grid::snapped(crs, [0.0, 0.0, 10.0, 10.0], (1.0, 1.0), (f64::NAN, 0.0)).is_err()); + } + + /// A raster projected to another system and back keeps its values in place, within a cell. + #[test] + fn a_projection_and_back_keeps_values_in_place() { + // 40 × 40 pixels of 0.001° near Teresina, value = column + 100 · row. + let gt = GeoTransform { origin_x: -42.82, origin_y: -5.08, pixel_width: 0.001, pixel_height: -0.001 }; + let v: Vec = (0..1600).map(|i| (i % 40 + 100 * (i / 40)) as f32).collect(); + let r = Raster::new("t", 40, 40, vec![v], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap(); + let utm23s = Crs::from_epsg(31983).unwrap(); + let (bounds, cell) = footprint_in(&r, utm23s).unwrap(); + assert!((cell - 110.0).abs() < 2.0, "{cell}"); + let g = Grid::snapped(utm23s, bounds, (cell, cell), (0.0, 0.0)).unwrap(); + let there = project(&r, &g, Method::Sample(Resampling::Nearest), "p".into()).unwrap(); + let back = project(&there, &Grid::of(&r), Method::Sample(Resampling::Nearest), "b".into()).unwrap(); + let (mut same, mut near) = (0, 0); + for (a, b) in values(&r).iter().zip(values(&back)) { + let (Some(a), Some(b)) = (a, b) else { continue }; + // Within a pixel each way: a column of 1, a row of 100. + let (dc, dr) = ((a % 100.0 - b % 100.0).abs(), ((a / 100.0).floor() - (b / 100.0).floor()).abs()); + same += usize::from(a == &b); + near += usize::from(dc <= 1.0 && dr <= 1.0); + } + assert_eq!(near, 1600, "every pixel within a pixel of its place"); + assert!(same > 1000, "{same} of 1600 in the very same pixel"); + // A longitude/latitude box placed in the projection covers its corners. + let b = bounds_in([-42.82, -5.12, -42.78, -5.08], utm23s).unwrap(); + assert!((b[0] - bounds[0]).abs() < 1.0 && (b[3] - bounds[3]).abs() < 1.0, "{b:?} {bounds:?}"); + } } diff --git a/crates/assistant/src/tests.rs b/crates/assistant/src/tests.rs index 31a41b55f..a315b89ec 100644 --- a/crates/assistant/src/tests.rs +++ b/crates/assistant/src/tests.rs @@ -284,7 +284,7 @@ fn gate_rules() { assert!(g("list_layers", json!({})).is_none()); assert!(g("run_command", json!({"command": "info.fileLayers", "params": {"path": "/etc/passwd"}})).is_some()); assert!(g("run_command", json!({"command": "info.fileColumns", "params": {"path": "/etc/passwd"}})).is_some()); - // Georeferencing is a file command only when it writes a file. + // Georeferencing and projecting are file commands only when they write a file. let extent = json!({"layer": "r", "extent": [10.0, 45.0, 10.04, 45.04]}); assert!(g("run_command", json!({"command": "raster.georeference", "params": extent})).is_none()); let pts = json!({"layer": "r", "points": [[0, 0, 1, 1], [1, 0, 2, 1], [0, 1, 1, 2]], "transform": "affine", "crs": 32613}); @@ -292,6 +292,12 @@ fn gate_rules() { let mut file = pts.clone(); file["path"] = json!("/tmp/x.tif"); assert!(g("run_command", json!({"command": "raster.georeference", "params": file})).is_some()); + // So is projecting a raster. + let project = json!({"layer": "r", "crs": 31983, "cellSize": 100, "extent": "Countries"}); + assert!(g("run_command", json!({"command": "raster.project", "params": project})).is_none()); + let mut file = project.clone(); + file["path"] = json!("/tmp/x.tif"); + assert!(g("run_command", json!({"command": "raster.project", "params": file})).is_some()); } /// What a command says beside its result reaches the model (issue #85). diff --git a/crates/engine/src/cmd/imagery.rs b/crates/engine/src/cmd/imagery.rs index b9493063d..3f442b7fa 100644 --- a/crates/engine/src/cmd/imagery.rs +++ b/crates/engine/src/cmd/imagery.rs @@ -41,6 +41,9 @@ pub fn commands() -> Vec { C::new("raster.georeference", "Georeference", "Analyze › Alignment", georeference).params( r#"{"layer": raster, "points"?: [[column, row, lon, lat], …] (two or more; with "transform": [column, row, x, y, on?]), "extent"?: [minLon, minLat, maxLon, maxLat], "transform"?: "similarity"|"affine"|"projective"|"polynomial2"|"polynomial3"|"spline" (fits the points and warps the raster onto a north-up grid covering the warped image; without it the points and extent only shift and scale the grid), "pointsCrs"?: EPSG code or "EPSG:n" of the points' x and y (default 4326), "crs"?: the output's (default pointsCrs), "resampling"?: "nearest"|"bilinear"|"cubic" (default bilinear), "cellSize"?: output units (default about the source pixel's size on the ground), "name"?: new layer for the result (default replaces the layer's raster), "path"?: "*.tif" (also writes the warp as a cloud-optimised GeoTIFF)}"#, ), + C::new("raster.project", "Project Raster", "Analyze › Alignment", project_raster).params( + r#"{"layer": raster, "crs"?: EPSG code or "EPSG:n" of the output (default the layer's), "cellSize"?: number or [x, y] in the output's units (default about the layer's pixel there), "origin"?: [x, y] a cell corner falls on (default [0, 0], so grids of the same cellSize line up), "extent"?: polygon layer (its selected or all features' bounds) | [west, south, east, north] in degrees | "view" (default the layer's footprint), "grid"?: raster layer whose coordinate system, cells and extent are taken exactly (instead of crs, cellSize, origin and extent), "resampling"?: "nearest" (default) | "bilinear" | "cubic" | "average" (the mean of the pixels whose centres fall in a cell) | "mode" (their most frequent value, the least of those tied), "name"?: new layer (default " projected"), "path"?: GeoTIFF to write it to}"#, + ), // A fit is a measurement of the points, not a project change: no undo step. C::new("raster.gcpFit", "Control Point Fit", "Analyze › Alignment", gcp_fit).view().params( r#"{"points": [[column, row, x, y, on?], …] (on: 1/0 or true/false, default on), "transform"?: "similarity"|"affine"|"projective"|"polynomial2"|"polynomial3"|"spline" (default: similarity under 3 points, else affine), "compare"?: bool (fit every transform and recommend one; transform and leaveOneOut are then ignored), "leaveOneOut"?: bool (default true)}"#, @@ -321,6 +324,114 @@ fn warp_to_fit(s: &mut Session, p: &Value, layer: &str, r: &Raster) -> CmdResult Ok(info) } +/// `raster.project`: the layer on a grid of another coordinate system, cell size, origin and +/// extent, or on another raster's grid, each cell sampled at its centre or taken from every +/// pixel whose centre falls in it (see `rustgis_analysis::warp::project`). A new layer. +fn project_raster(s: &mut Session, p: &Value) -> CmdResult { + let r = raster(s, p)?; + r.require_coordinates()?; + let method = match p.get("resampling") { + None | Some(Value::Null) => warp::Method::Sample(warp::Resampling::Nearest), + Some(v) => v.as_str().and_then(warp::Method::parse).ok_or("resampling is nearest, bilinear, cubic, average or mode")?, + }; + let grid = match ps(p, "grid") { + Some(template) => { + if ["crs", "cellSize", "origin", "extent"].iter().any(|k| p.get(*k).is_some_and(|v| !v.is_null())) { + return Err( + "grid takes that raster's coordinate system, cells and extent: leave out crs, cellSize, origin and extent".into() + ); + } + let t = raster_named(s, template)?; + t.require_coordinates()?; + warp::Grid::of(&t) + } + None => { + let crs = crs_param(p, "crs")?.unwrap_or(r.crs); + let (footprint, pixel) = warp::footprint_in(&r, crs)?; + let positive = |v: &Value| v.as_f64().filter(|c| c.is_finite() && *c > 0.0); + let bad_cell = "cellSize is a positive number, or [x, y], in the output's units"; + let cell = match p.get("cellSize") { + None | Some(Value::Null) => (pixel, pixel), + Some(Value::Array(a)) => match a.as_slice() { + [x, y] => (positive(x).ok_or(bad_cell)?, positive(y).ok_or(bad_cell)?), + _ => return Err(bad_cell.into()), + }, + Some(v) => positive(v).map(|c| (c, c)).ok_or(bad_cell)?, + }; + let origin = match p.get("origin") { + None | Some(Value::Null) => (0.0, 0.0), + Some(v) => { + let n: Vec = v.as_array().map(|a| a.iter().filter_map(Value::as_f64).collect()).unwrap_or_default(); + match n.as_slice() { + [x, y] if x.is_finite() && y.is_finite() => (*x, *y), + _ => return Err("origin is [x, y] in the output's units".into()), + } + } + }; + let bounds = match p.get("extent") { + None | Some(Value::Null) => footprint, + Some(Value::String(l)) if l != "view" => { + let layer = s.map().layer_by_name(l).filter(|l| l.raster.is_none()).ok_or_else(|| format!("no vector layer '{l}'"))?; + let rows: Vec = + if layer.selection.is_empty() { (0..layer.dataset.features.len()).collect() } else { layer.selection.clone() }; + let e = rows + .iter() + .filter_map(|row| layer.dataset.features.get(*row).and_then(|f| f.geometry.as_ref())) + .fold(Extent::EMPTY, |e, g| e.union(&g.extent())); + if e.is_empty() { + return Err(format!("{l} has no features to take an extent from")); + } + warp::bounds_in([e.min_x, e.min_y, e.max_x, e.max_y], crs)? + } + Some(v) => { + let area = super::catalogs::area(s, &json!({"bbox": v}))?; + if area[0] > area[2] { + return Err("an extent across the antimeridian: give it one side at a time".into()); + } + warp::bounds_in(area, crs)? + } + }; + warp::Grid::snapped(crs, bounds, cell, origin)? + } + }; + let path = match p.get("path") { + None | Some(Value::Null) => None, + Some(v) => Some(v.as_str().filter(|t| is_tif(t)).ok_or("path is a *.tif or *.tiff file")?.to_string()), + }; + let out_name = name(p, format!("{} projected", r.name)); + let out = warp::project(&r, &grid, method, out_name.clone())?; + let mut info = json!({ + "name": out_name, "crs": grid.crs.label(), "resampling": method.id(), "width": out.width, "height": out.height, "bands": out.band_count(), + "cellSize": [grid.cell_x, grid.cell_y], "corner": [grid.min_x, grid.max_y], + }); + if let Some((fw, fh)) = r.overview { + let note = + format!("{}: projected at reduced resolution, from its {} × {} overview of the {fw} × {fh} raster", r.name, r.width, r.height); + s.notify(note, false); + info["projectedFromOverview"] = json!({"width": r.width, "height": r.height, "of": [fw, fh]}); + } + let e = out.extent(); + info["extent"] = json!([e.min_x, e.min_y, e.max_x, e.max_y]); + // A file written is the layer's source; without one the pixels are kept in the project. + let rl = match &path { + Some(path) => { + let cog = + rustgis_data::raster::cogwrite::write_file(&out, 0, std::path::Path::new(path)).map_err(|e| format!("{path}: {e}"))?; + info["path"] = json!(path); + info["cog"] = json!({"overviews": cog.levels.len().saturating_sub(1), "bigTiff": cog.big, "bytes": cog.bytes}); + RasterLayer::from_path(out, path.clone()) + } + None => layer_of(out)?, + }; + let v = super::raster::add_raster(s, &out_name, rl)?; + if let (Some(o), Some(i)) = (info.as_object_mut(), v.as_object()) { + for (k, val) in i { + o.entry(k.clone()).or_insert_with(|| val.clone()); + } + } + Ok(info) +} + fn is_tif(path: &str) -> bool { std::path::Path::new(path) .extension() diff --git a/crates/engine/src/cmd/plugins.rs b/crates/engine/src/cmd/plugins.rs index 079c68cb0..77273a727 100644 --- a/crates/engine/src/cmd/plugins.rs +++ b/crates/engine/src/cmd/plugins.rs @@ -355,8 +355,9 @@ pub fn reaches(spec: &C, params: &Value) -> Option { if DATABASES.contains(&spec.id) { return Some(Reach::Files); } - // Georeferencing writes a file only when it is given a `path`; the rest of it stays in the project. - if spec.id == "raster.georeference" { + // Georeferencing and projecting write a file only when given a `path`; the rest of them stays + // in the project. + if spec.id == "raster.georeference" || spec.id == "raster.project" { return super::items::names_place(params, PLACE_KEYS, 0).then_some(Reach::Files); } let files = super::items::names_place(params, PLACE_KEYS, 0) diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index 2a2e386de..f5df479b4 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -16757,8 +16757,9 @@ fn scans_have_no_coordinates_to_give_away() { let e = s.run("share.exportRaster", &json!({"path": out.display().to_string(), "cog": true})).unwrap_err(); assert!(e.contains("georeference it first") && !out.exists(), "as a COG: {e}"); let before = s.map().layer_by_name("flat").unwrap().raster.clone().unwrap().placement; - let tools: [(&str, Value); 14] = [ + let tools: [(&str, Value); 15] = [ ("raster.function", json!({"layer": "flat", "function": "hillshade"})), + ("raster.project", json!({"layer": "flat", "crs": 32613})), ("raster.editPixels", json!({"layer": "flat", "polygon": [[0.0, 0.0], [0.1, 0.0], [0.0, 0.1]], "value": 1})), ("raster.process", json!({"layer": "flat", "operation": "resample", "factor": 2})), ("raster.process", json!({"layer": "flat", "operation": "clip"})), @@ -17618,3 +17619,84 @@ fn calculate_geometry_updates_ellipsoidal_lengths_in_one_undo_step() { let out = &s.map().layer_by_name("routes").unwrap().dataset; assert!((out.value(0, length).as_f64().unwrap() - 111.31949079327356).abs() < 1e-6); } + +/// `raster.project` puts a raster on a grid of another coordinate system, cell size, origin and +/// extent, or on another raster's grid: grids of the same cell and origin line up, a polygon +/// layer's selection gives the extent, and bad parameters are refused before anything is made. +#[test] +fn a_raster_projected_onto_a_chosen_grid() { + use rustgis_data::raster::{ColorHint, Crs, GeoTransform, Raster}; + let mut s = world(); + // 40 × 40 pixels of 0.001° near Teresina, value = column + 100 · row. + let gt = GeoTransform { origin_x: -42.82, origin_y: -5.08, pixel_width: 0.001, pixel_height: -0.001 }; + let v: Vec = (0..1600).map(|i| (i % 40 + 100 * (i / 40)) as f32).collect(); + let r = Raster::new("the", 40, 40, vec![v], gt, Crs::LonLat, None, ColorHint::Gray, "Float32").unwrap(); + let bytes = rustgis_data::raster::geotiff::write(&r, 0).unwrap(); + crate::cmd::raster::add_raster(&mut s, "the", rustgis_map::RasterLayer::from_bytes(r, "the.tif", bytes)).unwrap(); + + // SIRGAS 2000 / UTM 23S, 200 m cells averaged: about 22 × 22 cells, means within the values. + let v = s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 200, "resampling": "average"})).unwrap(); + assert_eq!( + (v["crs"].as_str(), v["resampling"].as_str(), v["name"].as_str()), + (Some("EPSG:31983"), Some("average"), Some("the projected")) + ); + let out = out_raster(&s, &v); + assert!((22..=24).contains(&out.width) && (22..=24).contains(&out.height), "{} × {}", out.width, out.height); + let mid = out.value(0, i64::from(out.width) / 2, i64::from(out.height) / 2).unwrap(); + assert!((0.0..=3939.0).contains(&mid) && mid.fract() != 0.0, "a mean of several pixels: {mid}"); + // Cell corners on whole multiples of the cell (origin [0, 0]), here and over a part of it. + let on_grid = |v: &Value| { + let c = &v["corner"]; + (c[0].as_f64().unwrap() / 200.0).fract() == 0.0 && (c[1].as_f64().unwrap() / 200.0).fract() == 0.0 + }; + assert!(on_grid(&v)); + let part = s + .run( + "raster.project", + &json!({"layer": "the", "crs": "EPSG:31983", "cellSize": 200, "extent": [-42.81, -5.11, -42.79, -5.09], "name": "part"}), + ) + .unwrap(); + assert!(on_grid(&part) && part["width"].as_u64() < v["width"].as_u64(), "{part}"); + // Another origin moves the corners with it. + let moved = + s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 200, "origin": [50, 50], "name": "moved"})).unwrap(); + assert_eq!((moved["corner"][0].as_f64().unwrap() - 50.0).rem_euclid(200.0), 0.0, "{moved}"); + + // Onto that grid exactly, by nearest neighbour (the default). + let again = s.run("raster.project", &json!({"layer": "the", "grid": "the projected", "name": "again"})).unwrap(); + let a = out_raster(&s, &again); + assert_eq!((a.width, a.height, a.transform, a.crs), (out.width, out.height, out.transform, out.crs)); + assert_eq!(again["resampling"], "nearest"); + + // The extent of a polygon layer's selection: Brazil, in 50 km cells. + s.run("select.byAttributes", &json!({"layer": "Countries", "where": "name = 'Brazil'"})).unwrap(); + let br = + s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 50_000, "extent": "Countries", "name": "br"})).unwrap(); + let (w, h) = (br["width"].as_u64().unwrap(), br["height"].as_u64().unwrap()); + assert!((70..=130).contains(&w) && (70..=130).contains(&h), "{w} × {h}"); + + // Refused: a grid with its own settings, bad cells, origins, extents and methods, a grid too + // large to hold, and a file that isn't a GeoTIFF. + let bad = [ + json!({"layer": "the", "grid": "the projected", "crs": 31983}), + json!({"layer": "the", "grid": "Countries"}), + json!({"layer": "the", "cellSize": -1}), + json!({"layer": "the", "cellSize": [1, 0]}), + json!({"layer": "the", "cellSize": "x"}), + json!({"layer": "the", "origin": "x"}), + json!({"layer": "the", "origin": [1, f64::NAN]}), + json!({"layer": "the", "extent": [170.0, -10.0, -170.0, 10.0]}), + json!({"layer": "the", "extent": "Nowhere"}), + json!({"layer": "the", "resampling": "lanczos"}), + json!({"layer": "the", "crs": 31983, "cellSize": 0.001}), + json!({"layer": "the", "crs": "nonsense"}), + json!({"layer": "the", "path": "out.png"}), + ]; + let before = s.map().layers.len(); + for p in bad { + assert!(s.run("raster.project", &p).is_err(), "{p}"); + } + assert_eq!(s.map().layers.len(), before, "nothing added by a refused run"); + let e = s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 0.001})).unwrap_err(); + assert!(e.contains("give a larger cellSize"), "{e}"); +} diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index 9e24ca0c9..e01780c4e 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -325,7 +325,8 @@ honour it, this restricts an unselected main feature input by bounding-box inter clipping its geometry; an existing selection takes precedence and overlay inputs stay whole. The Geostatistical Wizard uses the environment's extent and raster cell size (degrees). IDW Interpolation and Ordinary Kriging instead use their own `extent` and `pixelSize` parameters; -the environment does not set their output grids. Raster Functions › Resample matches a target +the environment does not set their output grids. Project Raster (`raster.project`) puts a raster +on a grid of any coordinate system, cell size and extent; Raster Functions › Resample matches a target raster's grid. The **Script** pane runs commands from text, one per line (`id {json}`); its box grows with the script down to the bottom of the pane, @@ -462,6 +463,36 @@ The Analyze tab's imagery groups open their tools in the Workflows pane: coordinate systems or cell sizes, or opened from an overview, are refused rather than quietly resampled, as are tiles whose cells don't line up; the mosaic has to fit in memory (120 million values: cells times bands). +- **Project Raster** (`raster.project`) puts a raster on a grid of your choosing, so rasters from + different sources line up cell for cell (for map algebra, the Suitability Modeler, a series for + Trend Analysis, or the inputs of a model): `{"layer": "ndvi", "crs": 31983, "cellSize": 100, + "extent": "Study area", "resampling": "average", "path": "ndvi_100m.tif"}`. It has no pane yet: + run it from the Script pane, the AI Assistant, the command line or Python. The result is a new + layer (`name`, by default the layer's name with "projected"), and with `path` a cloud-optimised + GeoTIFF too. + - **The grid:** `crs` (by default the layer's), `cellSize` in its units (one number, or + `[x, y]`; by default about the layer's pixel there), `origin`, a point a cell corner falls on + (`[0, 0]` by default, so cell corners are whole multiples of the cell size and any two rasters + projected with the same `crs` and `cellSize` line up), and `extent`: a polygon layer (its + selected features' bounds, or all of them), `[west, south, east, north]` in degrees or + `"view"`, by default the layer's footprint. The extent is widened to whole cells. + - **Another raster's grid:** `"grid": "dem_100m"` takes its coordinate system, cells and extent + exactly, instead of `crs`, `cellSize`, `origin` and `extent`. + - **Resampling:** `nearest` (the default: it never invents a value, which is right for classes), + `bilinear` or `cubic` sample each cell at its centre. `average` and `mode` take every pixel whose + centre falls in the cell (no-data skipped), its mean or its most frequent value (the least of + those tied), which is what a coarser grid of elevation or land-cover classes should hold. A + cell no pixel centre falls in, on a grid finer than the raster, takes the pixel under its + centre. Classes take `nearest` or `mode`, continuous values `bilinear` or `cubic`, or `average` + going coarser. A whole-number pixel type keeps its type, so its mean is rounded. + - Cells off the raster are its no-data value. A raster without one keeps its bands as they are, + with NaN there in a floating-point type and 0 in a whole-number one, as GDAL leaves them; a + colour image gains an alpha band instead, so they are see-through. + - Datums are shifted as everywhere in RustGIS, through each system's transformation to WGS 84 + (a Helmert step for SAD69 or Córrego Alegre, a few metres from the IBGE's grid shifts). + - The result has to fit in memory (120 million values, cells times bands; 32 million in the + web app): a grid past that is refused, naming the limit. A raster read from its overview is + projected at that resolution, and the result says so. - **Classification Wizard** (unsupervised k-means) and **Classification Tools** (minimum distance from training polygons). - **Mensuration**: distances along the surface, heights between two points, areas, and volumes @@ -511,6 +542,10 @@ 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. + + For the mean of a finer layer's pixels in each cell, or for a grid of your own, put the rasters + on one grid first with Project Raster (`"grid"`, the first one's, and `"resampling": + "average"`). - **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' From e720e0325b5472e41635d7d87c5780c8e1b85c75 Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Mon, 26 Oct 2026 12:42:47 -0300 Subject: [PATCH 2/7] Analysis: raster.project's review fixes (#821) An origin with anything but two numbers is refused, not read from the numbers left; a grid whose snapped corners overflow is refused. Average weighs a colour image's colours by alpha, as bilinear does, so see-through pixels take no colour into a cell. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/warp.rs | 46 ++++++++++++++++++++++++++------ crates/engine/src/cmd/imagery.rs | 5 ++-- crates/engine/src/tests.rs | 4 +++ 3 files changed, 45 insertions(+), 10 deletions(-) diff --git a/crates/analysis/src/warp.rs b/crates/analysis/src/warp.rs index 2fffd59f3..eed048b6b 100644 --- a/crates/analysis/src/warp.rs +++ b/crates/analysis/src/warp.rs @@ -341,6 +341,10 @@ impl Grid { let max_x = origin.0 + ((x1 - origin.0) / cx).ceil() * cx; let min_y = origin.1 + ((y0 - origin.1) / cy).floor() * cy; let max_y = origin.1 + ((y1 - origin.1) / cy).ceil() * cy; + // An origin or a cell so far from the area that the snapped corners overflow. + if !([min_x, max_x, min_y, max_y].iter().all(|v| v.is_finite()) && min_x < max_x && min_y < max_y) { + return Err("the origin and cellSize can't make a grid over this area".into()); + } let (w, h) = (((max_x - min_x) / cx).round().max(1.0), ((max_y - min_y) / cy).round().max(1.0)); if w * h > MAX_SAMPLES as f64 { return Err(format!( @@ -634,6 +638,9 @@ fn aggregated( let mut bottom: Vec> = vec![None; ow + 1]; row_positions(&x_at, g.max_y, source, &mut nodes, &mut top); let mut values: Vec> = vec![Vec::new(); canvas.nb]; + // Each value's weight: a colour counts as much as its pixel is opaque, as bilinear weighs it. + let mut weights: Vec> = vec![Vec::new(); canvas.nb]; + let alpha = canvas.alpha_src; let mut cell = vec![f32::NAN; canvas.nb]; for oy in 0..oh { row_positions(&x_at, g.max_y - (oy + 1) as f64 * g.cell_y, source, &mut nodes, &mut bottom); @@ -651,6 +658,7 @@ fn aggregated( let past = |hi: f64, n: usize| ((hi - 0.5).floor() + 1.0).clamp(0.0, n as f64) as usize; let (i0, i1, j0, j1) = (first(c0), past(c1, w), first(r0), past(r1, h)); values.iter_mut().for_each(Vec::clear); + weights.iter_mut().for_each(Vec::clear); let mut inside = false; for j in j0..j1 { for i in i0..i1 { @@ -659,17 +667,19 @@ fn aggregated( } inside = true; let idx = j * w + i; - for (b, v) in values.iter_mut().enumerate() { + let opacity = alpha.and_then(|a| r.bands.get(a)).and_then(|band| band.data.get(idx)).copied().unwrap_or(0.0); + for (b, (v, k)) in values.iter_mut().zip(weights.iter_mut()).enumerate() { if let Some(x) = r.bands.get(b).and_then(|band| band.data.get(idx)).filter(|x| r.is_valid(**x)) { v.push(*x); + k.push(if alpha.is_some_and(|a| a != b) { opacity } else { 1.0 }); } } } } let at = oy * ow + ox; if inside { - for (o, v) in cell.iter_mut().zip(values.iter_mut()) { - *o = if mode { mode_of(v) } else { mean_of(v, round) }; + for ((o, v), k) in cell.iter_mut().zip(values.iter_mut()).zip(&weights) { + *o = if mode { mode_of(v) } else { mean_of(v, k, round) }; } } else { // A grid finer than the raster: the pixel under the cell's centre. @@ -701,12 +711,14 @@ fn holds(q: &[(f64, f64); 4], p: (f64, f64)) -> bool { side(a, b) >= 0.0 && side(b, c) > 0.0 && side(c, d) > 0.0 && side(d, a) >= 0.0 } -/// The mean of `v`, rounded for a whole-number pixel type; NaN when it is empty. -fn mean_of(v: &[f32], round: bool) -> f32 { - if v.is_empty() { +/// The mean of `v`, each value weighed by its `k` (a colour by its pixel's alpha), rounded for a +/// whole-number pixel type; NaN when nothing weighs (empty, or all see-through). +fn mean_of(v: &[f32], k: &[f32], round: bool) -> f32 { + let (sum, weight) = v.iter().zip(k).fold((0.0f64, 0.0f64), |(s, w), (x, k)| (s + f64::from(*x) * f64::from(*k), w + f64::from(*k))); + if weight <= 1e-12 { return f32::NAN; } - let m = (v.iter().map(|x| f64::from(*x)).sum::() / v.len() as f64) as f32; + let m = (sum / weight) as f32; if round { m.round() } else { m } } @@ -716,7 +728,7 @@ fn mode_of(v: &mut [f32]) -> f32 { let (mut best, mut best_n) = (f32::NAN, 0usize); let mut k = 0; while let Some(x) = v.get(k).copied() { - let n = v.get(k..).map_or(0, |rest| rest.iter().take_while(|y| **y == x).count()).max(1); + let n = v.iter().skip(k).take_while(|y| **y == x).count(); // Sorted ascending: a later value only wins with more. if n > best_n { (best, best_n) = (x, n); @@ -942,6 +954,21 @@ mod tests { assert_eq!(values(&avg)[..2], [Some(4.0), Some(6.0)]); } + /// Average weighs a colour image's colours by alpha, as bilinear does: a see-through pixel + /// takes no colour into its cell; the alpha itself is the plain mean. + #[test] + fn average_weighs_colours_by_alpha() { + // 2 × 2 pixels into one cell: opaque red on the left, see-through black on the right. + let gt = GeoTransform { origin_x: 0.0, origin_y: 2.0, pixel_width: 1.0, pixel_height: -1.0 }; + let half = |a: f32, b: f32| vec![a, b, a, b]; + let bands = vec![half(255.0, 0.0), half(0.0, 0.0), half(0.0, 0.0), half(255.0, 0.0)]; + let crs = Crs::from_epsg(32613).unwrap(); + let r = Raster::new("rgba", 2, 2, bands, gt, crs, None, ColorHint::Rgb { alpha: true }, "UInt8").unwrap(); + let g = Grid::snapped(r.crs, [0.0, 0.0, 2.0, 2.0], (2.0, 2.0), (0.0, 0.0)).unwrap(); + let out = project(&r, &g, Method::Average, "a".into()).unwrap(); + assert_eq!([0, 1, 2, 3].map(|b| out.value(b, 0, 0)), [Some(255.0), Some(0.0), Some(0.0), Some(128.0)]); + } + /// A pixel centre on the edge between cells is in exactly one of them, and a grid finer than /// the raster takes the pixel under each cell's centre. #[test] @@ -996,6 +1023,9 @@ mod tests { assert!(Grid::snapped(crs, [0.0, 0.0, 10.0, 10.0], bad, (0.0, 0.0)).is_err(), "{bad:?}"); } assert!(Grid::snapped(crs, [10.0, 0.0, 0.0, 10.0], (1.0, 1.0), (0.0, 0.0)).is_err()); + // An origin so far away that the snapped corners overflow is refused, not a 1 × 1 grid. + let far = Grid::snapped(crs, [-1e308, 0.0, 0.0, 10.0], (1e300, 1.0), (1e308, 0.0)).unwrap_err(); + assert!(far.contains("can't make a grid"), "{far}"); assert!(Grid::snapped(crs, [0.0, 0.0, 10.0, 10.0], (1.0, 1.0), (f64::NAN, 0.0)).is_err()); } diff --git a/crates/engine/src/cmd/imagery.rs b/crates/engine/src/cmd/imagery.rs index 3f442b7fa..7c2ac9f71 100644 --- a/crates/engine/src/cmd/imagery.rs +++ b/crates/engine/src/cmd/imagery.rs @@ -361,8 +361,9 @@ fn project_raster(s: &mut Session, p: &Value) -> CmdResult { let origin = match p.get("origin") { None | Some(Value::Null) => (0.0, 0.0), Some(v) => { - let n: Vec = v.as_array().map(|a| a.iter().filter_map(Value::as_f64).collect()).unwrap_or_default(); - match n.as_slice() { + // Every entry a number: one that isn't refuses the origin rather than drop out. + let n: Option> = v.as_array().map(|a| a.iter().map(Value::as_f64).collect()).unwrap_or_default(); + match n.as_deref().unwrap_or_default() { [x, y] if x.is_finite() && y.is_finite() => (*x, *y), _ => return Err("origin is [x, y] in the output's units".into()), } diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index f5df479b4..6f0ac43c7 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -17660,6 +17660,10 @@ fn a_raster_projected_onto_a_chosen_grid() { // Another origin moves the corners with it. let moved = s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 200, "origin": [50, 50], "name": "moved"})).unwrap(); + // An origin with anything but two numbers is refused, not read from the numbers left. + for bad in [json!(["x", 1, 2]), json!([0, 0, "EPSG"]), json!([1])] { + assert!(s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 200, "origin": bad})).is_err(), "{bad}"); + } assert_eq!((moved["corner"][0].as_f64().unwrap() - 50.0).rem_euclid(200.0), 0.0, "{moved}"); // Onto that grid exactly, by nearest neighbour (the default). From 9191cdbbad383a76eb578ce6c7d33852c7ccf995 Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Mon, 26 Oct 2026 12:55:50 -0300 Subject: [PATCH 3/7] Analysis: raster.project refuses an extent or grid of the wrong kind; reports the method used; mode counts the pixels that show (#821) Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/warp.rs | 27 ++++++++++++++++++++++++++- crates/engine/src/cmd/imagery.rs | 18 ++++++++++++++---- crates/engine/src/tests.rs | 7 +++++++ 3 files changed, 47 insertions(+), 5 deletions(-) diff --git a/crates/analysis/src/warp.rs b/crates/analysis/src/warp.rs index eed048b6b..40c3d37c2 100644 --- a/crates/analysis/src/warp.rs +++ b/crates/analysis/src/warp.rs @@ -413,6 +413,16 @@ impl Method { Method::Mode => "mode", } } + + /// The method `r` is projected with: a palette image's classes are never mixed, so it is + /// sampled at the nearest pixel, or takes the mode, whatever was asked. + pub fn used_on(self, r: &Raster) -> Method { + match (self, &r.color) { + (Method::Average | Method::Mode, ColorHint::Palette(_)) => Method::Mode, + (Method::Sample(_), ColorHint::Palette(_)) => Method::Sample(Resampling::Nearest), + (m, _) => m, + } + } } /// Projects `r` onto the grid `g`, through longitude/latitude when their coordinate systems @@ -679,7 +689,14 @@ fn aggregated( let at = oy * ow + ox; if inside { for ((o, v), k) in cell.iter_mut().zip(values.iter_mut()).zip(&weights) { - *o = if mode { mode_of(v) } else { mean_of(v, k, round) }; + *o = if mode { + // The most frequent colour among the pixels that show: see-through ones + // count for nothing, as in the mean. + let mut shown: Vec = v.iter().zip(k).filter(|(_, k)| **k > 0.0).map(|(x, _)| *x).collect(); + mode_of(&mut shown) + } else { + mean_of(v, k, round) + }; } } else { // A grid finer than the raster: the pixel under the cell's centre. @@ -967,6 +984,14 @@ mod tests { let g = Grid::snapped(r.crs, [0.0, 0.0, 2.0, 2.0], (2.0, 2.0), (0.0, 0.0)).unwrap(); let out = project(&r, &g, Method::Average, "a".into()).unwrap(); assert_eq!([0, 1, 2, 3].map(|b| out.value(b, 0, 0)), [Some(255.0), Some(0.0), Some(0.0), Some(128.0)]); + // Mode counts the pixels that show: three see-through black ones lose to one opaque red. + let wide = GeoTransform { origin_x: 0.0, origin_y: 1.0, pixel_width: 1.0, pixel_height: -1.0 }; + let px = |a: f32, b: f32| vec![a, b, b, b]; + let bands = vec![px(255.0, 0.0), px(0.0, 0.0), px(0.0, 0.0), px(255.0, 0.0)]; + let r = Raster::new("rgba", 4, 1, bands, wide, crs, None, ColorHint::Rgb { alpha: true }, "UInt8").unwrap(); + let g = Grid::snapped(r.crs, [0.0, 0.0, 4.0, 1.0], (4.0, 1.0), (0.0, 0.0)).unwrap(); + let out = project(&r, &g, Method::Mode, "m".into()).unwrap(); + assert_eq!([0, 1, 2].map(|b| out.value(b, 0, 0)), [Some(255.0), Some(0.0), Some(0.0)]); } /// A pixel centre on the edge between cells is in exactly one of them, and a grid finer than diff --git a/crates/engine/src/cmd/imagery.rs b/crates/engine/src/cmd/imagery.rs index 7c2ac9f71..b2f1456c6 100644 --- a/crates/engine/src/cmd/imagery.rs +++ b/crates/engine/src/cmd/imagery.rs @@ -327,6 +327,9 @@ fn warp_to_fit(s: &mut Session, p: &Value, layer: &str, r: &Raster) -> CmdResult /// `raster.project`: the layer on a grid of another coordinate system, cell size, origin and /// extent, or on another raster's grid, each cell sampled at its centre or taken from every /// pixel whose centre falls in it (see `rustgis_analysis::warp::project`). A new layer. +/// What an `extent` that isn't one is told. +const EXTENT: &str = "extent is a polygon layer's name, [west, south, east, north] in degrees, or \"view\""; + fn project_raster(s: &mut Session, p: &Value) -> CmdResult { let r = raster(s, p)?; r.require_coordinates()?; @@ -334,7 +337,12 @@ fn project_raster(s: &mut Session, p: &Value) -> CmdResult { None | Some(Value::Null) => warp::Method::Sample(warp::Resampling::Nearest), Some(v) => v.as_str().and_then(warp::Method::parse).ok_or("resampling is nearest, bilinear, cubic, average or mode")?, }; - let grid = match ps(p, "grid") { + // A grid named as anything but a raster layer's name is refused rather than ignored. + let template = match p.get("grid") { + None | Some(Value::Null) => None, + Some(v) => Some(v.as_str().ok_or("grid is a raster layer's name")?), + }; + let grid = match template { Some(template) => { if ["crs", "cellSize", "origin", "extent"].iter().any(|k| p.get(*k).is_some_and(|v| !v.is_null())) { return Err( @@ -384,13 +392,15 @@ fn project_raster(s: &mut Session, p: &Value) -> CmdResult { } warp::bounds_in([e.min_x, e.min_y, e.max_x, e.max_y], crs)? } - Some(v) => { - let area = super::catalogs::area(s, &json!({"bbox": v}))?; + Some(v @ (Value::Array(_) | Value::String(_))) => { + let area = super::catalogs::area(s, &json!({"bbox": v})).map_err(|_| EXTENT)?; if area[0] > area[2] { return Err("an extent across the antimeridian: give it one side at a time".into()); } warp::bounds_in(area, crs)? } + // A number, a flag or an object isn't an extent: refused, not read as the view. + Some(_) => return Err(EXTENT.into()), }; warp::Grid::snapped(crs, bounds, cell, origin)? } @@ -402,7 +412,7 @@ fn project_raster(s: &mut Session, p: &Value) -> CmdResult { let out_name = name(p, format!("{} projected", r.name)); let out = warp::project(&r, &grid, method, out_name.clone())?; let mut info = json!({ - "name": out_name, "crs": grid.crs.label(), "resampling": method.id(), "width": out.width, "height": out.height, "bands": out.band_count(), + "name": out_name, "crs": grid.crs.label(), "resampling": method.used_on(&r).id(), "width": out.width, "height": out.height, "bands": out.band_count(), "cellSize": [grid.cell_x, grid.cell_y], "corner": [grid.min_x, grid.max_y], }); if let Some((fw, fh)) = r.overview { diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index 6f0ac43c7..9bdf196b1 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -17664,6 +17664,13 @@ fn a_raster_projected_onto_a_chosen_grid() { for bad in [json!(["x", 1, 2]), json!([0, 0, "EPSG"]), json!([1])] { assert!(s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 200, "origin": bad})).is_err(), "{bad}"); } + // An extent or a grid of the wrong kind is refused, not read as the view or left out. + for bad in [json!({"extent": 5}), json!({"extent": {}}), json!({"extent": true}), json!({"grid": 5}), json!({"grid": ["the"]})] { + let mut p = bad.clone(); + p["layer"] = json!("the"); + let e = s.run("raster.project", &p).unwrap_err(); + assert!(e.contains("extent is") || e.contains("grid is"), "{bad}: {e}"); + } assert_eq!((moved["corner"][0].as_f64().unwrap() - 50.0).rem_euclid(200.0), 0.0, "{moved}"); // Onto that grid exactly, by nearest neighbour (the default). From f2accc44fa831e56e158e5040d8deabc1bc47e03 Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Mon, 26 Oct 2026 13:13:06 -0300 Subject: [PATCH 4/7] Analysis: raster.project's fourth review round (#821) An alpha with no value weighs nothing; the method that runs is resolved once; the extent error keeps its reason; the command's doc comment is back above it; the docs say how average and mode treat a colour image. Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/warp.rs | 11 +++++++++-- crates/engine/src/cmd/imagery.rs | 8 ++++---- docs/user-guide/analysis.md | 4 +++- 3 files changed, 16 insertions(+), 7 deletions(-) diff --git a/crates/analysis/src/warp.rs b/crates/analysis/src/warp.rs index 40c3d37c2..e0e530ae5 100644 --- a/crates/analysis/src/warp.rs +++ b/crates/analysis/src/warp.rs @@ -449,7 +449,8 @@ pub fn project(r: &Raster, g: &Grid, method: Method, name: String) -> Result sampled(r, g, canvas, m, &source, name), Method::Average | Method::Mode => aggregated(r, g, canvas, method == Method::Mode, &source, name), } @@ -677,7 +678,13 @@ fn aggregated( } inside = true; let idx = j * w + i; - let opacity = alpha.and_then(|a| r.bands.get(a)).and_then(|band| band.data.get(idx)).copied().unwrap_or(0.0); + // An alpha with no value (the raster's no-data) or below 0 weighs nothing. + let opacity = alpha + .and_then(|a| r.bands.get(a)) + .and_then(|band| band.data.get(idx)) + .copied() + .filter(|a| r.is_valid(*a) && *a > 0.0) + .unwrap_or(0.0); for (b, (v, k)) in values.iter_mut().zip(weights.iter_mut()).enumerate() { if let Some(x) = r.bands.get(b).and_then(|band| band.data.get(idx)).filter(|x| r.is_valid(**x)) { v.push(*x); diff --git a/crates/engine/src/cmd/imagery.rs b/crates/engine/src/cmd/imagery.rs index b2f1456c6..d30718e88 100644 --- a/crates/engine/src/cmd/imagery.rs +++ b/crates/engine/src/cmd/imagery.rs @@ -324,12 +324,12 @@ fn warp_to_fit(s: &mut Session, p: &Value, layer: &str, r: &Raster) -> CmdResult Ok(info) } -/// `raster.project`: the layer on a grid of another coordinate system, cell size, origin and -/// extent, or on another raster's grid, each cell sampled at its centre or taken from every -/// pixel whose centre falls in it (see `rustgis_analysis::warp::project`). A new layer. /// What an `extent` that isn't one is told. const EXTENT: &str = "extent is a polygon layer's name, [west, south, east, north] in degrees, or \"view\""; +/// `raster.project`: the layer on a grid of another coordinate system, cell size, origin and +/// extent, or on another raster's grid, each cell sampled at its centre or taken from every +/// pixel whose centre falls in it (see `rustgis_analysis::warp::project`). A new layer. fn project_raster(s: &mut Session, p: &Value) -> CmdResult { let r = raster(s, p)?; r.require_coordinates()?; @@ -393,7 +393,7 @@ fn project_raster(s: &mut Session, p: &Value) -> CmdResult { warp::bounds_in([e.min_x, e.min_y, e.max_x, e.max_y], crs)? } Some(v @ (Value::Array(_) | Value::String(_))) => { - let area = super::catalogs::area(s, &json!({"bbox": v})).map_err(|_| EXTENT)?; + let area = super::catalogs::area(s, &json!({"bbox": v})).map_err(|e| e.replacen("bbox", "extent", 1))?; if area[0] > area[2] { return Err("an extent across the antimeridian: give it one side at a time".into()); } diff --git a/docs/user-guide/analysis.md b/docs/user-guide/analysis.md index e01780c4e..1be3e3266 100644 --- a/docs/user-guide/analysis.md +++ b/docs/user-guide/analysis.md @@ -484,7 +484,9 @@ The Analyze tab's imagery groups open their tools in the Workflows pane: those tied), which is what a coarser grid of elevation or land-cover classes should hold. A cell no pixel centre falls in, on a grid finer than the raster, takes the pixel under its centre. Classes take `nearest` or `mode`, continuous values `bilinear` or `cubic`, or `average` - going coarser. A whole-number pixel type keeps its type, so its mean is rounded. + going coarser. A whole-number pixel type keeps its type, so its mean is rounded. On a colour + image the colours count as much as their pixels are opaque, and `mode` is taken band by band + (each channel's most frequent value among the pixels that show). - Cells off the raster are its no-data value. A raster without one keeps its bands as they are, with NaN there in a floating-point type and 0 in a whole-number one, as GDAL leaves them; a colour image gains an alpha band instead, so they are see-through. From e173c3b333ef8021d6046956edcee4307365692b Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Mon, 26 Oct 2026 13:15:38 -0300 Subject: [PATCH 5/7] Analysis: mode takes a colour image's alpha among the pixels that show (#821) Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/warp.rs | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/crates/analysis/src/warp.rs b/crates/analysis/src/warp.rs index e0e530ae5..d442b4ada 100644 --- a/crates/analysis/src/warp.rs +++ b/crates/analysis/src/warp.rs @@ -695,12 +695,14 @@ fn aggregated( } let at = oy * ow + ox; if inside { - for ((o, v), k) in cell.iter_mut().zip(values.iter_mut()).zip(&weights) { + for (b, ((o, v), k)) in cell.iter_mut().zip(values.iter_mut()).zip(&weights).enumerate() { *o = if mode { - // The most frequent colour among the pixels that show: see-through ones - // count for nothing, as in the mean. - let mut shown: Vec = v.iter().zip(k).filter(|(_, k)| **k > 0.0).map(|(x, _)| *x).collect(); - mode_of(&mut shown) + // The most frequent colour, and alpha, among the pixels that show: + // see-through ones count for nothing, as in the mean. A cell where none + // shows takes the mode of all (see-through, alpha 0). + let shows = |(x, k): (&f32, &f32)| if alpha == Some(b) { *x > 0.0 } else { *k > 0.0 }; + let mut shown: Vec = v.iter().zip(k).filter(|p| shows(*p)).map(|(x, _)| *x).collect(); + if shown.is_empty() { mode_of(v) } else { mode_of(&mut shown) } } else { mean_of(v, k, round) }; @@ -998,7 +1000,7 @@ mod tests { let r = Raster::new("rgba", 4, 1, bands, wide, crs, None, ColorHint::Rgb { alpha: true }, "UInt8").unwrap(); let g = Grid::snapped(r.crs, [0.0, 0.0, 4.0, 1.0], (4.0, 1.0), (0.0, 0.0)).unwrap(); let out = project(&r, &g, Method::Mode, "m".into()).unwrap(); - assert_eq!([0, 1, 2].map(|b| out.value(b, 0, 0)), [Some(255.0), Some(0.0), Some(0.0)]); + assert_eq!([0, 1, 2, 3].map(|b| out.value(b, 0, 0)), [Some(255.0), Some(0.0), Some(0.0), Some(255.0)], "red, and showing"); } /// A pixel centre on the edge between cells is in exactly one of them, and a grid finer than From d042fd8806d9f199ccc344d19b01911516be6b45 Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Mon, 26 Oct 2026 13:29:44 -0300 Subject: [PATCH 6/7] Analysis: raster.project reports the name its layer got; mode takes no colour from a pixel that doesn't show (#821) Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/warp.rs | 26 ++++++++++++++++++-------- crates/engine/src/cmd/imagery.rs | 5 +++++ crates/engine/src/tests.rs | 7 +++++++ 3 files changed, 30 insertions(+), 8 deletions(-) diff --git a/crates/analysis/src/warp.rs b/crates/analysis/src/warp.rs index d442b4ada..a4dd7fdbe 100644 --- a/crates/analysis/src/warp.rs +++ b/crates/analysis/src/warp.rs @@ -702,7 +702,13 @@ fn aggregated( // shows takes the mode of all (see-through, alpha 0). let shows = |(x, k): (&f32, &f32)| if alpha == Some(b) { *x > 0.0 } else { *k > 0.0 }; let mut shown: Vec = v.iter().zip(k).filter(|p| shows(*p)).map(|(x, _)| *x).collect(); - if shown.is_empty() { mode_of(v) } else { mode_of(&mut shown) } + match (shown.is_empty(), alpha) { + (false, _) => mode_of(&mut shown), + // No pixel that shows has a colour here: none is taken from one that + // doesn't. + (true, Some(a)) if a != b => f32::NAN, + (true, _) => mode_of(v), + } } else { mean_of(v, k, round) }; @@ -752,14 +758,12 @@ fn mean_of(v: &[f32], k: &[f32], round: bool) -> f32 { fn mode_of(v: &mut [f32]) -> f32 { v.sort_unstable_by(f32::total_cmp); let (mut best, mut best_n) = (f32::NAN, 0usize); - let mut k = 0; - while let Some(x) = v.get(k).copied() { - let n = v.iter().skip(k).take_while(|y| **y == x).count(); - // Sorted ascending: a later value only wins with more. - if n > best_n { - (best, best_n) = (x, n); + // Sorted ascending: a later run only wins with more, so ties keep the least value. Each run + // holds one value at least (a NaN is a run of its own), so this always moves on. + for run in v.chunk_by(|a, b| a == b) { + if run.len() > best_n { + (best, best_n) = (run[0], run.len()); } - k += n; } best } @@ -1001,6 +1005,12 @@ mod tests { let g = Grid::snapped(r.crs, [0.0, 0.0, 4.0, 1.0], (4.0, 1.0), (0.0, 0.0)).unwrap(); let out = project(&r, &g, Method::Mode, "m".into()).unwrap(); assert_eq!([0, 1, 2, 3].map(|b| out.value(b, 0, 0)), [Some(255.0), Some(0.0), Some(0.0), Some(255.0)], "red, and showing"); + // The opaque pixel has no red: the see-through pixels' red isn't taken instead. + let bands = vec![px(-1.0, 90.0), px(0.0, 0.0), px(0.0, 0.0), px(255.0, 0.0)]; + let r = Raster::new("rgba", 4, 1, bands, wide, crs, Some(-1.0), ColorHint::Rgb { alpha: true }, "Int16").unwrap(); + let out = project(&r, &g, Method::Mode, "m".into()).unwrap(); + // (A cell without a value is cleared to 0 in an image with alpha, as outside the raster.) + assert_eq!([0, 3].map(|b| out.value(b, 0, 0)), [Some(0.0), Some(255.0)], "not the see-through 90"); } /// A pixel centre on the edge between cells is in exactly one of them, and a grid finer than diff --git a/crates/engine/src/cmd/imagery.rs b/crates/engine/src/cmd/imagery.rs index d30718e88..fcf7d3784 100644 --- a/crates/engine/src/cmd/imagery.rs +++ b/crates/engine/src/cmd/imagery.rs @@ -440,6 +440,11 @@ fn project_raster(s: &mut Session, p: &Value) -> CmdResult { o.entry(k.clone()).or_insert_with(|| val.clone()); } } + // The name the layer got (another may hold the one asked for: "… projected (2)"). + let id = v.get("layer").and_then(Value::as_u64).and_then(|id| u32::try_from(id).ok()); + if let Some(name) = id.and_then(|id| s.map().layer(id)).map(|l| l.name.clone()) { + info["name"] = json!(name); + } Ok(info) } diff --git a/crates/engine/src/tests.rs b/crates/engine/src/tests.rs index 9bdf196b1..9eb750f5b 100644 --- a/crates/engine/src/tests.rs +++ b/crates/engine/src/tests.rs @@ -17664,6 +17664,13 @@ fn a_raster_projected_onto_a_chosen_grid() { for bad in [json!(["x", 1, 2]), json!([0, 0, "EPSG"]), json!([1])] { assert!(s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 200, "origin": bad})).is_err(), "{bad}"); } + // Projected again with the default name: the name reported is the one the layer got. + let first = s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 500})).unwrap(); + let again = s.run("raster.project", &json!({"layer": "the", "crs": 31983, "cellSize": 500})).unwrap(); + assert_ne!(first["name"], again["name"]); + for r in [&first, &again] { + assert!(s.map().layer_by_name(r["name"].as_str().unwrap()).is_some(), "{r}"); + } // An extent or a grid of the wrong kind is refused, not read as the view or left out. for bad in [json!({"extent": 5}), json!({"extent": {}}), json!({"extent": true}), json!({"grid": 5}), json!({"grid": ["the"]})] { let mut p = bad.clone(); From 2253d8e2ea99c8f6edd36c1940730c4c35ad18b2 Mon Sep 17 00:00:00 2001 From: elpidioassinaturas Date: Mon, 26 Oct 2026 13:39:57 -0300 Subject: [PATCH 7/7] Analysis: average sums a cell as its pixels come, keeping no copy of them (#821) Co-Authored-By: Claude Opus 5.5 --- crates/analysis/src/warp.rs | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/crates/analysis/src/warp.rs b/crates/analysis/src/warp.rs index a4dd7fdbe..5d394515e 100644 --- a/crates/analysis/src/warp.rs +++ b/crates/analysis/src/warp.rs @@ -651,6 +651,7 @@ fn aggregated( let mut values: Vec> = vec![Vec::new(); canvas.nb]; // Each value's weight: a colour counts as much as its pixel is opaque, as bilinear weighs it. let mut weights: Vec> = vec![Vec::new(); canvas.nb]; + let mut sums: Vec<(f64, f64)> = vec![(0.0, 0.0); canvas.nb]; let alpha = canvas.alpha_src; let mut cell = vec![f32::NAN; canvas.nb]; for oy in 0..oh { @@ -670,6 +671,8 @@ fn aggregated( let (i0, i1, j0, j1) = (first(c0), past(c1, w), first(r0), past(r1, h)); values.iter_mut().for_each(Vec::clear); weights.iter_mut().for_each(Vec::clear); + // The mean is summed as the pixels come (no copy of the cell's values); mode keeps them. + sums.iter_mut().for_each(|s| *s = (0.0, 0.0)); let mut inside = false; for j in j0..j1 { for i in i0..i1 { @@ -685,17 +688,22 @@ fn aggregated( .copied() .filter(|a| r.is_valid(*a) && *a > 0.0) .unwrap_or(0.0); - for (b, (v, k)) in values.iter_mut().zip(weights.iter_mut()).enumerate() { + for (b, ((v, k), sum)) in values.iter_mut().zip(weights.iter_mut()).zip(sums.iter_mut()).enumerate() { if let Some(x) = r.bands.get(b).and_then(|band| band.data.get(idx)).filter(|x| r.is_valid(**x)) { - v.push(*x); - k.push(if alpha.is_some_and(|a| a != b) { opacity } else { 1.0 }); + let weight = if alpha.is_some_and(|a| a != b) { opacity } else { 1.0 }; + if mode { + v.push(*x); + k.push(weight); + } else { + *sum = (sum.0 + f64::from(*x) * f64::from(weight), sum.1 + f64::from(weight)); + } } } } } let at = oy * ow + ox; if inside { - for (b, ((o, v), k)) in cell.iter_mut().zip(values.iter_mut()).zip(&weights).enumerate() { + for (b, (((o, v), k), sum)) in cell.iter_mut().zip(values.iter_mut()).zip(&weights).zip(&sums).enumerate() { *o = if mode { // The most frequent colour, and alpha, among the pixels that show: // see-through ones count for nothing, as in the mean. A cell where none @@ -710,7 +718,7 @@ fn aggregated( (true, _) => mode_of(v), } } else { - mean_of(v, k, round) + mean_of(*sum, round) }; } } else { @@ -743,10 +751,10 @@ fn holds(q: &[(f64, f64); 4], p: (f64, f64)) -> bool { side(a, b) >= 0.0 && side(b, c) > 0.0 && side(c, d) > 0.0 && side(d, a) >= 0.0 } -/// The mean of `v`, each value weighed by its `k` (a colour by its pixel's alpha), rounded for a -/// whole-number pixel type; NaN when nothing weighs (empty, or all see-through). -fn mean_of(v: &[f32], k: &[f32], round: bool) -> f32 { - let (sum, weight) = v.iter().zip(k).fold((0.0f64, 0.0f64), |(s, w), (x, k)| (s + f64::from(*x) * f64::from(*k), w + f64::from(*k))); +/// The mean of a cell's values from their weighed sum and total weight (a colour weighs as its +/// pixel's alpha), rounded for a whole-number pixel type; NaN when nothing weighs (no value, or +/// all see-through). +fn mean_of((sum, weight): (f64, f64), round: bool) -> f32 { if weight <= 1e-12 { return f32::NAN; } @@ -761,8 +769,8 @@ fn mode_of(v: &mut [f32]) -> f32 { // Sorted ascending: a later run only wins with more, so ties keep the least value. Each run // holds one value at least (a NaN is a run of its own), so this always moves on. for run in v.chunk_by(|a, b| a == b) { - if run.len() > best_n { - (best, best_n) = (run[0], run.len()); + if let Some(&x) = run.first().filter(|_| run.len() > best_n) { + (best, best_n) = (x, run.len()); } } best