trailgen_core/

raster.rs

use crate::crs::{UtmCrs, geographic_to_utm, wgs84_to_web_mercator};
use crate::enrich::{ElevationSample, ElevationSampler};
use crate::geo::Coord;
use crate::model::Provenance;
use crate::{Result, TrailgenError};
use num_traits::ToPrimitive;
use serde::{Deserialize, Serialize};
use std::fs::File;
use std::iter::Peekable;
use std::path::{Path, PathBuf};
use std::str::FromStr;
use tiff::decoder::{Decoder, DecodingResult};
use tiff::tags::Tag;
use tiff::{TiffError, TiffFormatError};

#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
pub struct ArcAsciiGrid {
    pub ncols: usize,
    pub nrows: usize,
    pub xllcorner: f64,
    pub yllcorner: f64,
    pub cellsize: f64,
    pub nodata_value: f64,
    pub confidence: f64,
    pub provenance: Provenance,
    values: Vec<f64>,
}

#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
pub struct GeoTiffDem {
    pub width: usize,
    pub height: usize,
    pub crs: RasterCrs,
    #[serde(alias = "origin_lon")]
    pub origin_x: f64,
    #[serde(alias = "origin_lat")]
    pub origin_y: f64,
    #[serde(alias = "pixel_width_deg")]
    pub pixel_width: f64,
    #[serde(alias = "pixel_height_deg")]
    pub pixel_height: f64,
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub transform: Option<RasterTransform>,
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub nodata_value: Option<f64>,
    pub confidence: f64,
    pub provenance: Provenance,
    values: Vec<f64>,
}

#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
pub struct RasterTransform {
    pub x0: f64,
    pub y0: f64,
    pub dx_col: f64,
    pub dy_col: f64,
    pub dx_row: f64,
    pub dy_row: f64,
}

impl RasterTransform {
    fn north_up(origin_x: f64, origin_y: f64, pixel_width: f64, pixel_height: f64) -> Self {
        Self {
            x0: origin_x,
            y0: origin_y,
            dx_col: pixel_width,
            dy_col: 0.0,
            dx_row: 0.0,
            dy_row: -pixel_height,
        }
    }

    fn from_model_transformation(xs: &[f64]) -> Result<Self> {
        if xs.len() < 16 {
            return Err(TrailgenError::InvalidData(
                "GeoTIFF ModelTransformationTag must contain sixteen numbers".to_owned(),
            ));
        }
        let transform = Self {
            x0: xs[3],
            y0: xs[7],
            dx_col: xs[0],
            dy_col: xs[4],
            dx_row: xs[1],
            dy_row: xs[5],
        };
        transform.validate("GeoTIFF ModelTransformationTag")?;
        Ok(transform)
    }

    fn validate(self, label: &str) -> Result<()> {
        let finite = [
            self.x0,
            self.y0,
            self.dx_col,
            self.dy_col,
            self.dx_row,
            self.dy_row,
        ]
        .into_iter()
        .all(f64::is_finite);
        if !finite || self.det().abs() <= f64::EPSILON {
            return Err(TrailgenError::InvalidData(format!(
                "{label} must be finite and invertible"
            )));
        }
        Ok(())
    }

    fn pixel_xy(self, x: f64, y: f64) -> (f64, f64) {
        let dx = x - self.x0;
        let dy = y - self.y0;
        let det = self.det();
        (
            (dx.mul_add(self.dy_row, -self.dx_row * dy)) / det,
            (self.dx_col.mul_add(dy, -dx * self.dy_col)) / det,
        )
    }

    const fn det(self) -> f64 {
        self.dx_col * self.dy_row - self.dx_row * self.dy_col
    }

    fn pixel_width(self) -> f64 {
        self.dx_col.hypot(self.dy_col)
    }

    fn pixel_height(self) -> f64 {
        self.dx_row.hypot(self.dy_row)
    }
}

#[derive(Clone, Copy, Debug, Eq, PartialEq, Serialize, Deserialize)]
#[serde(rename_all = "kebab-case")]
pub enum RasterCrs {
    Wgs84Degrees,
    Nad83Degrees,
    WebMercatorMeters,
    #[serde(alias = "wgs84-utm-meters")]
    UtmMeters(UtmCrs),
}

impl RasterCrs {
    fn xy(self, coord: Coord) -> Option<(f64, f64)> {
        match self {
            Self::Wgs84Degrees | Self::Nad83Degrees => Some((coord.lon, coord.lat)),
            Self::WebMercatorMeters => Some(wgs84_to_web_mercator(coord)),
            Self::UtmMeters(crs) => geographic_to_utm(coord, crs),
        }
    }
}

const fn default_raster_crs() -> RasterCrs {
    RasterCrs::Wgs84Degrees
}

#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
pub struct VrtDem {
    pub width: usize,
    pub height: usize,
    #[serde(default = "default_raster_crs")]
    pub crs: RasterCrs,
    #[serde(alias = "origin_lon")]
    pub origin_x: f64,
    #[serde(alias = "origin_lat")]
    pub origin_y: f64,
    #[serde(alias = "pixel_width_deg")]
    pub pixel_width: f64,
    #[serde(alias = "pixel_height_deg")]
    pub pixel_height: f64,
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub transform: Option<RasterTransform>,
    #[serde(default, skip_serializing_if = "Option::is_none")]
    pub nodata_value: Option<f64>,
    pub source_filename: String,
    pub confidence: f64,
    pub provenance: Provenance,
    values: Vec<f64>,
}

#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
#[serde(tag = "kind", rename_all = "kebab-case")]
pub enum RasterDem {
    ArcAscii(ArcAsciiGrid),
    GeoTiff(GeoTiffDem),
    Vrt(VrtDem),
}

impl ArcAsciiGrid {
    pub fn parse(raw: &str, provenance: Provenance, confidence: f64) -> Result<Self> {
        let mut lines = raw
            .lines()
            .map(str::trim)
            .filter(|line| !line.is_empty())
            .peekable();
        let ncols: usize = header(&mut lines, "ncols")?;
        let nrows: usize = header(&mut lines, "nrows")?;
        let xllcorner: f64 = header(&mut lines, "xllcorner")?;
        let yllcorner: f64 = header(&mut lines, "yllcorner")?;
        let cellsize: f64 = header(&mut lines, "cellsize")?;
        let nodata_value: f64 = optional_header(&mut lines, "nodata_value")?.unwrap_or(-9999.0);
        if ncols == 0 || nrows == 0 {
            return Err(TrailgenError::InvalidData(
                "raster dimensions must be nonzero".to_owned(),
            ));
        }
        if !xllcorner.is_finite()
            || !yllcorner.is_finite()
            || !cellsize.is_finite()
            || cellsize <= 0.0
            || !nodata_value.is_finite()
        {
            return Err(TrailgenError::InvalidData(
                "raster georeference headers must be finite and cellsize positive".to_owned(),
            ));
        }
        let values = lines
            .flat_map(str::split_whitespace)
            .map(|x| {
                x.parse::<f64>()
                    .map_err(|e| TrailgenError::InvalidData(format!("invalid raster cell: {e}")))
            })
            .collect::<Result<Vec<_>>>()?;
        let expected = ncols
            .checked_mul(nrows)
            .ok_or_else(|| TrailgenError::InvalidData("raster dimensions overflow".to_owned()))?;
        if values.len() != expected {
            return Err(TrailgenError::InvalidData(format!(
                "raster has {} cells, expected {expected}",
                values.len()
            )));
        }
        Ok(Self {
            ncols,
            nrows,
            xllcorner,
            yllcorner,
            cellsize,
            nodata_value,
            confidence: confidence.clamp(0.0, 1.0),
            provenance,
            values,
        })
    }

    #[must_use]
    pub fn contains(&self, coord: Coord) -> bool {
        let xmax = self
            .cellsize
            .mul_add(usize_to_f64(self.ncols), self.xllcorner);
        let ymax = self
            .cellsize
            .mul_add(usize_to_f64(self.nrows), self.yllcorner);
        (self.xllcorner..=xmax).contains(&coord.lon) && (self.yllcorner..=ymax).contains(&coord.lat)
    }

    fn interpolated_elevation_m(&self, coord: Coord) -> Option<f64> {
        if !self.contains(coord) {
            return None;
        }
        let col = self.x_center_index(coord.lon);
        let row = self.y_center_index(coord.lat);
        let col0 = f64_to_index(col.floor(), self.ncols)?;
        let col1 = f64_to_index(col.ceil(), self.ncols)?;
        let row0 = f64_to_index(row.floor(), self.nrows)?;
        let row1 = f64_to_index(row.ceil(), self.nrows)?;
        let tx = col - usize_to_f64(col0);
        let ty = row - usize_to_f64(row0);
        let z00 = self.cell(row0, col0)?;
        let z01 = self.cell(row0, col1)?;
        let z10 = self.cell(row1, col0)?;
        let z11 = self.cell(row1, col1)?;
        let top = tx.mul_add(z01 - z00, z00);
        let bottom = tx.mul_add(z11 - z10, z10);
        Some(ty.mul_add(bottom - top, top))
    }

    fn x_center_index(&self, lon: f64) -> f64 {
        ((lon - self.xllcorner) / self.cellsize - 0.5)
            .clamp(0.0, usize_to_f64(self.ncols.saturating_sub(1)))
    }

    fn y_center_index(&self, lat: f64) -> f64 {
        let ymax = self
            .cellsize
            .mul_add(usize_to_f64(self.nrows), self.yllcorner);
        ((ymax - lat) / self.cellsize - 0.5).clamp(0.0, usize_to_f64(self.nrows.saturating_sub(1)))
    }

    fn cell(&self, row: usize, col: usize) -> Option<f64> {
        let value = *self
            .values
            .get(row.checked_mul(self.ncols)?.checked_add(col)?)?;
        if (value - self.nodata_value).abs() <= f64::EPSILON {
            return None;
        }
        Some(value)
    }
}

impl GeoTiffDem {
    pub fn from_path(path: &Path, provenance: Provenance, confidence: f64) -> Result<Self> {
        let file = File::open(path)
            .map_err(|error| TrailgenError::InvalidData(format!("open GeoTIFF: {error}")))?;
        let mut decoder = Decoder::new(file)
            .map_err(|error| TrailgenError::InvalidData(format!("decode GeoTIFF: {error}")))?;
        let (width, height) = decoder.dimensions().map_err(|error| {
            TrailgenError::InvalidData(format!("read GeoTIFF dimensions: {error}"))
        })?;
        let width = usize::try_from(width)
            .map_err(|_| TrailgenError::InvalidData("GeoTIFF width overflow".to_owned()))?;
        let height = usize::try_from(height)
            .map_err(|_| TrailgenError::InvalidData("GeoTIFF height overflow".to_owned()))?;
        let expected = width
            .checked_mul(height)
            .ok_or_else(|| TrailgenError::InvalidData("GeoTIFF dimensions overflow".to_owned()))?;
        let georef = GeoTiffGeoref::read(&mut decoder)?;
        let nodata_value = optional_ascii(&mut decoder, Tag::GdalNodata)?
            .and_then(|value| value.trim_matches('\0').trim().parse::<f64>().ok());
        let image = decoder
            .read_image()
            .map_err(|error| TrailgenError::InvalidData(format!("read GeoTIFF pixels: {error}")))?;
        let values = decoding_result_to_f64(image);
        if values.len() != expected {
            return Err(TrailgenError::UnsupportedFormat(format!(
                "GeoTIFF DEM must be single-band; decoded {} sample(s), expected {expected}",
                values.len()
            )));
        }
        Ok(Self {
            width,
            height,
            crs: georef.crs,
            origin_x: georef.origin_x,
            origin_y: georef.origin_y,
            pixel_width: georef.pixel_width,
            pixel_height: georef.pixel_height,
            transform: georef.transform,
            nodata_value,
            confidence: confidence.clamp(0.0, 1.0),
            provenance,
            values,
        })
    }

    #[must_use]
    pub fn contains(&self, coord: Coord) -> bool {
        raster_contains(self.width, self.height, self.crs, self.transform(), coord)
    }

    fn interpolated_elevation_m(&self, coord: Coord) -> Option<f64> {
        interpolated_raster_value(
            self.width,
            self.height,
            &self.values,
            self.nodata_value,
            self.crs,
            self.transform(),
            coord,
        )
    }

    fn transform(&self) -> RasterTransform {
        self.transform.unwrap_or_else(|| {
            RasterTransform::north_up(
                self.origin_x,
                self.origin_y,
                self.pixel_width,
                self.pixel_height,
            )
        })
    }
}

impl VrtDem {
    pub fn from_path(path: &Path, provenance: Provenance, confidence: f64) -> Result<Self> {
        let raw = std::fs::read_to_string(path)
            .map_err(|error| TrailgenError::InvalidData(format!("read VRT: {error}")))?;
        let spec = VrtSpec::parse(path, &raw)?;
        let source = GeoTiffDem::from_path(
            &spec.source_path,
            Provenance {
                source: "vrt-source-geotiff".to_owned(),
                layer: Some("vrt-source".to_owned()),
                source_id: spec
                    .source_path
                    .file_name()
                    .and_then(|name| name.to_str())
                    .map(str::to_owned),
                license: None,
            },
            confidence,
        )?;
        let crs = spec.crs.unwrap_or(source.crs);
        if spec.crs.is_some_and(|spec_crs| spec_crs != source.crs) {
            return Err(TrailgenError::UnsupportedFormat(format!(
                "VRT SRS {:?} does not match source GeoTIFF CRS {:?}; reproject or materialize a consistent VRT",
                crs, source.crs
            )));
        }
        if source.width != spec.width || source.height != spec.height {
            return Err(TrailgenError::UnsupportedFormat(format!(
                "VRT source raster dimensions {}x{} do not match VRT {}x{}",
                source.width, source.height, spec.width, spec.height
            )));
        }
        Ok(Self {
            width: spec.width,
            height: spec.height,
            crs,
            origin_x: spec.transform.x0,
            origin_y: spec.transform.y0,
            pixel_width: spec.transform.pixel_width(),
            pixel_height: spec.transform.pixel_height(),
            transform: Some(spec.transform),
            nodata_value: spec.nodata_value.or(source.nodata_value),
            source_filename: spec.source_path.display().to_string(),
            confidence: confidence.clamp(0.0, 1.0),
            provenance,
            values: source.values,
        })
    }

    pub fn referenced_sources(path: &Path) -> Result<Vec<PathBuf>> {
        let raw = std::fs::read_to_string(path)
            .map_err(|error| TrailgenError::InvalidData(format!("read VRT: {error}")))?;
        Ok(vec![VrtSpec::parse(path, &raw)?.source_path])
    }

    #[must_use]
    pub fn contains(&self, coord: Coord) -> bool {
        raster_contains(self.width, self.height, self.crs, self.transform(), coord)
    }

    fn interpolated_elevation_m(&self, coord: Coord) -> Option<f64> {
        interpolated_raster_value(
            self.width,
            self.height,
            &self.values,
            self.nodata_value,
            self.crs,
            self.transform(),
            coord,
        )
    }

    fn transform(&self) -> RasterTransform {
        self.transform.unwrap_or_else(|| {
            RasterTransform::north_up(
                self.origin_x,
                self.origin_y,
                self.pixel_width,
                self.pixel_height,
            )
        })
    }
}

impl ElevationSampler for ArcAsciiGrid {
    fn sample(&self, coord: Coord) -> Option<ElevationSample> {
        let value = self.interpolated_elevation_m(coord)?;
        Some(ElevationSample {
            ele_m: value,
            confidence: self.confidence,
            provenance: self.provenance.clone(),
        })
    }
}

impl ElevationSampler for RasterDem {
    fn sample(&self, coord: Coord) -> Option<ElevationSample> {
        match self {
            Self::ArcAscii(raster) => raster.sample(coord),
            Self::GeoTiff(raster) => raster.sample(coord),
            Self::Vrt(raster) => raster.sample(coord),
        }
    }
}

impl ElevationSampler for GeoTiffDem {
    fn sample(&self, coord: Coord) -> Option<ElevationSample> {
        let value = self.interpolated_elevation_m(coord)?;
        Some(ElevationSample {
            ele_m: value,
            confidence: self.confidence,
            provenance: self.provenance.clone(),
        })
    }
}

impl ElevationSampler for VrtDem {
    fn sample(&self, coord: Coord) -> Option<ElevationSample> {
        let value = self.interpolated_elevation_m(coord)?;
        Some(ElevationSample {
            ele_m: value,
            confidence: self.confidence,
            provenance: self.provenance.clone(),
        })
    }
}

struct VrtSpec {
    width: usize,
    height: usize,
    crs: Option<RasterCrs>,
    transform: RasterTransform,
    nodata_value: Option<f64>,
    source_path: PathBuf,
}

impl VrtSpec {
    fn parse(path: &Path, raw: &str) -> Result<Self> {
        let doc = roxmltree::Document::parse(raw)
            .map_err(|error| TrailgenError::InvalidData(format!("parse VRT XML: {error}")))?;
        let root = doc.root_element();
        if root.tag_name().name() != "VRTDataset" {
            return Err(TrailgenError::UnsupportedFormat(
                "raster VRT must have <VRTDataset> root".to_owned(),
            ));
        }
        let width = required_attr::<usize>(root, "rasterXSize")?;
        let height = required_attr::<usize>(root, "rasterYSize")?;
        if width == 0 || height == 0 {
            return Err(TrailgenError::InvalidData(
                "VRT raster dimensions must be nonzero".to_owned(),
            ));
        }
        let geotransform = child_text(root, "GeoTransform")
            .ok_or_else(|| TrailgenError::InvalidData("VRT missing GeoTransform".to_owned()))?
            .split(',')
            .map(str::trim)
            .map(str::parse::<f64>)
            .collect::<std::result::Result<Vec<_>, _>>()
            .map_err(|error| {
                TrailgenError::InvalidData(format!("invalid VRT GeoTransform: {error}"))
            })?;
        if geotransform.len() != 6 {
            return Err(TrailgenError::InvalidData(
                "VRT GeoTransform must contain six numbers".to_owned(),
            ));
        }
        let transform = RasterTransform {
            x0: geotransform[0],
            y0: geotransform[3],
            dx_col: geotransform[1],
            dy_col: geotransform[4],
            dx_row: geotransform[2],
            dy_row: geotransform[5],
        };
        transform.validate("VRT GeoTransform")?;
        let crs = child_text(root, "SRS").map(parse_vrt_srs).transpose()?;
        let band = root
            .children()
            .find(|node| node.has_tag_name("VRTRasterBand"))
            .ok_or_else(|| TrailgenError::InvalidData("VRT missing VRTRasterBand".to_owned()))?;
        let simple = band
            .children()
            .find(|node| node.has_tag_name("SimpleSource"))
            .ok_or_else(|| {
                TrailgenError::UnsupportedFormat("VRT DEM requires SimpleSource".to_owned())
            })?;
        ensure_identity_rect(simple, "SrcRect", width, height)?;
        ensure_identity_rect(simple, "DstRect", width, height)?;
        let source_filename = simple
            .children()
            .find(|node| node.has_tag_name("SourceFilename"))
            .ok_or_else(|| {
                TrailgenError::InvalidData("VRT SimpleSource missing SourceFilename".to_owned())
            })?;
        let source_text = source_filename
            .text()
            .map(str::trim)
            .filter(|text| !text.is_empty())
            .ok_or_else(|| TrailgenError::InvalidData("VRT SourceFilename is empty".to_owned()))?;
        let source_path = resolve_vrt_source(path, source_filename, source_text);
        let source_band = child_text(simple, "SourceBand").unwrap_or("1").trim();
        if source_band != "1" {
            return Err(TrailgenError::UnsupportedFormat(
                "VRT DEM only supports SourceBand 1".to_owned(),
            ));
        }
        let nodata_value = child_text(band, "NoDataValue")
            .or_else(|| child_text(root, "NoDataValue"))
            .map(str::trim)
            .map(str::parse::<f64>)
            .transpose()
            .map_err(|error| {
                TrailgenError::InvalidData(format!("invalid VRT NoDataValue: {error}"))
            })?;
        Ok(Self {
            width,
            height,
            crs,
            transform,
            nodata_value,
            source_path,
        })
    }
}

fn raster_contains(
    width: usize,
    height: usize,
    crs: RasterCrs,
    transform: RasterTransform,
    coord: Coord,
) -> bool {
    let Some((x, y)) = crs.xy(coord) else {
        return false;
    };
    let (col, row) = transform.pixel_xy(x, y);
    (0.0..=usize_to_f64(width)).contains(&col) && (0.0..=usize_to_f64(height)).contains(&row)
}

fn interpolated_raster_value(
    width: usize,
    height: usize,
    values: &[f64],
    nodata_value: Option<f64>,
    crs: RasterCrs,
    transform: RasterTransform,
    coord: Coord,
) -> Option<f64> {
    if !raster_contains(width, height, crs, transform, coord) {
        return None;
    }
    let (x, y) = crs.xy(coord)?;
    let (col_px, row_px) = transform.pixel_xy(x, y);
    let col = (col_px - 0.5).clamp(0.0, usize_to_f64(width.saturating_sub(1)));
    let row = (row_px - 0.5).clamp(0.0, usize_to_f64(height.saturating_sub(1)));
    let col0 = f64_to_index(col.floor(), width)?;
    let col1 = f64_to_index(col.ceil(), width)?;
    let row0 = f64_to_index(row.floor(), height)?;
    let row1 = f64_to_index(row.ceil(), height)?;
    let tx = col - usize_to_f64(col0);
    let ty = row - usize_to_f64(row0);
    let z00 = raster_cell(width, values, nodata_value, row0, col0)?;
    let z01 = raster_cell(width, values, nodata_value, row0, col1)?;
    let z10 = raster_cell(width, values, nodata_value, row1, col0)?;
    let z11 = raster_cell(width, values, nodata_value, row1, col1)?;
    let top = tx.mul_add(z01 - z00, z00);
    let bottom = tx.mul_add(z11 - z10, z10);
    Some(ty.mul_add(bottom - top, top))
}

fn raster_cell(
    width: usize,
    values: &[f64],
    nodata_value: Option<f64>,
    row: usize,
    col: usize,
) -> Option<f64> {
    let value = *values.get(row.checked_mul(width)?.checked_add(col)?)?;
    if nodata_value.is_some_and(|nodata| (value - nodata).abs() <= f64::EPSILON) {
        return None;
    }
    value.is_finite().then_some(value)
}

fn parse_vrt_srs(srs: &str) -> Result<RasterCrs> {
    let normalized: String = srs
        .chars()
        .filter(char::is_ascii_alphanumeric)
        .flat_map(char::to_uppercase)
        .collect();
    if normalized.contains("EPSG3857")
        || normalized.contains("EPSG900913")
        || normalized.contains("WEBMERCATOR")
        || normalized.contains("PSEUDOMERCATOR")
    {
        Ok(RasterCrs::WebMercatorMeters)
    } else if let Some(crs) = UtmCrs::from_normalized_srs(&normalized) {
        Ok(RasterCrs::UtmMeters(crs))
    } else if normalized.contains("PROJCS")
        || normalized.contains("PROJCRS")
        || normalized.contains("PROJECTION")
    {
        Err(TrailgenError::UnsupportedFormat(
            "projected VRT SRS must be EPSG:3857 Web Mercator or WGS84/NAD83 UTM (EPSG:326xx/327xx/269xx), or reprojected before ingestion"
                .to_owned(),
        ))
    } else if normalized.contains("EPSG4326")
        || normalized.contains("OGC13CRS84")
        || normalized.contains("OGC14CRS84")
        || normalized.contains("CRS84")
        || normalized.contains("WGS84")
        || normalized.contains("WGS1984")
    {
        Ok(RasterCrs::Wgs84Degrees)
    } else if normalized.contains("EPSG4269")
        || normalized.contains("NAD83")
        || normalized.contains("NORTHAMERICANDATUM1983")
    {
        Ok(RasterCrs::Nad83Degrees)
    } else {
        Err(TrailgenError::UnsupportedFormat(
            "VRT SRS must declare WGS84/NAD83/CRS84 geographic, EPSG:3857 Web Mercator, or WGS84/NAD83 UTM (EPSG:326xx/327xx/269xx)"
                .to_owned(),
        ))
    }
}

fn required_attr<T>(node: roxmltree::Node<'_, '_>, key: &str) -> Result<T>
where
    T::Err: std::fmt::Display,
    T: FromStr,
{
    node.attribute(key)
        .ok_or_else(|| TrailgenError::InvalidData(format!("VRT missing {key} attribute")))?
        .parse()
        .map_err(|error| TrailgenError::InvalidData(format!("invalid VRT {key}: {error}")))
}

fn child_text<'a>(node: roxmltree::Node<'a, 'a>, tag: &str) -> Option<&'a str> {
    node.children()
        .find(|child| child.has_tag_name(tag))
        .and_then(|child| child.text())
}

fn ensure_identity_rect(
    simple: roxmltree::Node<'_, '_>,
    tag: &str,
    width: usize,
    height: usize,
) -> Result<()> {
    let Some(rect) = simple.children().find(|node| node.has_tag_name(tag)) else {
        return Ok(());
    };
    let x_off = required_attr::<usize>(rect, "xOff")?;
    let y_off = required_attr::<usize>(rect, "yOff")?;
    let x_size = required_attr::<usize>(rect, "xSize")?;
    let y_size = required_attr::<usize>(rect, "ySize")?;
    if x_off == 0 && y_off == 0 && x_size == width && y_size == height {
        Ok(())
    } else {
        Err(TrailgenError::UnsupportedFormat(format!(
            "VRT DEM only supports full-raster identity {tag}"
        )))
    }
}

fn resolve_vrt_source(
    vrt_path: &Path,
    source_filename: roxmltree::Node<'_, '_>,
    source_text: &str,
) -> PathBuf {
    let path = PathBuf::from(source_text);
    if path.is_absolute() {
        return path;
    }
    let relative = source_filename
        .attribute("relativeToVRT")
        .is_some_and(|value| matches!(value, "1" | "true" | "TRUE" | "True"));
    if relative || vrt_path.parent().is_some() {
        vrt_path
            .parent()
            .unwrap_or_else(|| Path::new("."))
            .join(path)
    } else {
        path
    }
}

#[derive(Clone, Copy)]
struct GeoTiffGeoref {
    crs: RasterCrs,
    origin_x: f64,
    origin_y: f64,
    pixel_width: f64,
    pixel_height: f64,
    transform: Option<RasterTransform>,
}

impl GeoTiffGeoref {
    fn read<R: std::io::Read + std::io::Seek>(decoder: &mut Decoder<R>) -> Result<Self> {
        let crs = validate_geokeys(decoder)?;
        if let Some(xs) = optional_f64_vec(decoder, Tag::ModelTransformationTag)? {
            let transform = RasterTransform::from_model_transformation(&xs)?;
            return Ok(Self {
                crs,
                origin_x: transform.x0,
                origin_y: transform.y0,
                pixel_width: transform.pixel_width(),
                pixel_height: transform.pixel_height(),
                transform: Some(transform),
            });
        }
        let scale = required_f64_vec(decoder, Tag::ModelPixelScaleTag)?;
        let tiepoint = required_f64_vec(decoder, Tag::ModelTiepointTag)?;
        if scale.len() < 2 || tiepoint.len() < 6 {
            return Err(TrailgenError::InvalidData(
                "GeoTIFF DEM requires ModelPixelScaleTag[0..2] and ModelTiepointTag[0..6]"
                    .to_owned(),
            ));
        }
        let pixel_width = scale[0].abs();
        let pixel_height = scale[1].abs();
        if !pixel_width.is_finite()
            || !pixel_height.is_finite()
            || pixel_width <= 0.0
            || pixel_height <= 0.0
        {
            return Err(TrailgenError::InvalidData(
                "GeoTIFF DEM pixel scale must be finite and positive".to_owned(),
            ));
        }
        let origin_x = tiepoint[0].mul_add(-pixel_width, tiepoint[3]);
        let origin_y = tiepoint[1].mul_add(pixel_height, tiepoint[4]);
        if !origin_x.is_finite() || !origin_y.is_finite() {
            return Err(TrailgenError::InvalidData(
                "GeoTIFF DEM tiepoint georeference must be finite".to_owned(),
            ));
        }
        Ok(Self {
            crs,
            origin_x,
            origin_y,
            pixel_width,
            pixel_height,
            transform: None,
        })
    }
}

fn validate_geokeys<R: std::io::Read + std::io::Seek>(
    decoder: &mut Decoder<R>,
) -> Result<RasterCrs> {
    let Some(keys) = optional_u16_vec(decoder, Tag::GeoKeyDirectoryTag)? else {
        return Err(TrailgenError::InvalidData(
            "GeoTIFF DEM requires GeoKeyDirectoryTag declaring WGS84/NAD83 geographic, EPSG:3857 Web Mercator, or WGS84/NAD83 UTM (EPSG:326xx/327xx/269xx)"
                .to_owned(),
        ));
    };
    if keys.len() < 4 {
        return Err(TrailgenError::InvalidData(
            "GeoTIFF GeoKeyDirectoryTag header is truncated".to_owned(),
        ));
    }
    let count = usize::from(keys[3]);
    if keys.len() < 4 + count.saturating_mul(4) {
        return Err(TrailgenError::InvalidData(
            "GeoTIFF GeoKeyDirectoryTag entries are truncated".to_owned(),
        ));
    }
    let mut model_type = None;
    let mut angular_units = None;
    let mut geographic_crs = None;
    let mut projected_crs = None;
    let mut linear_units = None;
    for entry in keys[4..][..count * 4].chunks_exact(4) {
        match entry[0] {
            1024 if entry[1] == 0 => model_type = Some(entry[3]),
            2048 if entry[1] == 0 => geographic_crs = Some(entry[3]),
            2054 if entry[1] == 0 => angular_units = Some(entry[3]),
            3072 if entry[1] == 0 => projected_crs = Some(entry[3]),
            3076 if entry[1] == 0 => linear_units = Some(entry[3]),
            _ => {}
        }
    }
    if model_type == Some(1) {
        if let Some(crs) = projected_crs.and_then(projected_raster_crs) {
            if linear_units.is_some_and(|unit| unit != 9001) {
                return Err(TrailgenError::UnsupportedFormat(
                    "projected GeoTIFF DEM linear units must be metres".to_owned(),
                ));
            }
            return Ok(crs);
        }
        return Err(TrailgenError::UnsupportedFormat(
            "projected GeoTIFF DEM must be EPSG:3857 Web Mercator or WGS84/NAD83 UTM (EPSG:326xx/327xx/269xx), or reprojected before ingestion"
                .to_owned(),
        ));
    }
    if model_type == Some(2) {
        if angular_units.is_some_and(|unit| unit != 9102) {
            return Err(TrailgenError::UnsupportedFormat(
                "GeoTIFF DEM angular units must be degrees".to_owned(),
            ));
        }
        return geographic_crs.map_or(Ok(RasterCrs::Wgs84Degrees), geographic_raster_crs);
    }
    Err(TrailgenError::UnsupportedFormat(
        "GeoTIFF DEM must declare GTModelTypeGeoKey=geographic, EPSG:3857 projected, or WGS84/NAD83 UTM projected".to_owned(),
    ))
}

fn geographic_raster_crs(epsg: u16) -> Result<RasterCrs> {
    match epsg {
        4326 => Ok(RasterCrs::Wgs84Degrees),
        4269 => Ok(RasterCrs::Nad83Degrees),
        _ => Err(TrailgenError::UnsupportedFormat(
            "GeoTIFF DEM geographic CRS must be WGS84 EPSG:4326 or NAD83 EPSG:4269".to_owned(),
        )),
    }
}

fn projected_raster_crs(epsg: u16) -> Option<RasterCrs> {
    match epsg {
        3857 => Some(RasterCrs::WebMercatorMeters),
        32601..=32760 | 26901..=26923 => Some(RasterCrs::UtmMeters(UtmCrs::from_epsg(epsg)?)),
        _ => None,
    }
}

fn decoding_result_to_f64(image: DecodingResult) -> Vec<f64> {
    match image {
        DecodingResult::U8(xs) => xs.into_iter().map(f64::from).collect(),
        DecodingResult::U16(xs) => xs.into_iter().map(f64::from).collect(),
        DecodingResult::U32(xs) => xs.into_iter().map(f64::from).collect(),
        DecodingResult::U64(xs) => xs.into_iter().map(|x| numeric_to_f64(&x)).collect(),
        DecodingResult::F16(xs) => xs.into_iter().map(|x| f64::from(x.to_f32())).collect(),
        DecodingResult::F32(xs) => xs.into_iter().map(f64::from).collect(),
        DecodingResult::F64(xs) => xs,
        DecodingResult::I8(xs) => xs.into_iter().map(f64::from).collect(),
        DecodingResult::I16(xs) => xs.into_iter().map(f64::from).collect(),
        DecodingResult::I32(xs) => xs.into_iter().map(f64::from).collect(),
        DecodingResult::I64(xs) => xs.into_iter().map(|x| numeric_to_f64(&x)).collect(),
    }
}

fn required_f64_vec<R: std::io::Read + std::io::Seek>(
    decoder: &mut Decoder<R>,
    tag: Tag,
) -> Result<Vec<f64>> {
    optional_f64_vec(decoder, tag)?
        .ok_or_else(|| TrailgenError::InvalidData(format!("GeoTIFF DEM missing {tag:?}")))
}

fn optional_f64_vec<R: std::io::Read + std::io::Seek>(
    decoder: &mut Decoder<R>,
    tag: Tag,
) -> Result<Option<Vec<f64>>> {
    match decoder.get_tag(tag) {
        Ok(value) => value
            .into_f64_vec()
            .map(Some)
            .map_err(|error| TrailgenError::InvalidData(format!("read GeoTIFF {tag:?}: {error}"))),
        Err(error) if missing_tag(&error, tag) => Ok(None),
        Err(error) => Err(TrailgenError::InvalidData(format!(
            "read GeoTIFF {tag:?}: {error}"
        ))),
    }
}

fn optional_u16_vec<R: std::io::Read + std::io::Seek>(
    decoder: &mut Decoder<R>,
    tag: Tag,
) -> Result<Option<Vec<u16>>> {
    match decoder.get_tag(tag) {
        Ok(value) => value
            .into_u16_vec()
            .map(Some)
            .map_err(|error| TrailgenError::InvalidData(format!("read GeoTIFF {tag:?}: {error}"))),
        Err(error) if missing_tag(&error, tag) => Ok(None),
        Err(error) => Err(TrailgenError::InvalidData(format!(
            "read GeoTIFF {tag:?}: {error}"
        ))),
    }
}

fn optional_ascii<R: std::io::Read + std::io::Seek>(
    decoder: &mut Decoder<R>,
    tag: Tag,
) -> Result<Option<String>> {
    match decoder.get_tag_ascii_string(tag) {
        Ok(value) => Ok(Some(value)),
        Err(error) if missing_tag(&error, tag) => Ok(None),
        Err(error) => Err(TrailgenError::InvalidData(format!(
            "read GeoTIFF {tag:?}: {error}"
        ))),
    }
}

fn missing_tag(error: &TiffError, tag: Tag) -> bool {
    matches!(
        error,
        TiffError::FormatError(TiffFormatError::RequiredTagNotFound(missing)) if *missing == tag
    )
}

fn header<'a, I, T>(lines: &mut I, key: &str) -> Result<T>
where
    I: Iterator<Item = &'a str>,
    T::Err: std::fmt::Display,
    T: FromStr,
{
    let line = lines
        .next()
        .ok_or_else(|| TrailgenError::InvalidData(format!("raster missing {key} header")))?;
    let (actual, value) = split_header(line)?;
    if !actual.eq_ignore_ascii_case(key) {
        return Err(TrailgenError::InvalidData(format!(
            "raster expected {key} header, found {actual}"
        )));
    }
    value
        .parse()
        .map_err(|e| TrailgenError::InvalidData(format!("invalid raster {key}: {e}")))
}

fn optional_header<'a, I, T>(lines: &mut Peekable<I>, key: &str) -> Result<Option<T>>
where
    I: Iterator<Item = &'a str>,
    T::Err: std::fmt::Display,
    T: FromStr,
{
    let Some(line) = lines.peek() else {
        return Ok(None);
    };
    let (actual, _) = split_header(line)?;
    if !actual.eq_ignore_ascii_case(key) {
        return Ok(None);
    }
    let Some(line) = lines.next() else {
        return Ok(None);
    };
    let (_, value) = split_header(line)?;
    value
        .parse()
        .map(Some)
        .map_err(|e| TrailgenError::InvalidData(format!("invalid raster {key}: {e}")))
}

fn split_header(line: &str) -> Result<(&str, &str)> {
    let mut parts = line.split_whitespace();
    let key = parts
        .next()
        .ok_or_else(|| TrailgenError::InvalidData("empty raster header".to_owned()))?;
    let value = parts
        .next()
        .ok_or_else(|| TrailgenError::InvalidData(format!("raster header {key} has no value")))?;
    Ok((key, value))
}

fn f64_to_index(value: f64, upper: usize) -> Option<usize> {
    if value.is_nan() || value < 0.0 || upper == 0 {
        return None;
    }
    let floored = value.floor().min(usize_to_f64(upper.saturating_sub(1)));
    floored.to_usize()
}

fn usize_to_f64(value: usize) -> f64 {
    value.to_f64().unwrap_or(f64::INFINITY)
}

fn numeric_to_f64<T: ToPrimitive + ?Sized>(value: &T) -> f64 {
    value.to_f64().unwrap_or(f64::NAN)
}