oxiproj-transformations 0.1.1

#![forbid(unsafe_code)]
//! Velocity-grid deformation (`deformation`) — port of PROJ `src/transformations/deformation.cpp`.
//!
//! Applies a time-dependent coordinate shift derived from a velocity grid:
//!   shift = velocity × dt, where dt = coordinate_time − t_epoch
//!
//! Parameters:
//!   `+grids=<name>` — name of the velocity grid in the registry (2 or 3 bands)
//!   `+t_epoch=<f64>` — reference epoch (decimal year)
//!
//! Grid band convention (same units as the input coordinate space):
//!   band 0 = eastward velocity (degrees/year for geographic, meters/year for geocentric)
//!   band 1 = northward velocity (degrees/year for geographic, meters/year for geocentric)
//!   band 2 = upward velocity (meters/year), optional
//!
//! For this implementation we assume a geographic input (radians in, radians out) with
//! horizontal velocities in degrees/year and vertical velocity in meters/year.

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

#[derive(Debug)]
struct Deformation {
    grid: GridSet,
    t_epoch: f64,
}

impl Operation for Deformation {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let t = v[3];
        let dt = t - self.t_epoch;
        let lon_deg = v[0].to_degrees();
        let lat_deg = v[1].to_degrees();
        let velocities = sample_grid(&self.grid, lat_deg, lon_deg).ok_or(ProjError::OutsideGrid)?;
        let dvlon_deg = velocities[0] * dt;
        let dvlat_deg = if velocities.len() >= 2 {
            velocities[1] * dt
        } else {
            0.0
        };
        let dvz = if velocities.len() >= 3 {
            velocities[2] * dt
        } else {
            0.0
        };
        Ok(Coord::new(
            v[0] + dvlon_deg.to_radians(),
            v[1] + dvlat_deg.to_radians(),
            v[2] + dvz,
            v[3],
        ))
    }

    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        // Inverse: subtract the velocity displacement.
        let v = c.v();
        let t = v[3];
        let dt = t - self.t_epoch;
        let lon_deg = v[0].to_degrees();
        let lat_deg = v[1].to_degrees();
        let velocities = sample_grid(&self.grid, lat_deg, lon_deg).ok_or(ProjError::OutsideGrid)?;
        let dvlon_deg = velocities[0] * dt;
        let dvlat_deg = if velocities.len() >= 2 {
            velocities[1] * dt
        } else {
            0.0
        };
        let dvz = if velocities.len() >= 3 {
            velocities[2] * dt
        } else {
            0.0
        };
        Ok(Coord::new(
            v[0] - dvlon_deg.to_radians(),
            v[1] - dvlat_deg.to_radians(),
            v[2] - dvz,
            v[3],
        ))
    }

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

/// Construct a `deformation` transform from parsed parameters.
///
/// Requires:
///   `+grids=<name>` — velocity grid (NTv2 or GTX format, tried in order)
///   `+t_epoch=<f64>` — reference epoch in decimal years
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let grid_name = p.params.get_str("grids").ok_or(ProjError::MissingArg)?;
    let t_epoch = p.params.get_f64("t_epoch").ok_or(ProjError::MissingArg)?;
    let registry = p.registry.ok_or(ProjError::FileNotFound)?;
    let data = registry
        .get_grid(grid_name)
        .ok_or(ProjError::FileNotFound)?;

    // Try NTv2 first (supports multi-band grids).
    if let Ok(grids) = read_ntv2(data, grid_name) {
        if let Some(grid) = grids.into_iter().next() {
            return Ok(TransBuild::new(
                Box::new(Deformation { grid, t_epoch }),
                IoUnits::Radians,
                IoUnits::Radians,
            ));
        }
    }

    // Fall back to GTX (single band).
    let grid = read_gtx(data, grid_name)?;
    Ok(TransBuild::new(
        Box::new(Deformation { grid, t_epoch }),
        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 given band values (all uniform per band).
    fn make_velocity_grid(band_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 = band_values
            .iter()
            .map(|&v| GridBand {
                extent: extent.clone(),
                values: vec![v; 4],
            })
            .collect();
        GridSet {
            bands,
            source: GridSource::Memory,
        }
    }

    #[test]
    fn test_deformation_forward_2band() {
        // Velocity: vlon=0.1 deg/yr, vlat=0.2 deg/yr
        let grid = make_velocity_grid(&[0.1f32, 0.2]);
        let op = Deformation {
            grid,
            t_epoch: 2000.0,
        };
        // At t=2010 → dt=10 → shift=(1 deg, 2 deg)
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 2010.0);
        let out = op.forward_4d(input).unwrap();
        let v = out.v();
        let expected_lon = (0.5 + 1.0) * DEG_TO_RAD;
        let expected_lat = (0.5 + 2.0) * DEG_TO_RAD;
        assert!((v[0] - expected_lon).abs() < 1e-9, "lon");
        assert!((v[1] - expected_lat).abs() < 1e-9, "lat");
    }

    #[test]
    fn test_deformation_forward_3band() {
        // Band 2: vertical velocity 0.5 m/yr → dt=10 → dz=5 m
        let grid = make_velocity_grid(&[0.0f32, 0.0, 0.5]);
        let op = Deformation {
            grid,
            t_epoch: 2000.0,
        };
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 100.0, 2010.0);
        let out = op.forward_4d(input).unwrap();
        assert!((out.v()[2] - 105.0).abs() < 1e-6, "z should be 105");
    }

    #[test]
    fn test_deformation_zero_dt() {
        // dt=0 → no shift
        let grid = make_velocity_grid(&[1.0f32, 1.0]);
        let op = Deformation {
            grid,
            t_epoch: 2000.0,
        };
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 2000.0);
        let out = op.forward_4d(input).unwrap();
        let v = out.v();
        assert!((v[0] - input.v()[0]).abs() < 1e-12, "no lon shift");
        assert!((v[1] - input.v()[1]).abs() < 1e-12, "no lat shift");
    }

    #[test]
    fn test_deformation_inverse() {
        // Use very small velocities so that dt=10 produces a tiny shift (0.001/0.002 deg)
        // that keeps the shifted point well inside the [0,1]×[0,1] grid extent.
        let grid = make_velocity_grid(&[0.0001f32, 0.0002]);
        let op = Deformation {
            grid,
            t_epoch: 2000.0,
        };
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 2010.0);
        let fwd = op.forward_4d(input).unwrap();
        // For inverse, we look up velocity at the SHIFTED position, not original.
        // Since our velocity is uniform (same everywhere), the round-trip is exact.
        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");
    }

    #[test]
    fn test_deformation_outside_returns_err() {
        let grid = make_velocity_grid(&[0.1f32, 0.2]);
        let op = Deformation {
            grid,
            t_epoch: 2000.0,
        };
        let far = Coord::new(45.0 * DEG_TO_RAD, 45.0 * DEG_TO_RAD, 0.0, 2010.0);
        assert_eq!(op.forward_4d(far).err(), Some(ProjError::OutsideGrid));
    }
}