oxiproj-transformations 0.1.2

//! `geocent` conversion. Ported from PROJ `src/conversions/geocent.cpp`.
//!
//! In PROJ, `geocent` is declared with `left = RADIANS`, `right = CARTESIAN`
//! and `P->is_geocent = 1`. Its 2D `forward`/`inverse` callbacks are a
//! pass-through stub (`xy.x = lp.lam`, `xy.y = lp.phi` and back); the genuine
//! geodetic <-> geocentric Cartesian conversion is performed *around* that stub
//! by `fwd_finalize`/`inv_prepare`, which — because `P->is_geocent` is set —
//! run the coordinate through an internal `cart` object (`src/fwd.cpp`
//! `if (P->is_geocent) coo = proj_trans(P->cart, PJ_FWD, coo)` and the mirror in
//! `src/inv.cpp`).
//!
//! OxiProj drives conversions with prepare/finalize bypassed
//! (`registry.rs`: transformations set `bypass_prepare_finalize = true`), so the
//! `is_geocent` finalize/prepare hooks never fire for a standalone
//! `+proj=geocent`. To match PROJ end-to-end we therefore fold the `is_geocent`
//! Cartesian round-trip directly into the operation: the forward direction maps
//! geodetic (lon, lat, h) in radians/metres to geocentric Cartesian (X, Y, Z) in
//! metres, and the inverse maps them back — identical math to the `cart`
//! conversion, sharing the helpers in the crate-internal `internal` module.

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

/// The `geocent` conversion operator.
///
/// PROJ's `geocent` is a 2D stub whose real geodetic <-> geocentric Cartesian
/// conversion is supplied by the `is_geocent` branch of `fwd_finalize` /
/// `inv_prepare` via an internal `cart` object. OxiProj bypasses
/// prepare/finalize for conversions, so this operation performs that Cartesian
/// conversion itself (equivalent to the stub plus the `is_geocent` cart
/// round-trip), storing a copy of the ellipsoid (`Ellipsoid` is `Copy`).
#[derive(Debug)]
struct Geocent {
    ell: Ellipsoid,
}

impl Operation for Geocent {
    /// Geographic 2D (lon, lat) -> geocentric (X, Y), assuming zero height.
    ///
    /// Equivalent to PROJ's `geocent` 2D stub followed by the `is_geocent`
    /// `cart` forward with `h = 0` in `fwd_finalize`.
    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.
    ///
    /// Equivalent to the `is_geocent` `cart` inverse with `Z = 0` in
    /// `inv_prepare` followed by PROJ's `geocent` 2D stub inverse.
    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 the `is_geocent` `cart` forward in `src/fwd.cpp`
    /// (`fwd_finalize`), sharing [`crate::internal::geodetic_to_cartesian`] with
    /// the `cart` conversion.
    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 the `is_geocent` `cart` inverse in `src/inv.cpp`
    /// (`inv_prepare`), sharing [`crate::internal::cartesian_to_geodetic`] with
    /// the `cart` conversion.
    fn inverse_3d(&self, xyz: Xyz) -> ProjResult<Lpz> {
        Ok(crate::internal::cartesian_to_geodetic(xyz, &self.ell))
    }

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

/// Construct the `geocent` conversion.
///
/// Ported from `src/conversions/geocent.cpp` (`PJ_CONVERSION(geocent)`). The
/// geodetic input is in radians (with ellipsoidal height in metres) and the
/// output is geocentric Cartesian (metres). The real conversion is carried out
/// by the operation itself, matching PROJ's `is_geocent` `cart` round-trip that
/// `fwd_finalize`/`inv_prepare` apply around the historical 2D stub.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    Ok(TransBuild::new(
        Box::new(Geocent { 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
        }
    }

    fn build_wgs84() -> TransBuild {
        let ell = Ellipsoid::named("WGS84").unwrap();
        let np = NoParams;
        let tp = TransParams {
            ellipsoid: &ell,
            params: &np,
            registry: None,
        };
        new(&tp).unwrap()
    }

    #[test]
    fn geocent_build_units() {
        let b = build_wgs84();
        assert_eq!(b.left, IoUnits::Radians);
        assert_eq!(b.right, IoUnits::Cartesian);
    }

    #[test]
    fn geocent_forward_matches_proj_ecef() {
        // Ground truth from Homebrew PROJ 9.x:
        //   echo "12 55 100" | cct +proj=geocent +ellps=WGS84
        //   -> 3586525.7610  762339.5841  5201465.4384
        let op = build_wgs84().operation;
        let out = op
            .forward_3d(Lpz::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 100.0))
            .unwrap();
        assert!((out.x - 3586525.761017917).abs() < 1e-4, "X = {}", out.x);
        assert!((out.y - 762339.584102928).abs() < 1e-4, "Y = {}", out.y);
        assert!((out.z - 5201465.438406702).abs() < 1e-4, "Z = {}", out.z);
    }

    #[test]
    fn geocent_inverse_matches_proj_geodetic() {
        // Ground truth from Homebrew PROJ 9.x:
        //   echo "3586525.7610 762339.5841 5201465.4384" | cct -I +proj=geocent +ellps=WGS84
        //   -> 12  55  100
        let op = build_wgs84().operation;
        let lpz = op
            .inverse_3d(Xyz::new(
                3586525.761017917,
                762339.584102928,
                5201465.438406702,
            ))
            .unwrap();
        assert!(
            (lpz.lam - 12.0 * DEG_TO_RAD).abs() < 1e-9,
            "lam = {}",
            lpz.lam
        );
        assert!(
            (lpz.phi - 55.0 * DEG_TO_RAD).abs() < 1e-9,
            "phi = {}",
            lpz.phi
        );
        assert!((lpz.z - 100.0).abs() < 1e-4, "z = {}", lpz.z);
    }

    #[test]
    fn geocent_round_trip_preserves_time() {
        let op = build_wgs84().operation;
        let input = Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 100.0, 2020.5);
        let ecef = op.forward_4d(input).unwrap();
        // Real ECEF, not a pass-through of the input radians.
        assert!(
            ecef.v()[0] > 3.0e6,
            "X should be ECEF metres, got {}",
            ecef.v()[0]
        );
        assert_eq!(ecef.v()[3], 2020.5, "time must pass through");
        let back = op.inverse_4d(ecef).unwrap();
        assert!((back.v()[0] - input.v()[0]).abs() < 1e-9);
        assert!((back.v()[1] - input.v()[1]).abs() < 1e-9);
        assert!((back.v()[2] - input.v()[2]).abs() < 1e-4);
        assert_eq!(back.v()[3], 2020.5);
    }
}