oxiproj-transformations 0.1.2

//! `cart` conversion: geographic (lon, lat, h) <-> geocentric Cartesian (X, Y, Z).
//!
//! Ported from PROJ src/conversions/cart.cpp. PROJ declares this operation with
//! `left = RADIANS, right = CARTESIAN`, so the forward direction maps geodetic
//! coordinates (in radians, with ellipsoidal height) to geocentric Cartesian
//! coordinates (in metres), and the inverse maps them back. The underlying math
//! lives in `crate::internal` (`geodetic_to_cartesian` / `cartesian_to_geodetic`),
//! shared with the `geocent` and `topocentric` conversions.

use crate::{TransBuild, TransParams};
use oxiproj_core::{Ellipsoid, IoUnits, Lp, Lpz, Operation, ProjResult, Xy, Xyz};

/// The `cart` conversion operator. Ported from src/conversions/cart.cpp.
///
/// Stores a copy of the ellipsoid (`Ellipsoid` is `Copy`) so the operation is
/// self-contained.
#[derive(Debug)]
struct Cart {
    ell: Ellipsoid,
}

impl Operation for Cart {
    /// Geographic 2D (lon, lat) -> geocentric (X, Y), assuming zero height.
    /// Ported from src/conversions/cart.cpp (`forward`, height-0 case).
    fn forward_2d(&self, lp: Lp) -> ProjResult<Xy> {
        let xyz = crate::internal::geodetic_to_cartesian(Lpz::new(lp.lam, lp.phi, 0.0), &self.ell);
        Ok(Xy::new(xyz.x, xyz.y))
    }

    /// Geocentric (X, Y) -> geographic 2D (lon, lat), assuming zero Z.
    /// Ported from src/conversions/cart.cpp (`reverse`, Z-0 case).
    fn inverse_2d(&self, xy: Xy) -> ProjResult<Lp> {
        let lpz = crate::internal::cartesian_to_geodetic(Xyz::new(xy.x, xy.y, 0.0), &self.ell);
        Ok(Lp::new(lpz.lam, lpz.phi))
    }

    /// Geographic (lon, lat, h) -> geocentric Cartesian (X, Y, Z).
    /// Ported from src/conversions/cart.cpp (`cartesian` / `forward`).
    fn forward_3d(&self, lpz: Lpz) -> ProjResult<Xyz> {
        Ok(crate::internal::geodetic_to_cartesian(lpz, &self.ell))
    }

    /// Geocentric Cartesian (X, Y, Z) -> geographic (lon, lat, h).
    /// Ported from src/conversions/cart.cpp (`geodetic` / `reverse`).
    fn inverse_3d(&self, xyz: Xyz) -> ProjResult<Lpz> {
        Ok(crate::internal::cartesian_to_geodetic(xyz, &self.ell))
    }

    /// The `cart` conversion is always invertible.
    fn has_inverse(&self) -> bool {
        true
    }
}

/// Construct the `cart` conversion. Ported from src/conversions/cart.cpp
/// (`PJ_CONVERSION(cart)`): `left = RADIANS`, `right = CARTESIAN`.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    Ok(TransBuild::new(
        Box::new(Cart { ell: *p.ellipsoid }),
        IoUnits::Radians,
        IoUnits::Cartesian,
    ))
}

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

    struct NoParams;
    impl TransParamLookup for NoParams {
        fn get_dms(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_f64(&self, _key: &str) -> Option<f64> {
            None
        }
        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 {
            false
        }
    }

    #[test]
    fn cart_forward_inverse_wgs84() {
        let ell = Ellipsoid::named("WGS84").unwrap();
        let np = NoParams;
        let b = new(&TransParams {
            ellipsoid: &ell,
            params: &np,
            registry: None,
        })
        .unwrap();
        let op = b.operation;

        let lam = 12.0 * DEG_TO_RAD;
        let phi = 55.0 * DEG_TO_RAD;
        let h = 100.0;

        let xyz = op.forward_3d(Lpz::new(lam, phi, h)).unwrap();
        assert!((xyz.x - 3586525.761017917).abs() < 1e-4, "X = {}", xyz.x);
        assert!((xyz.y - 762339.584102928).abs() < 1e-4, "Y = {}", xyz.y);
        assert!((xyz.z - 5201465.438406702).abs() < 1e-4, "Z = {}", xyz.z);

        let back = op.inverse_3d(xyz).unwrap();
        assert!((back.lam - lam).abs() < 1e-9, "lam = {}", back.lam);
        assert!((back.phi - phi).abs() < 1e-9, "phi = {}", back.phi);
        assert!((back.z - h).abs() < 1e-4, "z = {}", back.z);
    }

    #[test]
    fn cart_forward_4d_passes_t() {
        let ell = Ellipsoid::named("WGS84").unwrap();
        let np = NoParams;
        let b = new(&TransParams {
            ellipsoid: &ell,
            params: &np,
            registry: None,
        })
        .unwrap();
        let op = b.operation;

        let out = op
            .forward_4d(Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 100.0, 7.0))
            .unwrap();
        assert_eq!(out.v()[3], 7.0);
    }
}