oxiproj-transformations 0.1.2

//! `vertoffset` transformation ("Vertical Offset and Slope").
//!
//! Ported from PROJ 9.8.0 `src/transformations/vertoffset.cpp`.
//!
//! Implements EPSG Dataset coordinate operation method code 1046 "Vertical
//! Offset and Slope". The transformation modifies only the vertical (`z`)
//! component; longitude and latitude pass through unchanged. The applied offset
//! is a base vertical offset (`dh`) plus latitude- and longitude-dependent slope
//! terms scaled by the meridian (`rho0`) and prime-vertical (`nu0`) radii of
//! curvature evaluated at the origin latitude `phi0`.
//!
//! Note on `lam0`: PROJ's transformation engine subtracts `lam0` before the
//! forward step and adds it back afterwards (desirable for map projections but
//! not for this z-only method — see `forward_3d` / `reverse_3d` in the original
//! source). OxiProj's engine does not perform that adjustment here, so to
//! reproduce `cct` faithfully the longitude slope term uses `(lam - lam0)`
//! directly.

use crate::{TransBuild, TransParams};
use oxiproj_core::{Coord, IoUnits, Lpz, Operation, ProjResult, Xyz, DEG_TO_RAD};

/// Arcseconds to radians. Ported from PROJ `src/transformations/vertoffset.cpp`
/// (`#define ARCSEC_TO_RAD (DEG_TO_RAD / 3600.0)`).
const ARCSEC_TO_RAD: f64 = DEG_TO_RAD / 3600.0;

/// The `vertoffset` transformation. Ported from PROJ
/// `src/transformations/vertoffset.cpp` (`struct pj_opaque_vertoffset`).
#[derive(Debug)]
struct VertOffset {
    /// Origin longitude (`lon_0`), radians.
    lam0: f64,
    /// Origin latitude (`lat_0`), radians.
    phi0: f64,
    /// Latitude slope (`slope_lat`), radians per metre of along-meridian arc.
    slope_lat: f64,
    /// Longitude slope (`slope_lon`), radians per metre of along-parallel arc.
    slope_lon: f64,
    /// Base vertical offset (`dh`), metres.
    zoff: f64,
    /// Meridian radius of curvature at `phi0`.
    rho0: f64,
    /// Prime-vertical radius of curvature at `phi0`.
    nu0: f64,
}

impl VertOffset {
    /// Compute the vertical offset to add for a given geodetic position.
    ///
    /// Ported from PROJ `src/transformations/vertoffset.cpp`
    /// (`get_forward_offset`), using `(lam - lam0)` for the longitude term (see
    /// the module-level note on `lam0`).
    fn forward_offset(&self, phi: f64, lam: f64) -> f64 {
        let dlam = lam - self.lam0;
        self.zoff
            + self.slope_lat * self.rho0 * (phi - self.phi0)
            + self.slope_lon * self.nu0 * dlam * phi.cos()
    }
}

impl Operation for VertOffset {
    /// Ported from PROJ `src/transformations/vertoffset.cpp` (`forward_3d`).
    /// Longitude (`x = lam`) and latitude (`y = phi`) pass through unchanged;
    /// only `z` is modified.
    fn forward_3d(&self, lpz: Lpz) -> ProjResult<Xyz> {
        let off = self.forward_offset(lpz.phi, lpz.lam);
        Ok(Xyz::new(lpz.lam, lpz.phi, lpz.z + off))
    }

    /// Ported from PROJ `src/transformations/vertoffset.cpp` (`reverse_3d`).
    /// Here the incoming `x` is longitude and `y` is latitude; only `z` is
    /// modified (the offset is subtracted).
    fn inverse_3d(&self, xyz: Xyz) -> ProjResult<Lpz> {
        let off = self.forward_offset(xyz.y, xyz.x);
        Ok(Lpz::new(xyz.x, xyz.y, xyz.z - off))
    }

    /// Route the 4D forward through the 3D operation, preserving the time slot.
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let o = self.forward_3d(Lpz::new(v[0], v[1], v[2]))?;
        Ok(Coord::new(o.x, o.y, o.z, v[3]))
    }

    /// Route the 4D inverse through the 3D operation, preserving the time slot.
    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let o = self.inverse_3d(Xyz::new(v[0], v[1], v[2]))?;
        Ok(Coord::new(o.lam, o.phi, o.z, v[3]))
    }

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

/// Construct the `vertoffset` transformation. Ported from PROJ
/// `src/transformations/vertoffset.cpp` (`PJ_TRANSFORMATION(vertoffset, 1)`).
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let pa = p.params;
    let phi0 = pa.get_dms("lat_0").unwrap_or(0.0);
    let lam0 = pa.get_dms("lon_0").unwrap_or(0.0);
    let zoff = pa.get_f64("dh").unwrap_or(0.0);
    let slope_lat = pa.get_f64("slope_lat").unwrap_or(0.0) * ARCSEC_TO_RAD;
    let slope_lon = pa.get_f64("slope_lon").unwrap_or(0.0) * ARCSEC_TO_RAD;

    let a = p.ellipsoid.a;
    let es = p.ellipsoid.es;
    let sinlat0 = phi0.sin();
    let one_minus_es_sinlat0_square = 1.0 - es * (sinlat0 * sinlat0);
    let rho0 = a * (1.0 - es) / (one_minus_es_sinlat0_square * one_minus_es_sinlat0_square.sqrt());
    let nu0 = a / one_minus_es_sinlat0_square.sqrt();

    let op = VertOffset {
        lam0,
        phi0,
        slope_lat,
        slope_lon,
        zoff,
        rho0,
        nu0,
    };
    Ok(TransBuild::new(
        Box::new(op),
        IoUnits::Radians,
        IoUnits::Radians,
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::TransParamLookup;
    use oxiproj_core::{Ellipsoid, DEG_TO_RAD};
    use std::collections::HashMap;

    /// Parameter lookup backed by two maps: `dms` holds already-radian angle
    /// values (for `lat_0` / `lon_0`); `f64` holds plain numbers (`dh`,
    /// `slope_lat`, `slope_lon`).
    struct ProjParamLookup {
        dms: HashMap<&'static str, f64>,
        f64s: HashMap<&'static str, f64>,
    }

    impl TransParamLookup for ProjParamLookup {
        fn get_dms(&self, key: &str) -> Option<f64> {
            self.dms.get(key).copied()
        }
        fn get_f64(&self, key: &str) -> Option<f64> {
            self.f64s.get(key).copied()
        }
        fn get_int(&self, _key: &str) -> Option<i64> {
            None
        }
        fn get_str(&self, _key: &str) -> Option<&str> {
            None
        }
        fn get_bool(&self, _key: &str) -> bool {
            false
        }
        fn exists(&self, key: &str) -> bool {
            self.dms.contains_key(key) || self.f64s.contains_key(key)
        }
    }

    fn build(dms: &[(&'static str, f64)], f64s: &[(&'static str, f64)]) -> TransBuild {
        let ell = Ellipsoid::named("WGS84").unwrap();
        let params = ProjParamLookup {
            dms: dms.iter().copied().collect(),
            f64s: f64s.iter().copied().collect(),
        };
        let pp = TransParams {
            ellipsoid: &ell,
            params: &params,
            registry: None,
        };
        new(&pp).unwrap()
    }

    #[test]
    fn cct_origin_point() {
        // lat_0=55, lon_0=12, dh=1.5, slope_lat=0.1, slope_lon=0.2.
        let b = build(
            &[("lat_0", 55.0 * DEG_TO_RAD), ("lon_0", 12.0 * DEG_TO_RAD)],
            &[("dh", 1.5), ("slope_lat", 0.1), ("slope_lon", 0.2)],
        );
        // Input at the origin: phi == phi0 and lam == lam0, so both slope terms
        // vanish and only `dh` is applied: z = 100 + 1.5 = 101.5.
        let out = b
            .operation
            .forward_3d(Lpz::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 100.0))
            .unwrap();
        assert!((out.z - 101.5).abs() < 1e-6, "z = {}", out.z);
        // Longitude (x) and latitude (y) pass through unchanged.
        assert!((out.x - 12.0 * DEG_TO_RAD).abs() < 1e-12, "x = {}", out.x);
        assert!((out.y - 55.0 * DEG_TO_RAD).abs() < 1e-12, "y = {}", out.y);
    }

    #[test]
    fn slope_round_trip() {
        // Same params, but latitude 56 deg (lon stays 12 -> dlam == 0), so only
        // the latitude slope term is active.
        let b = build(
            &[("lat_0", 55.0 * DEG_TO_RAD), ("lon_0", 12.0 * DEG_TO_RAD)],
            &[("dh", 1.5), ("slope_lat", 0.1), ("slope_lon", 0.2)],
        );
        let input = Lpz::new(12.0 * DEG_TO_RAD, 56.0 * DEG_TO_RAD, 100.0);
        let fwd = b.operation.forward_3d(input).unwrap();
        // The applied offset must differ from the base 1.5: the latitude slope
        // term is non-trivial.
        let offset = fwd.z - 100.0;
        assert!(
            (offset - 1.5).abs() > 1e-9,
            "slope_lat term not exercised: offset = {offset}"
        );
        // forward then inverse recovers the input height.
        let back = b
            .operation
            .inverse_3d(Xyz::new(fwd.x, fwd.y, fwd.z))
            .unwrap();
        assert!((back.z - input.z).abs() < 1e-9, "z = {}", back.z);
    }
}