oxiproj-transformations 0.1.2

#![forbid(unsafe_code)]
//! XYZ grid shift (`xyzgridshift`) — port of PROJ `src/transformations/xyzgridshift.cpp`.
//!
//! This operation performs a *geocentric* grid shift. Both the input and the
//! output are geocentric Cartesian coordinates `X, Y, Z` expressed in **metres**
//! (`left = right = CARTESIAN`, matching PROJ). The 3-band grid stores the
//! translation to apply, expressed as `x_translation` / `y_translation` /
//! `z_translation` in metres (default band order `0`/`1`/`2`).
//!
//! For each point the Cartesian coordinate is converted to geodetic
//! `lon`/`lat` (via the operation's ellipsoid) purely to locate the grid cell;
//! the three interpolated offsets `dX, dY, dZ` are then added **directly to the
//! Cartesian coordinate** (never to `lon`/`lat`/`h`). An optional `+multiplier`
//! scales the offsets.
//!
//! Two conventions are supported through `+grid_ref`:
//!   * `input_crs` (default): the grid is referenced in the *input* CRS, so the
//!     forward direction is a single direct lookup and the inverse is solved by
//!     a 10-iteration fixed-point.
//!   * `output_crs`: the grid is referenced in the *output* CRS (e.g. the
//!     classic NTF→RGF93 grid, whose deltas are referenced in RGF93). The roles
//!     of direct/iterative are swapped.

use crate::internal::cartesian_to_geodetic;
use crate::{TransBuild, TransParams};
use oxiproj_core::{Ellipsoid, IoUnits, Lpz, Operation, ProjError, ProjResult, Xyz};
use oxiproj_grids::{sample_grid, BandUnit, GridSet};

/// Maximum fixed-point iterations for the iterative adjustment, per PROJ
/// (`iterative_adjustment` loops `for i in 0..10`).
const MAX_ITER: usize = 10;
/// Convergence threshold on the squared residual, in metres² — matches PROJ's
/// `if (err < 1e-10) break;`.
const ITER_TOL: f64 = 1e-10;

/// State for the `xyzgridshift` transformation.
///
/// Mirrors PROJ's `xyzgridshiftData`: the resolved grid, the ellipsoid used to
/// convert Cartesian → geodetic for the grid lookup (PROJ builds an internal
/// `+proj=cart` sub-transform that inherits `P`'s ellipsoid), the `grid_ref`
/// convention flag, and the `+multiplier` applied to the deltas.
#[derive(Debug)]
struct XyzGridShift {
    grid: GridSet,
    ell: Ellipsoid,
    grid_ref_is_input: bool,
    multiplier: f64,
}

impl XyzGridShift {
    /// Look up the `(dX, dY, dZ)` translation (metres) for the geocentric point
    /// `(x, y, z)`. Ported from PROJ `get_grid_values`: the Cartesian point is
    /// converted to geodetic `lon`/`lat` (the grid cell locator), the three
    /// bands are bilinearly interpolated in default `x/y/z_translation` order
    /// (`0`/`1`/`2`), and the `+multiplier` is applied.
    fn grid_values(&self, x: f64, y: f64, z: f64) -> ProjResult<(f64, f64, f64)> {
        let geod = cartesian_to_geodetic(Xyz::new(x, y, z), &self.ell);
        let lon_deg = geod.lam.to_degrees();
        let lat_deg = geod.phi.to_degrees();
        let s = sample_grid(&self.grid, lat_deg, lon_deg).ok_or(ProjError::OutsideGrid)?;
        // A 3-band grid is guaranteed by `new()`; `sample_grid` returns one
        // value per band in band order (default x/y/z_translation = 0/1/2).
        let dx = s.first().copied().ok_or(ProjError::OutsideGrid)?;
        let dy = s.get(1).copied().ok_or(ProjError::OutsideGrid)?;
        let dz = s.get(2).copied().ok_or(ProjError::OutsideGrid)?;
        Ok((
            dx * self.multiplier,
            dy * self.multiplier,
            dz * self.multiplier,
        ))
    }

    /// Direct adjustment: one grid lookup, offsets added to the Cartesian
    /// coordinate scaled by `factor` (`+1` forward, `-1` inverse). Ported from
    /// PROJ `direct_adjustment`.
    fn direct_adjustment(
        &self,
        x: f64,
        y: f64,
        z: f64,
        factor: f64,
    ) -> ProjResult<(f64, f64, f64)> {
        let (dx, dy, dz) = self.grid_values(x, y, z)?;
        Ok((x + factor * dx, y + factor * dy, z + factor * dz))
    }

    /// Iterative (fixed-point) adjustment used when the grid is referenced in
    /// the CRS *opposite* to the direction being solved. Ported from PROJ
    /// `iterative_adjustment` (10 iterations, squared-residual break `< 1e-10`).
    fn iterative_adjustment(
        &self,
        x0: f64,
        y0: f64,
        z0: f64,
        factor: f64,
    ) -> ProjResult<(f64, f64, f64)> {
        let (mut x, mut y, mut z) = (x0, y0, z0);
        for _ in 0..MAX_ITER {
            let (mut dx, mut dy, mut dz) = self.grid_values(x, y, z)?;
            dx *= factor;
            dy *= factor;
            dz *= factor;

            let ex = (x - x0) - dx;
            let ey = (y - y0) - dy;
            let ez = (z - z0) - dz;
            let err = ex * ex + ey * ey + ez * ez;

            x = x0 + dx;
            y = y0 + dy;
            z = z0 + dz;

            if err < ITER_TOL {
                break;
            }
        }
        Ok((x, y, z))
    }
}

impl Operation for XyzGridShift {
    /// Forward: geocentric Cartesian `(X, Y, Z)` in → shifted `(X, Y, Z)` out.
    ///
    /// With `left = CARTESIAN` the geocentric coordinates arrive in the `Lpz`
    /// slots (`lam`, `phi`, `z`), exactly as for `topocentric`. Ported from
    /// PROJ `pj_xyzgridshift_forward_3d`.
    fn forward_3d(&self, lpz: Lpz) -> ProjResult<Xyz> {
        let (x, y, z) = (lpz.lam, lpz.phi, lpz.z);
        let (rx, ry, rz) = if self.grid_ref_is_input {
            self.direct_adjustment(x, y, z, 1.0)?
        } else {
            self.iterative_adjustment(x, y, z, 1.0)?
        };
        Ok(Xyz::new(rx, ry, rz))
    }

    /// Inverse: shifted geocentric Cartesian `(X, Y, Z)` in → original
    /// `(X, Y, Z)` out (returned in the `Lpz` slots to preserve the `CARTESIAN`
    /// convention). Ported from PROJ `pj_xyzgridshift_reverse_3d`.
    fn inverse_3d(&self, xyz: Xyz) -> ProjResult<Lpz> {
        let (x, y, z) = (xyz.x, xyz.y, xyz.z);
        let (rx, ry, rz) = if self.grid_ref_is_input {
            self.iterative_adjustment(x, y, z, -1.0)?
        } else {
            self.direct_adjustment(x, y, z, -1.0)?
        };
        Ok(Lpz::new(rx, ry, rz))
    }

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

/// Construct an `xyzgridshift` transform from parsed parameters.
///
/// Requires `+grids=<name>`. The named grid is resolved registry-first, then
/// from the local `PROJ_DATA`/`PROJ_LIB` directories and the unified disk
/// cache, then (with the `network` feature) the PROJ CDN — see
/// `super::gridshift::resolve_grid_bytes`. The grid must be a GeoTIFF with at
/// least 3 bands (`x_translation`, `y_translation`, `z_translation`) whose unit
/// is metre; it is parsed with [`oxiproj_grids::read_geotiff`].
///
/// Optional parameters, mirroring PROJ:
///   * `+grid_ref=input_crs` (default) or `+grid_ref=output_crs`.
///   * `+multiplier=<f>` scales the deltas (default `1.0`).
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let grid_name = p.params.get_str("grids").ok_or(ProjError::MissingArg)?;

    // `+grid_ref`: input_crs (default) vs output_crs.
    let grid_ref_is_input = match p.params.get_str("grid_ref") {
        None | Some("input_crs") => true,
        Some("output_crs") => false,
        Some(_) => return Err(ProjError::IllegalArgValue),
    };

    // `+multiplier`, default 1.0.
    let multiplier = p.params.get_f64("multiplier").unwrap_or(1.0);

    let data = super::gridshift::resolve_grid_bytes(p.registry, grid_name)?;
    let grid = oxiproj_grids::read_geotiff(&data, grid_name)?;

    // PROJ: "xyzgridshift: grid has not enough samples" if samplesPerPixel < 3.
    if grid.bands.len() < 3 {
        return Err(ProjError::IllegalArgValue);
    }

    // PROJ: "Only unit=metre currently handled" — an empty/unknown unit is
    // accepted (real translation grids frequently omit the tag).
    if let Some(band0) = grid.bands.first() {
        match band0.semantics.unit {
            BandUnit::ArcSecond | BandUnit::Degree | BandUnit::Radian => {
                return Err(ProjError::IllegalArgValue);
            }
            BandUnit::Metre | BandUnit::Unknown => {}
        }
    }

    Ok(TransBuild::new(
        Box::new(XyzGridShift {
            grid,
            ell: *p.ellipsoid,
            grid_ref_is_input,
            multiplier,
        }),
        IoUnits::Cartesian,
        IoUnits::Cartesian,
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::internal::geodetic_to_cartesian;
    use oxiproj_core::{Coord, Ellipsoid, DEG_TO_RAD};
    use oxiproj_grids::{GridBand, GridExtent, GridSource};

    fn wgs84() -> Ellipsoid {
        Ellipsoid::named("WGS84").expect("WGS84 ellipsoid")
    }

    /// Build a 2×2 GridSet spanning lon/lat `0..1` degrees with the given
    /// per-band uniform values (one entry per band).
    fn make_uniform_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],
                semantics: Default::default(),
            })
            .collect();
        GridSet {
            bands,
            source: GridSource::Memory,
        }
    }

    fn op(values: &[f32], grid_ref_is_input: bool, multiplier: f64) -> XyzGridShift {
        XyzGridShift {
            grid: make_uniform_gridset(values),
            ell: wgs84(),
            grid_ref_is_input,
            multiplier,
        }
    }

    /// Geocentric coordinate for a point inside the test grid (lon/lat ≈ 0.5°).
    fn point_in_grid() -> Xyz {
        geodetic_to_cartesian(
            Lpz::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 100.0),
            &wgs84(),
        )
    }

    #[test]
    fn forward_adds_deltas_to_cartesian() {
        let o = op(&[1000.0, 2000.0, 3000.0], true, 1.0);
        let p = point_in_grid();
        let out = o.forward_3d(Lpz::new(p.x, p.y, p.z)).expect("forward");
        assert!((out.x - (p.x + 1000.0)).abs() < 1e-6, "X = {}", out.x);
        assert!((out.y - (p.y + 2000.0)).abs() < 1e-6, "Y = {}", out.y);
        assert!((out.z - (p.z + 3000.0)).abs() < 1e-6, "Z = {}", out.z);
    }

    #[test]
    fn multiplier_scales_deltas() {
        let o = op(&[1000.0, 2000.0, 3000.0], true, 2.0);
        let p = point_in_grid();
        let out = o.forward_3d(Lpz::new(p.x, p.y, p.z)).expect("forward");
        assert!((out.x - (p.x + 2000.0)).abs() < 1e-6, "X = {}", out.x);
        assert!((out.y - (p.y + 4000.0)).abs() < 1e-6, "Y = {}", out.y);
        assert!((out.z - (p.z + 6000.0)).abs() < 1e-6, "Z = {}", out.z);
    }

    #[test]
    fn round_trip_input_crs() {
        let o = op(&[123.0, -456.0, 789.0], true, 1.0);
        let p = point_in_grid();
        let fwd = o.forward_3d(Lpz::new(p.x, p.y, p.z)).expect("forward");
        let inv = o.inverse_3d(fwd).expect("inverse");
        assert!((inv.lam - p.x).abs() < 1e-6, "X = {}", inv.lam);
        assert!((inv.phi - p.y).abs() < 1e-6, "Y = {}", inv.phi);
        assert!((inv.z - p.z).abs() < 1e-6, "Z = {}", inv.z);
    }

    #[test]
    fn round_trip_output_crs() {
        let o = op(&[123.0, -456.0, 789.0], false, 1.0);
        let p = point_in_grid();
        let fwd = o.forward_3d(Lpz::new(p.x, p.y, p.z)).expect("forward");
        let inv = o.inverse_3d(fwd).expect("inverse");
        assert!((inv.lam - p.x).abs() < 1e-6, "X = {}", inv.lam);
        assert!((inv.phi - p.y).abs() < 1e-6, "Y = {}", inv.phi);
        assert!((inv.z - p.z).abs() < 1e-6, "Z = {}", inv.z);
    }

    #[test]
    fn input_and_output_crs_match_for_uniform_grid() {
        // For a spatially uniform grid the direct and iterative adjustments
        // yield the same result, so the two grid_ref conventions agree.
        let p = point_in_grid();
        let a = op(&[10.0, 20.0, 30.0], true, 1.0)
            .forward_3d(Lpz::new(p.x, p.y, p.z))
            .expect("input_crs");
        let b = op(&[10.0, 20.0, 30.0], false, 1.0)
            .forward_3d(Lpz::new(p.x, p.y, p.z))
            .expect("output_crs");
        assert!((a.x - b.x).abs() < 1e-6);
        assert!((a.y - b.y).abs() < 1e-6);
        assert!((a.z - b.z).abs() < 1e-6);
    }

    #[test]
    fn forward_4d_passes_time_through() {
        let o = op(&[1.0, 2.0, 3.0], true, 1.0);
        let p = point_in_grid();
        let out = o
            .forward_4d(Coord::new(p.x, p.y, p.z, 2025.5))
            .expect("forward_4d");
        let v = out.v();
        assert!((v[3] - 2025.5).abs() < 1e-9, "t = {}", v[3]);
        assert!((v[0] - (p.x + 1.0)).abs() < 1e-6, "X = {}", v[0]);
    }

    #[test]
    fn outside_grid_errors() {
        let o = op(&[1.0, 2.0, 3.0], true, 1.0);
        // A geocentric point whose geodetic lon/lat (≈ 45°/45°) is outside the
        // 0..1° test grid.
        let far = geodetic_to_cartesian(
            Lpz::new(45.0 * DEG_TO_RAD, 45.0 * DEG_TO_RAD, 0.0),
            &wgs84(),
        );
        assert_eq!(
            o.forward_3d(Lpz::new(far.x, far.y, far.z)).err(),
            Some(ProjError::OutsideGrid)
        );
    }

    #[test]
    fn missing_grids_errors() {
        use crate::TransParamLookup;
        struct Empty;
        impl TransParamLookup for Empty {
            fn get_dms(&self, _: &str) -> Option<f64> {
                None
            }
            fn get_f64(&self, _: &str) -> Option<f64> {
                None
            }
            fn get_int(&self, _: &str) -> Option<i64> {
                None
            }
            fn get_str(&self, _: &str) -> Option<&str> {
                None
            }
            fn get_bool(&self, _: &str) -> bool {
                false
            }
            fn exists(&self, _: &str) -> bool {
                false
            }
        }
        let ell = wgs84();
        let err = new(&TransParams {
            ellipsoid: &ell,
            params: &Empty,
            registry: None,
        })
        .err();
        assert_eq!(err, Some(ProjError::MissingArg));
    }
}