oxiproj-transformations 0.1.2

#![forbid(unsafe_code)]
//! Kinematic velocity-grid deformation (`deformation`) — port of PROJ
//! `src/transformations/deformation.cpp`.
//!
//! Performs a datum shift by means of a deformation / velocity model:
//!
//! ```text
//! X_out = X_in + (t_obs − t_epoch) · DX
//! Y_out = Y_in + (t_obs − t_epoch) · DY
//! Z_out = Z_in + (t_obs − t_epoch) · DZ
//! ```
//!
//! The operation takes **geocentric Cartesian** coordinates as input and
//! returns geocentric Cartesian coordinates as well (`left = right =
//! CARTESIAN`, matching PROJ). The gridded corrections are stored in
//! **east / north / up (ENU)** space, in units of **millimetres per year**.
//! For each point the Cartesian coordinate is converted to geodetic
//! `lon`/`lat` (purely to locate the grid cell), the three ENU velocities are
//! bilinearly interpolated, divided by 1000 to obtain metres per year, rotated
//! from ENU into geocentric `PJ_XYZ` space, scaled by the time span `dt`, and
//! added to the Cartesian coordinate.
//!
//! The full model is a 3-band Geodetic TIFF Grid whose bands carry the GDAL
//! metadata descriptions `east_velocity` / `north_velocity` / `up_velocity`
//! and the unit `millimetres per year`.
//!
//! Time-span parameters (mirroring PROJ, plus the historical `+t_obs`):
//!   * `+dt=<f64>` — a fixed time span (years) applied to every coordinate.
//!   * `+t_epoch=<f64>` — reference epoch; the span is computed per-coordinate
//!     as `dt = t − t_epoch`, where `t` is the coordinate's 4th (time) slot.
//!   * `+t_obs=<f64>` — a fixed observation epoch used together with
//!     `+t_epoch` to form a fixed span `dt = t_obs − t_epoch`.
//!
//! `+dt` and `+t_epoch` are mutually exclusive; at least one time-span source
//! must be supplied.

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

/// Maximum fixed-point iterations for the reverse shift, matching PROJ's
/// `#define MAX_ITERATIONS 10`.
const MAX_ITERATIONS: usize = 10;
/// Horizontal convergence threshold (metres) for the reverse shift, matching
/// PROJ's `#define TOL 1e-8`.
const TOL: f64 = 1e-8;

/// State for the `deformation` transformation.
///
/// Mirrors PROJ's `deformationData`: the resolved velocity 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
/// resolved E/N/U band indices, and the time-span source (`+dt`/`+t_obs`
/// yield a fixed span; `+t_epoch` yields a per-coordinate span).
#[derive(Debug)]
struct Deformation {
    grid: GridSet,
    ell: Ellipsoid,
    /// Band index carrying the eastward velocity (default 0).
    sample_e: usize,
    /// Band index carrying the northward velocity (default 1).
    sample_n: usize,
    /// Band index carrying the upward velocity (default 2).
    sample_u: usize,
    /// Fixed time span (years) from `+dt` or `+t_obs − t_epoch`; when set it
    /// takes precedence and the coordinate's time slot is ignored.
    fixed_dt: Option<f64>,
    /// Reference epoch from `+t_epoch`; when `fixed_dt` is `None` the span is
    /// computed per-coordinate as `t − t_epoch`.
    t_epoch: Option<f64>,
}

impl Deformation {
    /// Resolve the time span `dt` (years) for a coordinate whose time slot is
    /// `t`. A fixed `+dt`/`+t_obs` span always wins; otherwise the span is
    /// `t − t_epoch`, which requires a finite `t`.
    fn dt_for(&self, t: f64) -> ProjResult<f64> {
        if let Some(dt) = self.fixed_dt {
            return Ok(dt);
        }
        match self.t_epoch {
            Some(epoch) if t.is_finite() => Ok(t - epoch),
            Some(_) => Err(ProjError::MissingTime),
            None => Err(ProjError::MissingArg),
        }
    }

    /// Read the ENU velocity correction at the geocentric point `(x, y, z)` and
    /// return it rotated into geocentric `PJ_XYZ` space, in metres per year.
    ///
    /// Ported from PROJ `pj_deformation_get_grid_shift` /
    /// `pj_deformation_get_grid_values`: convert Cartesian → geodetic to locate
    /// the cell, bilinearly interpolate the three ENU bands (mm/yr), divide by
    /// 1000 to obtain m/yr, then rotate ENU → `PJ_XYZ` via the explicit
    /// sine/cosine matrix (Nørbech et al., 2003).
    fn grid_shift(&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)?;

        // ENU velocities, mm/yr → m/yr (PROJ divides by 1000).
        let e = s
            .get(self.sample_e)
            .copied()
            .ok_or(ProjError::OutsideGrid)?
            / 1000.0;
        let n = s
            .get(self.sample_n)
            .copied()
            .ok_or(ProjError::OutsideGrid)?
            / 1000.0;
        let u = s
            .get(self.sample_u)
            .copied()
            .ok_or(ProjError::OutsideGrid)?
            / 1000.0;

        // Pre-calc sines/cosines of the geodetic lon/lat.
        let sp = geod.phi.sin();
        let cp = geod.phi.cos();
        let sl = geod.lam.sin();
        let cl = geod.lam.cos();

        // ENU → PJ_XYZ.
        let tx = -sp * cl * n - sl * e + cp * cl * u;
        let ty = -sp * sl * n + cl * e + cp * sl * u;
        let tz = cp * n + sp * u;
        Ok((tx, ty, tz))
    }

    /// Forward correction: add `dt · shift` to the geocentric coordinate.
    /// Ported from PROJ `pj_deformation_forward_3d` / `_forward_4d`.
    fn forward_xyz(&self, x: f64, y: f64, z: f64, dt: f64) -> ProjResult<(f64, f64, f64)> {
        let (sx, sy, sz) = self.grid_shift(x, y, z)?;
        Ok((x + dt * sx, y + dt * sy, z + dt * sz))
    }

    /// Reverse correction: iteratively recover the coordinate that maps to the
    /// given (already-shifted) geocentric point. Ported verbatim from PROJ
    /// `pj_deformation_reverse_shift` (10 iterations, horizontal break at
    /// `hypot(dif.x, dif.y) <= 1e-8`, vertical component restored from the
    /// original grid shift `z0`).
    fn reverse_xyz(&self, ix: f64, iy: f64, iz: f64, dt: f64) -> ProjResult<(f64, f64, f64)> {
        let (dx, dy, dz) = self.grid_shift(ix, iy, iz)?;

        // Original z shift, re-applied after the horizontal loop converges.
        let z0 = dz;

        // The z-component is carried along the iteration (needed for the
        // Cartesian → geodetic conversion) and overwritten with z0 afterwards.
        let mut ox = ix - dt * dx;
        let mut oy = iy - dt * dy;
        let mut oz = iz + dt * dz;

        for _ in 0..MAX_ITERATIONS {
            let (ex, ey, ez) = match self.grid_shift(ox, oy, oz) {
                Ok(v) => v,
                // PROJ breaks out of the loop when the shift becomes invalid,
                // keeping the last iterate.
                Err(_) => break,
            };
            let difx = ox + dt * ex - ix;
            let dify = oy + dt * ey - iy;
            let difz = oz - dt * ez - iz;
            ox += difx;
            oy += dify;
            oz += difz;
            if difx.hypot(dify) <= TOL {
                break;
            }
        }

        oz = iz - dt * z0;
        Ok((ox, oy, oz))
    }
}

impl Operation for Deformation {
    /// Forward, 3D: geocentric Cartesian `(X, Y, Z)` in → shifted out.
    ///
    /// With `left = CARTESIAN` the geocentric coordinates arrive in the `Lpz`
    /// slots (`lam`, `phi`, `z`), matching `xyzgridshift`/`topocentric`. The 3D
    /// entry point requires a fixed span (`+dt`/`+t_obs`), mirroring PROJ's
    /// "`+dt` must be specified" guard.
    fn forward_3d(&self, lpz: Lpz) -> ProjResult<Xyz> {
        let dt = self.fixed_dt.ok_or(ProjError::MissingTime)?;
        let (x, y, z) = self.forward_xyz(lpz.lam, lpz.phi, lpz.z, dt)?;
        Ok(Xyz::new(x, y, z))
    }

    /// Inverse, 3D: shifted geocentric Cartesian `(X, Y, Z)` in → original out
    /// (returned in the `Lpz` slots to preserve the `CARTESIAN` convention).
    fn inverse_3d(&self, xyz: Xyz) -> ProjResult<Lpz> {
        let dt = self.fixed_dt.ok_or(ProjError::MissingTime)?;
        let (x, y, z) = self.reverse_xyz(xyz.x, xyz.y, xyz.z, dt)?;
        Ok(Lpz::new(x, y, z))
    }

    /// Forward, 4D: geocentric Cartesian `(X, Y, Z, t)`; the span is the fixed
    /// `+dt` if present, else `t − t_epoch`. Ported from
    /// `pj_deformation_forward_4d`.
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let dt = self.dt_for(v[3])?;
        let (x, y, z) = self.forward_xyz(v[0], v[1], v[2], dt)?;
        Ok(Coord::new(x, y, z, v[3]))
    }

    /// Inverse, 4D. Ported from `pj_deformation_reverse_4d`.
    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let dt = self.dt_for(v[3])?;
        let (x, y, z) = self.reverse_xyz(v[0], v[1], v[2], dt)?;
        Ok(Coord::new(x, y, z, v[3]))
    }

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

/// Construct a `deformation` transform from parsed parameters.
///
/// Requires `+grids=<name>` — a 3-band GeoTIFF velocity grid 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 have at least
/// three bands; the E/N/U band roles are recovered from the GDAL metadata
/// `DESCRIPTION` items (`east_velocity`/`north_velocity`/`up_velocity`,
/// defaulting to bands `0`/`1`/`2`), and the eastward band's unit — when
/// declared — must be `millimetres per year`.
///
/// Requires exactly one time-span source:
///   * `+dt=<f64>` (fixed span), or
///   * `+t_epoch=<f64>` (per-coordinate span `t − t_epoch`), or
///   * `+t_obs=<f64>` together with `+t_epoch` (fixed span `t_obs − t_epoch`).
///
/// `+dt` and `+t_epoch` are mutually exclusive.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let grid_name = p.params.get_str("grids").ok_or(ProjError::MissingArg)?;

    // Resolve the time-span source before touching the grid, mirroring PROJ's
    // parameter validation ordering.
    let dt = p.params.get_f64("dt");
    let t_epoch = p.params.get_f64("t_epoch");
    let t_obs = p.params.get_f64("t_obs");

    let (fixed_dt, epoch) = if let Some(dt_val) = dt {
        // PROJ: "+dt or +t_epoch are mutually exclusive."
        if t_epoch.is_some() {
            return Err(ProjError::MutuallyExclusiveArgs);
        }
        (Some(dt_val), None)
    } else if let Some(t_obs_val) = t_obs {
        // Historical `+t_obs`: a fixed observation epoch, meaningful only with
        // a reference epoch, yielding the fixed span `t_obs − t_epoch`.
        let epoch = t_epoch.ok_or(ProjError::MissingArg)?;
        (Some(t_obs_val - epoch), None)
    } else if let Some(epoch) = t_epoch {
        (None, Some(epoch))
    } else {
        // PROJ: "either +dt or +t_epoch needs to be set."
        return Err(ProjError::MissingArg);
    };

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

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

    // Resolve the E/N/U band roles from GDAL metadata, defaulting to 0/1/2, and
    // validate the eastward band's unit (PROJ: "Only unit=millimetres per year
    // currently handled" — an absent/empty unit is accepted).
    let mut sample_e = 0usize;
    let mut sample_n = 1usize;
    let mut sample_u = 2usize;
    if let Some(md) = read_geotiff_gdal_metadata(&data)? {
        for i in 0..grid.bands.len() {
            match md.metadata_item("DESCRIPTION", Some(i as i32)) {
                Some("east_velocity") => sample_e = i,
                Some("north_velocity") => sample_n = i,
                Some("up_velocity") => sample_u = i,
                _ => {}
            }
        }
        if let Some(unit) = md.metadata_item("UNITTYPE", Some(sample_e as i32)) {
            if !unit.is_empty() && unit != "millimetres per year" {
                return Err(ProjError::IllegalArgValue);
            }
        }
    }

    Ok(TransBuild::new(
        Box::new(Deformation {
            grid,
            ell: *p.ellipsoid,
            sample_e,
            sample_n,
            sample_u,
            fixed_dt,
            t_epoch: epoch,
        }),
        IoUnits::Cartesian,
        IoUnits::Cartesian,
    ))
}

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

    fn grs80() -> Ellipsoid {
        Ellipsoid::named("GRS80").expect("GRS80 ellipsoid")
    }

    /// Build a 2×2 GridSet spanning lon/lat `−1..1` degrees (so it covers the
    /// prime-meridian/equator point) with the given per-band uniform velocity
    /// values (mm/yr), one entry per band.
    fn make_velocity_grid(band_values: &[f32]) -> GridSet {
        let extent = GridExtent {
            ll_lat: -1.0,
            ll_lon: -1.0,
            ur_lat: 1.0,
            ur_lon: 1.0,
            lat_inc: 2.0,
            lon_inc: 2.0,
        };
        let bands = band_values
            .iter()
            .map(|&v| GridBand {
                extent: extent.clone(),
                values: vec![v; 4],
                semantics: Default::default(),
            })
            .collect();
        GridSet {
            bands,
            source: GridSource::Memory,
        }
    }

    /// Deformation with default band order (E=0, N=1, U=2) and a fixed span.
    fn op_fixed(band_values: &[f32], dt: f64) -> Deformation {
        Deformation {
            grid: make_velocity_grid(band_values),
            ell: grs80(),
            sample_e: 0,
            sample_n: 1,
            sample_u: 2,
            fixed_dt: Some(dt),
            t_epoch: None,
        }
    }

    #[test]
    fn enu_rotation_at_prime_meridian_equator() {
        // At lon=0, lat=0 the ENU→XYZ rotation reduces to:
        //   shift.x = u, shift.y = e, shift.z = n
        // Velocities: e=100, n=200, u=300 mm/yr → 0.1/0.2/0.3 m/yr.
        // With dt=10: dX=3.0 (u), dY=1.0 (e), dZ=2.0 (n).
        let op = op_fixed(&[100.0, 200.0, 300.0], 10.0);
        let p = geodetic_to_cartesian(Lpz::new(0.0, 0.0, 0.0), &grs80());
        let out = op.forward_3d(Lpz::new(p.x, p.y, p.z)).expect("forward");
        assert!((out.x - (p.x + 3.0)).abs() < 1e-9, "dX(u): {}", out.x - p.x);
        assert!((out.y - (p.y + 1.0)).abs() < 1e-9, "dY(e): {}", out.y - p.y);
        assert!((out.z - (p.z + 2.0)).abs() < 1e-9, "dZ(n): {}", out.z - p.z);
    }

    #[test]
    fn forward_round_trips_through_inverse() {
        // The reverse is PROJ's fixed-point iteration (horizontal break at
        // `hypot(dif) <= 1e-8`, vertical restored from the first shift). Because
        // the ENU→XYZ rotation makes the geocentric shift position-dependent
        // even for a spatially uniform ENU grid, the round-trip is not exact;
        // sub-millimetre closure matches PROJ's own reverse accuracy.
        let op = op_fixed(&[12.0, -34.0, 56.0], 25.0);
        let p = geodetic_to_cartesian(
            Lpz::new(0.3 * DEG_TO_RAD, 0.4 * DEG_TO_RAD, 100.0),
            &grs80(),
        );
        let fwd = op.forward_3d(Lpz::new(p.x, p.y, p.z)).expect("forward");
        let inv = op
            .inverse_3d(Xyz::new(fwd.x, fwd.y, fwd.z))
            .expect("inverse");
        assert!((inv.lam - p.x).abs() < 1e-3, "X {} vs {}", inv.lam, p.x);
        assert!((inv.phi - p.y).abs() < 1e-3, "Y {} vs {}", inv.phi, p.y);
        assert!((inv.z - p.z).abs() < 1e-3, "Z {} vs {}", inv.z, p.z);
    }

    #[test]
    fn zero_dt_is_identity() {
        let op = op_fixed(&[1000.0, 1000.0, 1000.0], 0.0);
        let p = geodetic_to_cartesian(Lpz::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0), &grs80());
        let out = op.forward_3d(Lpz::new(p.x, p.y, p.z)).expect("forward");
        assert!((out.x - p.x).abs() < 1e-12, "X");
        assert!((out.y - p.y).abs() < 1e-12, "Y");
        assert!((out.z - p.z).abs() < 1e-12, "Z");
    }

    #[test]
    fn forward_4d_uses_epoch_to_form_span() {
        // t_epoch mode: dt = t − t_epoch. At t=2010, epoch=2000 → dt=10, so the
        // result matches the fixed-dt=10 forward.
        let epoch_op = Deformation {
            grid: make_velocity_grid(&[100.0, 200.0, 300.0]),
            ell: grs80(),
            sample_e: 0,
            sample_n: 1,
            sample_u: 2,
            fixed_dt: None,
            t_epoch: Some(2000.0),
        };
        let p = geodetic_to_cartesian(Lpz::new(0.0, 0.0, 0.0), &grs80());
        let out = epoch_op
            .forward_4d(Coord::new(p.x, p.y, p.z, 2010.0))
            .expect("forward_4d");
        let v = out.v();
        assert!((v[3] - 2010.0).abs() < 1e-9, "time passthrough");
        assert!((v[0] - (p.x + 3.0)).abs() < 1e-9, "dX(u)");
        assert!((v[1] - (p.y + 1.0)).abs() < 1e-9, "dY(e)");
        assert!((v[2] - (p.z + 2.0)).abs() < 1e-9, "dZ(n)");
    }

    #[test]
    fn epoch_mode_requires_finite_time() {
        let epoch_op = Deformation {
            grid: make_velocity_grid(&[1.0, 2.0, 3.0]),
            ell: grs80(),
            sample_e: 0,
            sample_n: 1,
            sample_u: 2,
            fixed_dt: None,
            t_epoch: Some(2000.0),
        };
        // Missing/non-finite time in epoch mode is an error.
        let err = epoch_op
            .forward_4d(Coord::new(0.0, 0.0, 0.0, f64::INFINITY))
            .err();
        assert_eq!(err, Some(ProjError::MissingTime));
    }

    #[test]
    fn three_d_entry_requires_fixed_span() {
        // Epoch-only op has no fixed span, so the 3D entry point errors.
        let epoch_op = Deformation {
            grid: make_velocity_grid(&[1.0, 2.0, 3.0]),
            ell: grs80(),
            sample_e: 0,
            sample_n: 1,
            sample_u: 2,
            fixed_dt: None,
            t_epoch: Some(2000.0),
        };
        assert_eq!(
            epoch_op.forward_3d(Lpz::new(0.0, 0.0, 0.0)).err(),
            Some(ProjError::MissingTime)
        );
    }

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

    // --- constructor parameter validation (no grid resolution required) ---

    use crate::TransParamLookup;

    #[derive(Default)]
    struct Params {
        grids: Option<String>,
        dt: Option<f64>,
        t_epoch: Option<f64>,
        t_obs: Option<f64>,
    }

    impl TransParamLookup for Params {
        fn get_dms(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_f64(&self, key: &str) -> Option<f64> {
            match key {
                "dt" => self.dt,
                "t_epoch" => self.t_epoch,
                "t_obs" => self.t_obs,
                _ => None,
            }
        }
        fn get_int(&self, _key: &str) -> Option<i64> {
            None
        }
        fn get_str(&self, key: &str) -> Option<&str> {
            match key {
                "grids" => self.grids.as_deref(),
                _ => None,
            }
        }
        fn get_bool(&self, _key: &str) -> bool {
            false
        }
        fn exists(&self, key: &str) -> bool {
            match key {
                "grids" => self.grids.is_some(),
                "dt" => self.dt.is_some(),
                "t_epoch" => self.t_epoch.is_some(),
                "t_obs" => self.t_obs.is_some(),
                _ => false,
            }
        }
    }

    fn build_new(params: Params) -> ProjResult<TransBuild> {
        let ell = grs80();
        new(&TransParams {
            ellipsoid: &ell,
            params: &params,
            registry: None,
        })
    }

    #[test]
    fn missing_grids_errors() {
        let err = build_new(Params {
            dt: Some(1.0),
            ..Default::default()
        })
        .err();
        assert_eq!(err, Some(ProjError::MissingArg));
    }

    #[test]
    fn missing_time_span_errors() {
        // grids present but no dt/t_epoch/t_obs → MissingArg (before resolving).
        let err = build_new(Params {
            grids: Some("whatever.tif".to_string()),
            ..Default::default()
        })
        .err();
        assert_eq!(err, Some(ProjError::MissingArg));
    }

    #[test]
    fn dt_and_epoch_mutually_exclusive() {
        let err = build_new(Params {
            grids: Some("whatever.tif".to_string()),
            dt: Some(1.0),
            t_epoch: Some(2000.0),
            ..Default::default()
        })
        .err();
        assert_eq!(err, Some(ProjError::MutuallyExclusiveArgs));
    }

    #[test]
    fn t_obs_without_epoch_errors() {
        let err = build_new(Params {
            grids: Some("whatever.tif".to_string()),
            t_obs: Some(2020.0),
            ..Default::default()
        })
        .err();
        assert_eq!(err, Some(ProjError::MissingArg));
    }
}