oxiproj-transformations 0.1.1

#![forbid(unsafe_code)]
//! XYZ grid shift (`xyzgridshift`) — port of PROJ `src/transformations/xyzgridshift.cpp`.
//! Uses a 3-band grid where:
//!   band 0 = eastward (longitude) shift in degrees
//!   band 1 = northward (latitude) shift in degrees
//!   band 2 = upward (ellipsoidal height) shift in meters
//! Input/output in radians (horizontal) and meters (vertical).

use crate::{TransBuild, TransParams};
use oxiproj_core::{Coord, IoUnits, Operation, ProjError, ProjResult};
use oxiproj_grids::{sample_grid, GridSet};

const MAX_ITER: usize = 10;
const ITER_TOL: f64 = 1e-12;

#[derive(Debug)]
struct XyzGridShift {
    grid: GridSet,
}

impl Operation for XyzGridShift {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let lon_deg = v[0].to_degrees();
        let lat_deg = v[1].to_degrees();
        let shifts = sample_grid(&self.grid, lat_deg, lon_deg).ok_or(ProjError::OutsideGrid)?;
        let dx_rad = shifts[0].to_radians();
        let dy_rad = shifts[1].to_radians();
        let dz = if shifts.len() >= 3 { shifts[2] } else { 0.0 };
        Ok(Coord::new(v[0] + dx_rad, v[1] + dy_rad, v[2] + dz, v[3]))
    }

    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let mut lon_r = v[0];
        let mut lat_r = v[1];
        // Iterative inversion for horizontal.
        for _ in 0..MAX_ITER {
            let lon_deg_try = lon_r.to_degrees();
            let lat_deg_try = lat_r.to_degrees();
            let shifts =
                sample_grid(&self.grid, lat_deg_try, lon_deg_try).ok_or(ProjError::OutsideGrid)?;
            let new_lon = v[0] - shifts[0].to_radians();
            let new_lat = v[1] - shifts[1].to_radians();
            let dlon = (new_lon - lon_r).abs();
            let dlat = (new_lat - lat_r).abs();
            lon_r = new_lon;
            lat_r = new_lat;
            if dlon < ITER_TOL && dlat < ITER_TOL {
                break;
            }
        }
        // Vertical: look up shift at the (converged) source position and subtract.
        let shifts_at_src = sample_grid(&self.grid, lat_r.to_degrees(), lon_r.to_degrees())
            .ok_or(ProjError::OutsideGrid)?;
        let dz = if shifts_at_src.len() >= 3 {
            shifts_at_src[2]
        } else {
            0.0
        };
        Ok(Coord::new(lon_r, lat_r, v[2] - dz, v[3]))
    }

    fn has_inverse(&self) -> bool {
        true
    }
}

/// Construct an `xyzgridshift` transform from parsed parameters.
///
/// Requires `+grids=<name>` and a populated [`GridRegistry`](crate::GridRegistry).
/// The grid must have at least 2 bands (dx, dy); 3 bands adds dz.
/// The grid is parsed using [`oxiproj_grids::read_geotiff`] (GeoTIFF format).
/// If the grid is not a valid GeoTIFF, construction fails.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let grid_name = p.params.get_str("grids").ok_or(ProjError::MissingArg)?;
    let registry = p.registry.ok_or(ProjError::FileNotFound)?;
    let data = registry
        .get_grid(grid_name)
        .ok_or(ProjError::FileNotFound)?;
    let grid = oxiproj_grids::read_geotiff(data, grid_name)?;
    if grid.bands.len() < 2 {
        return Err(ProjError::FileNotFound);
    }
    Ok(TransBuild::new(
        Box::new(XyzGridShift { grid }),
        IoUnits::Radians,
        IoUnits::Radians,
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use oxiproj_core::DEG_TO_RAD;
    use oxiproj_grids::{GridBand, GridExtent, GridSet, GridSource};

    /// Build a 2×2 GridSet with a given number of bands, each with a uniform value.
    fn make_uniform_xyz_gridset(values: &[f32]) -> GridSet {
        let extent = GridExtent {
            ll_lat: 0.0,
            ll_lon: 0.0,
            ur_lat: 1.0,
            ur_lon: 1.0,
            lat_inc: 1.0,
            lon_inc: 1.0,
        };
        let bands = values
            .iter()
            .map(|&v| GridBand {
                extent: extent.clone(),
                values: vec![v; 4],
            })
            .collect();
        GridSet {
            bands,
            source: GridSource::Memory,
        }
    }

    #[test]
    fn test_xyzgridshift_forward_2band() {
        let gs = make_uniform_xyz_gridset(&[0.1f32, 0.2]);
        let op = XyzGridShift { grid: gs };
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 100.0, 0.0);
        let out = op.forward_4d(input).unwrap();
        let v = out.v();
        let expected_lon = (0.5 + 0.1) * DEG_TO_RAD;
        let expected_lat = (0.5 + 0.2) * DEG_TO_RAD;
        assert!((v[0] - expected_lon).abs() < 1e-9, "lon");
        assert!((v[1] - expected_lat).abs() < 1e-9, "lat");
        assert!((v[2] - 100.0).abs() < 1e-6, "z unchanged for 2-band");
    }

    #[test]
    fn test_xyzgridshift_forward_3band() {
        let gs = make_uniform_xyz_gridset(&[0.1f32, 0.2, 5.0]);
        let op = XyzGridShift { grid: gs };
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 100.0, 0.0);
        let out = op.forward_4d(input).unwrap();
        assert!((out.v()[2] - 105.0).abs() < 1e-6, "z += 5.0");
    }

    #[test]
    fn test_xyzgridshift_round_trip() {
        let gs = make_uniform_xyz_gridset(&[0.05f32, 0.03, 10.0]);
        let op = XyzGridShift { grid: gs };
        let input = Coord::new(0.4 * DEG_TO_RAD, 0.6 * DEG_TO_RAD, 50.0, 2020.0);
        let fwd = op.forward_4d(input).unwrap();
        let inv = op.inverse_4d(fwd).unwrap();
        let vi = inv.v();
        let vi0 = input.v();
        assert!((vi[0] - vi0[0]).abs() < 1e-9, "lon round-trip");
        assert!((vi[1] - vi0[1]).abs() < 1e-9, "lat round-trip");
        assert!((vi[2] - vi0[2]).abs() < 1e-6, "z round-trip");
    }

    #[test]
    fn test_xyzgridshift_outside_err() {
        let gs = make_uniform_xyz_gridset(&[0.1f32, 0.2]);
        let op = XyzGridShift { grid: gs };
        let far = Coord::new(45.0 * DEG_TO_RAD, 45.0 * DEG_TO_RAD, 0.0, 0.0);
        assert_eq!(op.forward_4d(far).err(), Some(ProjError::OutsideGrid));
    }
}