diff --git a/crates/analysis/src/warp.rs b/crates/analysis/src/warp.rs index 76db77772..5d394515e 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,524 @@ 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; + // 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!( + "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", + } + } + + /// 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 +/// 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)) + }; + // The method that runs, as the result reports it (a palette image's classes are never mixed). + match method.used_on(r) { + 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]; + // 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 { + 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); + 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 { + if !holds(&quad, (i as f64 + 0.5, j as f64 + 0.5)) { + continue; + } + inside = true; + let idx = j * w + i; + // 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), 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)) { + 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)); + } + } + } } } - 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 (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 + // 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(); + 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(*sum, 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) +} - 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()); +/// 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; } - Ok(out) + 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 +} + +/// 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; + } + let m = (sum / weight) 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); + // 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 let Some(&x) = run.first().filter(|_| run.len() > best_n) { + (best, best_n) = (x, run.len()); + } + } + best } #[cfg(test)] @@ -562,4 +957,155 @@ 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)]); + } + + /// 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)]); + // 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, 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 + /// 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()); + // 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()); + } + + /// 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..fcf7d3784 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,130 @@ fn warp_to_fit(s: &mut Session, p: &Value, layer: &str, r: &Raster) -> CmdResult Ok(info) } +/// 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()?; + 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")?, + }; + // 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( + "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) => { + // 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()), + } + } + }; + 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 @ (Value::Array(_) | Value::String(_))) => { + 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()); + } + 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)? + } + }; + 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.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 { + 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()); + } + } + // 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) +} + 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..9eb750f5b 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,102 @@ 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(); + // 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}"); + } + // 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(); + 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). + 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..1be3e3266 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,38 @@ 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. 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. + - 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 +544,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'