oxiproj_engine/

pj.rs

//! The `Pj` operation wrapper: prepare/finalize around an inner operation.
//!
//! Ported from PROJ 9.8.0 `src/fwd.cpp` (`fwd_prepare`/`fwd_finalize`) and
//! `src/inv.cpp` (`inv_prepare`/`inv_finalize`).

use std::sync::Arc;

use oxiproj_core::{
    Coord, Ellipsoid, IoUnits, Lp, Operation, ProjError, ProjResult, ProjectGenericBox, Xy,
};

/// Latitude tolerance for the prepare-stage domain check (radians).
const PJ_EPS_LAT: f64 = 1e-12;

/// A prepared coordinate operation with input/output unit handling.
///
/// Wraps an inner [`Operation`] and applies the unit conversions, prime-meridian
/// and central-meridian shifts, and axis bookkeeping that PROJ performs around
/// the raw projection math.
///
/// # Field visibility and invariants
///
/// A `Pj` is a *prepared* operation: its fields are mutually consistent by
/// construction (only [`crate::create`] and the crate-internal
/// `registry::build_single_op` produce a valid one). All fields are therefore
/// `pub(crate)` and reachable
/// from outside the crate only through the read accessors on [`Pj`]. Flipping,
/// for example, [`inverted`](Pj::inverted) or `bypass_prepare_finalize` after
/// construction would silently invalidate the prepare/finalize contract, so
/// those knobs are deliberately not externally writable.
///
/// The operation's *boundary* — its reference ellipsoid and forward-sense
/// input/output units — is exposed read-only through the
/// [`ellipsoid`](Pj::ellipsoid), [`left`](Pj::left) and [`right`](Pj::right)
/// accessors, which downstream crates (the CLI, parity harnesses) depend on.
/// Their invariant: `left`/`right` are the *forward-sense* input/output units of
/// the whole operation, and `ellipsoid` is the reference ellipsoid the boundary
/// math is expressed against. They are read-only; mutating them would not
/// re-derive any of the dependent fields.
#[derive(Debug)]
pub struct Pj {
    /// The inner operation (projection or conversion) driven via 4D.
    pub(crate) operation: Box<dyn Operation>,
    /// The reference ellipsoid.
    ///
    /// Stable read-only boundary field (see the type-level invariants);
    /// read from outside the crate via [`Pj::ellipsoid`].
    pub(crate) ellipsoid: Ellipsoid,
    /// Effective squared eccentricity used by the map-scale factor / Tissot
    /// metric ([`Pj::factors`] / [`Pj::factors_exact`]).
    ///
    /// PROJ forces `P->es = 0` in the setup of every sphere-only projection
    /// (`src/projections/*.cpp`), so its factors are computed on the sphere even
    /// when the user supplied an ellipsoid. OxiProj mirrors that: for a
    /// projection whose forward ignores `es` (detected empirically at build
    /// time), this is `0.0`; for a genuinely ellipsoidal projection it equals
    /// `ellipsoid.es`. The corresponding `one_es` is `1.0 - factors_es`.
    pub(crate) factors_es: f64,
    /// Central meridian (radians).
    pub(crate) lam0: f64,
    /// Latitude of origin (radians).
    pub(crate) phi0: f64,
    /// False easting (meters, pre-scale).
    pub(crate) x0: f64,
    /// False northing (meters, pre-scale).
    pub(crate) y0: f64,
    /// False up (meters, pre-scale).
    pub(crate) z0: f64,
    /// Scale factor at the central meridian.
    pub(crate) k0: f64,
    /// Horizontal input-unit-to-meter factor.
    pub(crate) to_meter: f64,
    /// Horizontal meter-to-output-unit factor.
    pub(crate) fr_meter: f64,
    /// Vertical input-unit-to-meter factor.
    pub(crate) vto_meter: f64,
    /// Vertical meter-to-output-unit factor.
    pub(crate) vfr_meter: f64,
    /// Prime-meridian offset from Greenwich (radians).
    pub(crate) from_greenwich: f64,
    /// Allow longitudes outside the normal range (`+over`).
    pub(crate) over: bool,
    /// Geocentric latitude flag (`+geoc`).
    pub(crate) geoc: bool,
    /// Longitude wrap center for `+lon_wrap` (radians), applied to angular
    /// (`IoUnits::Radians`) output on the operation's own forward pass —
    /// PROJ's `fwd_finalize` `is_long_wrap_set` branch (`src/fwd.cpp`):
    /// output longitude is normalized into `[center - pi, center + pi]` via
    /// `center + adjlon(lam - center)`, taking precedence over the plain
    /// `+over`-gated `adjlon` wrap. `None` when `+lon_wrap` was not supplied
    /// (the ordinary `+over` wrap applies instead). PROJ's `inv_finalize`
    /// never applies this wrap (only `fwd_finalize` does), so it is
    /// consulted only on the operation's true forward direction — see
    /// [`Pj::forward`]/[`Pj::inverse`]'s dispatch through `bypass_dispatch`.
    pub(crate) lon_wrap_center: Option<f64>,
    /// Whether this is a pass-through geographic (lat/long) operation.
    pub(crate) is_latlong: bool,
    /// Input units (forward direction).
    ///
    /// Stable read-only boundary field (see the type-level invariants);
    /// read from outside the crate via [`Pj::left`].
    pub(crate) left: IoUnits,
    /// Output units (forward direction).
    ///
    /// Stable read-only boundary field (see the type-level invariants);
    /// read from outside the crate via [`Pj::right`].
    pub(crate) right: IoUnits,
    /// Whether the operation runs inverted in its containing pipeline.
    pub(crate) inverted: bool,
    /// Skip prepare/finalize entirely (used by the top-level pipeline driver).
    pub(crate) bypass_prepare_finalize: bool,
    /// Skip this step when the pipeline runs in the forward direction.
    pub(crate) omit_fwd: bool,
    /// Skip this step when the pipeline runs in the inverse direction.
    pub(crate) omit_inv: bool,
    /// Optional AD-capable projection for exact Jacobian computation (T7.1).
    ///
    /// When `Some`, `factors_exact` uses this to drive forward-mode AD via
    /// `Dual1<2>`. Set to `None` for projections that do not implement
    /// `ProjectGeneric` or whose math is not differentiable.
    pub(crate) ad_proj: Option<Arc<dyn ProjectGenericBox>>,
    /// The operation name (e.g. `"helmert"`, `"axisswap"`, `"noop"`).
    ///
    /// Populated from the `+proj=` key at construction time. Used by the
    /// pipeline optimizer (`pipeline_opt`) to detect and eliminate
    /// redundant steps. Empty for the outer pipeline wrapper itself.
    pub(crate) op_name: String,
    /// Canonical signed axis permutation, populated only for `axisswap` steps.
    ///
    /// Encodes the operation's *non-inverted* forward mapping in PROJ
    /// `+order=` form: `axisswap_order[i] = sign * (source_axis + 1)`, so entry
    /// `i` (a value in `±1..=±4`) means output slot `i` takes input slot
    /// `|value| - 1` multiplied by `signum(value)`. It is derived by probing the
    /// built operation (see `registry::probe_axisswap_order`) rather than by
    /// re-parsing parameters, so it always reflects the operation's true
    /// behaviour. The pipeline optimizer uses this — together with the
    /// per-step [`inverted`](Pj::inverted) flag — to cancel a pair of adjacent
    /// `axisswap` steps only when they compose to the exact identity. `None`
    /// for every non-`axisswap` operation.
    pub(crate) axisswap_order: Option<[i8; 4]>,
}

impl Pj {
    /// The inner operation (projection or conversion) as a trait reference.
    pub fn operation(&self) -> &dyn Operation {
        &*self.operation
    }
    /// The operation name (e.g. `"merc"`, `"axisswap"`); empty for a pipeline.
    pub fn op_name(&self) -> &str {
        &self.op_name
    }
    /// Whether this operation runs inverted in its containing pipeline.
    pub fn inverted(&self) -> bool {
        self.inverted
    }
    /// Whether the wrapped inner operation provides an inverse.
    ///
    /// Unlike the [`Operation::has_inverse`] default (which the `Pj` wrapper
    /// inherits as always-`true`), this reports the *inner* operation's own
    /// answer — e.g. `false` for a forward-only projection like `urm5` or a
    /// singular `affine`. The pipeline builder uses it to reject, at creation
    /// time, a step whose required direction has no implementation, mirroring
    /// PROJ's `pj_create_pipeline` forward-path check (`src/pipeline.cpp`).
    pub fn inner_has_inverse(&self) -> bool {
        self.operation.has_inverse()
    }
    /// Whether this is a pass-through geographic (lat/long) operation.
    pub fn is_latlong(&self) -> bool {
        self.is_latlong
    }
    /// Central meridian (radians).
    pub fn lam0(&self) -> f64 {
        self.lam0
    }
    /// Latitude of origin (radians).
    pub fn phi0(&self) -> f64 {
        self.phi0
    }
    /// False easting (meters, pre-scale).
    pub fn x0(&self) -> f64 {
        self.x0
    }
    /// False northing (meters, pre-scale).
    pub fn y0(&self) -> f64 {
        self.y0
    }
    /// False up (meters, pre-scale).
    pub fn z0(&self) -> f64 {
        self.z0
    }
    /// Scale factor at the central meridian.
    pub fn k0(&self) -> f64 {
        self.k0
    }
    /// Horizontal input-unit-to-meter factor.
    pub fn to_meter(&self) -> f64 {
        self.to_meter
    }
    /// Horizontal meter-to-output-unit factor.
    pub fn fr_meter(&self) -> f64 {
        self.fr_meter
    }
    /// Vertical input-unit-to-meter factor.
    pub fn vto_meter(&self) -> f64 {
        self.vto_meter
    }
    /// Vertical meter-to-output-unit factor.
    pub fn vfr_meter(&self) -> f64 {
        self.vfr_meter
    }
    /// Prime-meridian offset from Greenwich (radians).
    pub fn from_greenwich(&self) -> f64 {
        self.from_greenwich
    }
    /// Allow longitudes outside the normal range (`+over`).
    pub fn over(&self) -> bool {
        self.over
    }
    /// Geocentric latitude flag (`+geoc`).
    pub fn geoc(&self) -> bool {
        self.geoc
    }
    /// The reference ellipsoid the boundary math is expressed against.
    pub fn ellipsoid(&self) -> Ellipsoid {
        self.ellipsoid
    }
    /// Input units of the whole operation (forward direction).
    pub fn left(&self) -> IoUnits {
        self.left
    }
    /// Output units of the whole operation (forward direction).
    pub fn right(&self) -> IoUnits {
        self.right
    }

    /// Forward transform with prepare/finalize.
    ///
    /// Ported from `src/fwd.cpp` (`fwd_prepare` -> operation -> `fwd_finalize`).
    fn forward_impl(&self, c: Coord) -> ProjResult<Coord> {
        // --- fwd_prepare ---
        let v = c.v();
        // Ported from src/fwd.cpp (`fwd_prepare`): `if (HUGE_VAL == coo.v[0] ||
        // HUGE_VAL == coo.v[1] || HUGE_VAL == coo.v[2])`. PROJ's guard is an
        // *exact* equality check against the `HUGE_VAL` (+infinity) error
        // sentinel, not a general non-finite test: NaN input is NOT caught
        // here and flows through into the operation ("NaN in -> NaN out"),
        // while a genuine +infinity input short-circuits to the error
        // coordinate. `t` (`v[3]`) is intentionally excluded, matching PROJ
        // (its own `t == HUGE_VAL` "unspecified time" sentinel is handled
        // separately downstream and must not be disturbed by this guard).
        if v[0] == f64::INFINITY || v[1] == f64::INFINITY || v[2] == f64::INFINITY {
            return Ok(Coord::error());
        }

        let prepared = if self.left == IoUnits::Radians {
            let mut lam = v[0];
            let mut phi = v[1];
            let z = v[2];
            let t = v[3];

            let tt = phi.abs() - oxiproj_core::M_HALFPI;
            if tt > PJ_EPS_LAT {
                return Err(ProjError::InvalidCoord);
            }
            // Ported from src/fwd.cpp (`fwd_prepare`): `if (coo.lp.lam > 10 ||
            // coo.lp.lam < -10) error;`. This must be written as a direct `>`
            // / `<` disjunction, NOT as a negated `RangeInclusive::contains`
            // check: for NaN, `!(-10.0..=10.0).contains(&lam)` is `true`
            // (`contains` on NaN is `false`, since NaN fails both internal
            // comparisons, so the negation flips it to `true`), which would
            // incorrectly reject a NaN longitude that PROJ's own `lam > 10 ||
            // lam < -10` lets through (both disjuncts are `false` for NaN).
            // NOT equivalent to `!(-10.0..=10.0).contains(&lam)` for NaN (see
            // the comment above) — clippy's suggested rewrite would silently
            // reintroduce the bug, so it is deliberately overridden here.
            #[allow(clippy::manual_range_contains)]
            if lam > 10.0 || lam < -10.0 {
                return Err(ProjError::InvalidCoord);
            }
            // Clamp latitude to ±π/2 (PROJ src/fwd.cpp).
            phi = phi.clamp(-oxiproj_core::M_HALFPI, oxiproj_core::M_HALFPI);

            // `+geoc`: the input latitude is geocentric, so fwd_prepare converts
            // it to geographic before the projection. PROJ src/fwd.cpp calls
            // `pj_geocentric_latitude(P, PJ_INV, ...)`, i.e. the geocentric ->
            // geographic direction (`rone_es`), NOT the forward `one_es`. The
            // near-pole guard (last cm) copies the latitude through unchanged,
            // mirroring `pj_geocentric_latitude` (`src/conversions/geoc.cpp`).
            if self.geoc && self.ellipsoid.es != 0.0 {
                let limit = oxiproj_core::M_HALFPI - 1e-9;
                if phi.abs() < limit {
                    phi = (self.ellipsoid.rone_es * phi.tan()).atan();
                }
            }
            if !self.over {
                lam = oxiproj_core::adjlon(lam);
            }
            // `+geoidgrids` vertical grid shift is applied as its own injected
            // `vgridshift` pipeline step (see `create::create_single`), before
            // the projection and thus before this central-meridian subtraction,
            // matching PROJ's `P->vgridshift` in fwd_prepare; it is not applied
            // inline here.
            lam = (lam - self.from_greenwich) - self.lam0;
            if !self.over {
                lam = oxiproj_core::adjlon(lam);
            }
            Coord::new(lam, phi, z, t)
        } else {
            c
        };

        // --- operation ---
        let mid = self.operation.forward_4d(prepared)?;

        // --- fwd_finalize ---
        let mv = mid.v();
        let mut x = mv[0];
        let mut y = mv[1];
        let mut z = mv[2];
        let t = mv[3];
        match self.right {
            IoUnits::Classic | IoUnits::Projected => {
                if self.right == IoUnits::Classic {
                    x *= self.ellipsoid.a;
                    y *= self.ellipsoid.a;
                }
                x = self.fr_meter * (x + self.x0);
                y = self.fr_meter * (y + self.y0);
                z = self.vfr_meter * (z + self.z0);
            }
            IoUnits::Radians => {
                z = self.vfr_meter * (z + self.z0);
                match self.lon_wrap_center {
                    // `+lon_wrap`: wrap into [center-pi, center+pi],
                    // unconditionally of `+over` (PROJ `fwd.cpp` has no
                    // `over` check in the `is_long_wrap_set` branch).
                    Some(center) => {
                        x = center + oxiproj_core::adjlon(x - center);
                    }
                    None => {
                        if !self.over {
                            x = oxiproj_core::adjlon(x);
                        }
                    }
                }
            }
            IoUnits::Cartesian => {
                x *= self.fr_meter;
                y *= self.fr_meter;
                z *= self.fr_meter;
                // PROJ's `is_geocent` branch runs the coordinate through an
                // internal `cart` op here; OxiProj instead folds that geodetic
                // <-> geocentric Cartesian math into the `geocent` operation
                // itself (see oxiproj-transformations `conversions/geocent.rs`),
                // which runs before this scaling, so no extra step is needed.
            }
            IoUnits::Whatever | IoUnits::Degrees => {}
        }
        // `+axis` output reordering is applied as its own injected `axisswap`
        // pipeline step after the projection (see `create::create_single`),
        // matching PROJ's `P->axisswap` in fwd_finalize; not applied inline.
        Ok(Coord::new(x, y, z, t))
    }

    /// Inverse transform with prepare/finalize.
    ///
    /// Ported from `src/inv.cpp` (`inv_prepare` -> operation -> `inv_finalize`).
    fn inverse_impl(&self, c: Coord) -> ProjResult<Coord> {
        // --- inv_prepare ---
        let v = c.v();
        // Ported from src/inv.cpp (`inv_prepare`): `if (coo.v[0] == HUGE_VAL ||
        // coo.v[1] == HUGE_VAL || coo.v[2] == HUGE_VAL)`. Like `forward_impl`,
        // this is an exact-equality check against the `HUGE_VAL` (+infinity)
        // error sentinel, not a general non-finite test — NaN input passes
        // through to the operation ("NaN in -> NaN out") instead of being
        // folded into an error here. `t` (`v[3]`) stays excluded, matching
        // PROJ (its own `t == HUGE_VAL` sentinel is handled separately and
        // must not be disturbed by this guard).
        //
        // PROJ's inv_prepare (src/inv.cpp) sets an errno but still yields the
        // error coordinate, so pj_inv/pj_inv4d degrade this single point to
        // `proj_coord_error()` rather than aborting the batch. Mirror the
        // forward path (`forward_impl`) and return the error coordinate so a
        // lone +infinity inverse input does not fail a whole array.
        if v[0] == f64::INFINITY || v[1] == f64::INFINITY || v[2] == f64::INFINITY {
            return Ok(Coord::error());
        }
        // Axisswap on the inverse path is applied as its own injected pipeline
        // step (see `create::create_single`), mirroring PROJ's `P->axisswap`
        // in inv_prepare; it is not re-applied inside this inline prepare.
        let t = v[3];
        let prepared = match self.right {
            IoUnits::Whatever | IoUnits::Degrees => c,
            IoUnits::Cartesian => {
                let x = self.to_meter * v[0];
                let y = self.to_meter * v[1];
                let z = self.to_meter * v[2];
                // PROJ's `is_geocent` branch runs the coordinate through the
                // internal `cart` inverse here; OxiProj folds that geocentric
                // Cartesian <-> geodetic math into the `geocent` operation
                // itself (see oxiproj-transformations `conversions/geocent.rs`),
                // which runs after this de-scaling, so no extra step is needed.
                Coord::new(x, y, z, t)
            }
            IoUnits::Projected | IoUnits::Classic => {
                let mut x = self.to_meter * v[0] - self.x0;
                let mut y = self.to_meter * v[1] - self.y0;
                let z = self.vto_meter * v[2] - self.z0;
                if self.right == IoUnits::Classic {
                    x *= self.ellipsoid.ra;
                    y *= self.ellipsoid.ra;
                }
                Coord::new(x, y, z, t)
            }
            IoUnits::Radians => {
                let z = self.vto_meter * v[2] - self.z0;
                Coord::new(v[0], v[1], z, t)
            }
        };

        // --- operation ---
        let mid = self.operation.inverse_4d(prepared)?;

        // --- inv_finalize ---
        if self.left == IoUnits::Radians {
            let mv = mid.v();
            let mut lam = mv[0];
            let phi = mv[1];
            let z = mv[2];
            let t = mv[3];
            lam = lam + self.from_greenwich + self.lam0;
            if !self.over {
                lam = oxiproj_core::adjlon(lam);
            }
            // `+geoc`: inv_finalize converts the geographic latitude back to
            // geocentric. PROJ src/inv.cpp calls
            // `pj_geocentric_latitude(P, PJ_FWD, ...)`, i.e. the geographic ->
            // geocentric direction (`one_es`) — the inverse of what fwd_prepare
            // applied. The near-pole guard copies the latitude through unchanged.
            let phi = if self.geoc && self.ellipsoid.es != 0.0 {
                let limit = oxiproj_core::M_HALFPI - 1e-9;
                if phi.abs() < limit {
                    (self.ellipsoid.one_es * phi.tan()).atan()
                } else {
                    phi
                }
            } else {
                phi
            };
            // `+geoidgrids` vgridshift and datum shifts are applied as their
            // own injected pipeline steps (see `create::create_single`): on the
            // inverse path the pipeline runs those steps in reverse (vgridshift
            // INVERSE before the horizontal shift), matching PROJ's
            // `P->vgridshift`/`P->helmert` order in inv_finalize.
            Ok(Coord::new(lam, phi, z, t))
        } else {
            Ok(mid)
        }
    }

    /// Apply the `+lon_wrap` output wrap (PROJ `fwd_finalize`'s
    /// `is_long_wrap_set` branch) to a coordinate whose first slot is an
    /// angular longitude, when this operation both carries a wrap center and
    /// has angular (`Radians`) forward-sense output. A no-op otherwise,
    /// including for a non-finite (error) longitude — mirrors PROJ's
    /// `if (coo.lpz.lam != HUGE_VAL)` guard around the wrap.
    fn apply_lon_wrap(&self, c: Coord) -> Coord {
        let Some(center) = self.lon_wrap_center else {
            return c;
        };
        if self.right != IoUnits::Radians {
            return c;
        }
        let v = c.v();
        if !v[0].is_finite() {
            return c;
        }
        Coord::new(
            center + oxiproj_core::adjlon(v[0] - center),
            v[1],
            v[2],
            v[3],
        )
    }

    /// Dispatch a bypassed (prepare/finalize-skipping) operation call.
    ///
    /// `run_inner_inverse` selects which of the inner operation's own
    /// forward/inverse maps actually runs. PROJ's `proj_trans` flips
    /// `PJ_FWD`/`PJ_INV` at the top level when `P->inverted`, and only the
    /// `PJ_FWD` path threads through `fwd_finalize`'s `is_long_wrap_set`
    /// branch (`inv_finalize` never wraps) — so the `+lon_wrap` output wrap
    /// is applied exactly when `run_inner_inverse` is `false`.
    fn bypass_dispatch(&self, c: Coord, run_inner_inverse: bool) -> ProjResult<Coord> {
        if run_inner_inverse {
            self.operation.inverse_4d(c)
        } else {
            let out = self.operation.forward_4d(c)?;
            Ok(self.apply_lon_wrap(out))
        }
    }

    /// Run the operation in the forward direction, honoring `inverted`.
    pub fn forward(&self, c: Coord) -> ProjResult<Coord> {
        if self.bypass_prepare_finalize {
            return self.bypass_dispatch(c, self.inverted);
        }
        if self.inverted {
            self.inverse_impl(c)
        } else {
            self.forward_impl(c)
        }
    }

    /// Run the operation in the inverse direction, honoring `inverted`.
    pub fn inverse(&self, c: Coord) -> ProjResult<Coord> {
        if self.bypass_prepare_finalize {
            return self.bypass_dispatch(c, !self.inverted);
        }
        if self.inverted {
            self.forward_impl(c)
        } else {
            self.inverse_impl(c)
        }
    }

    /// Compute map-scale factors (Tissot indicatrix) at a geographic coordinate.
    ///
    /// `coord` must be in geographic input units — radians if `self.left == IoUnits::Radians`,
    /// degrees if `self.left == IoUnits::Degrees`. The coordinate layout is
    /// `Coord::new(longitude, latitude, z, t)`.
    ///
    /// Returns the [`oxiproj_core::Factors`] struct describing Tissot indicatrix parameters
    /// at the given point.
    ///
    /// Ported from PROJ 9.8.0 `src/factors.cpp` (`proj_factors` → `pj_factors` → `pj_deriv`).
    pub fn factors(&self, coord: Coord) -> ProjResult<oxiproj_core::Factors> {
        let v = coord.v();
        // Geographic input must be angular. PROJ's `proj_factors`
        // (src/factors.cpp) only handles the RADIANS/DEGREES boundary units;
        // for any other left unit (Cartesian/Whatever — e.g. a pipeline `Pj` or
        // a transformation) the coordinate is not a lon/lat pair and has no
        // meaningful map-scale factors, so report the operation as unsupported
        // rather than silently misinterpreting the input as degrees.
        let (mut phi, lam_abs) = match self.left {
            IoUnits::Radians => (v[1], v[0]),
            IoUnits::Degrees => (v[1].to_radians(), v[0].to_radians()),
            _ => return Err(ProjError::UnsupportedOperation),
        };

        // Effective factors eccentricity: `0` for sphere-only projections
        // (PROJ forces `P->es = 0` in their setup), else the ellipsoid's `es`.
        let es = self.factors_es;
        let one_es = 1.0 - es;

        // If the latitude is geocentric (`+geoc`), convert to geographic before
        // taking the derivative, mirroring `fwd_prepare` and PROJ `pj_factors`
        // (`if (internal->geoc) lp = pj_geocentric_latitude(internal, PJ_INV, ...)`,
        // src/factors.cpp:50): PJ_INV is the geocentric -> geographic direction,
        // i.e. `rone_es`, not the forward `one_es`.
        if self.geoc && es != 0.0 {
            phi = (self.ellipsoid.rone_es * phi.tan()).atan();
        }

        // Longitude relative to central meridian (as done in fwd_prepare)
        let lam = lam_abs - self.from_greenwich - self.lam0;

        // The closure calls the raw inner operation's forward in normalized space.
        let op = &*self.operation;

        oxiproj_core::Factors::compute(phi, lam, es, one_es, |lam_f, phi_f| {
            let c = Coord::new(lam_f, phi_f, 0.0, 0.0);
            let r = op.forward_4d(c)?;
            let rv = r.v();
            Ok((rv[0], rv[1]))
        })
    }

    /// Compute exact distortion factors via automatic differentiation.
    ///
    /// Returns `Err(ProjError::UnsupportedOperation)` when:
    /// - `ad_proj` is `None` (no AD-capable projection registered), or
    /// - the projection's `project_fwd_generic` returns `UnsupportedOperation`.
    ///
    /// `coord` must be in the same geographic input units as for [`Pj::factors`].
    pub fn factors_exact(&self, coord: Coord) -> ProjResult<oxiproj_core::FactorsExact> {
        use oxiproj_core::autodiff::Dual1;

        let ad = match &self.ad_proj {
            None => return Err(ProjError::UnsupportedOperation),
            Some(p) => p,
        };

        let v = coord.v();
        // Same angular-input contract as [`Pj::factors`]: only RADIANS/DEGREES
        // boundary units carry lon/lat; reject anything else.
        let (mut phi, lam_abs) = match self.left {
            IoUnits::Radians => (v[1], v[0]),
            IoUnits::Degrees => (v[1].to_radians(), v[0].to_radians()),
            _ => return Err(ProjError::UnsupportedOperation),
        };
        // Effective factors eccentricity: `0` for sphere-only projections
        // (PROJ forces `P->es = 0` in their setup), else the ellipsoid's `es`.
        let es = self.factors_es;
        // Geocentric-latitude (`+geoc`) adjustment, mirroring `fwd_prepare` /
        // PROJ `pj_factors` (PJ_INV, geocentric -> geographic, `rone_es`).
        if self.geoc && es != 0.0 {
            phi = (self.ellipsoid.rone_es * phi.tan()).atan();
        }
        let lam = lam_abs - self.from_greenwich - self.lam0;

        let lam_d = Dual1::<2>::variable(lam, 0);
        let phi_d = Dual1::<2>::variable(phi, 1);
        let (x, y) = ad.project_fwd_dual2(lam_d, phi_d)?;

        // The `ProjectGeneric` forward runs in the engine's normalised (a = 1)
        // space — the semi-major-axis scaling is applied later by
        // `fwd_finalize`, not by `project_fwd_generic`. So the Jacobian is the
        // normalised-space Jacobian and MUST be reduced with `a = 1.0` (passing
        // `self.ellipsoid.a` would wrongly shrink every magnitude by the radius
        // whenever a ≠ 1). The squared eccentricity threads the ellipsoid
        // radius-of-curvature correction (PROJ `factors.cpp`); `es == 0`
        // recovers the spherical factors exactly.
        Ok(oxiproj_core::FactorsExact::from_jacobian_es(
            x.d[0], x.d[1], y.d[0], y.d[1], phi, 1.0, es,
        ))
    }
}

impl Operation for Pj {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        self.forward(c)
    }
    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        self.inverse(c)
    }
    fn forward_2d(&self, lp: Lp) -> ProjResult<Xy> {
        let c = self.forward(Coord::new(lp.lam, lp.phi, 0.0, 0.0))?;
        Ok(c.xy())
    }
    fn inverse_2d(&self, xy: Xy) -> ProjResult<Lp> {
        let c = self.inverse(Coord::new(xy.x, xy.y, 0.0, 0.0))?;
        Ok(c.lp())
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[derive(Debug)]
    struct Id;
    impl Operation for Id {
        fn forward_2d(&self, lp: Lp) -> ProjResult<Xy> {
            Ok(Xy::new(lp.lam, lp.phi))
        }
        fn inverse_2d(&self, xy: Xy) -> ProjResult<Lp> {
            Ok(Lp::new(xy.x, xy.y))
        }
    }

    fn unit_latlong() -> Pj {
        let ellipsoid = Ellipsoid::named("WGS84").unwrap();
        Pj {
            operation: Box::new(Id),
            ellipsoid,
            factors_es: ellipsoid.es,
            lam0: 0.0,
            phi0: 0.0,
            x0: 0.0,
            y0: 0.0,
            z0: 0.0,
            k0: 1.0,
            to_meter: 1.0,
            fr_meter: 1.0,
            vto_meter: 1.0,
            vfr_meter: 1.0,
            from_greenwich: 0.0,
            over: false,
            geoc: false,
            lon_wrap_center: None,
            is_latlong: true,
            left: IoUnits::Radians,
            right: IoUnits::Radians,
            inverted: false,
            bypass_prepare_finalize: false,
            omit_fwd: false,
            omit_inv: false,
            ad_proj: None,
            op_name: String::new(),
            axisswap_order: None,
        }
    }

    #[test]
    fn boundary_accessors_expose_fields() {
        // Regression for the Pj field-hygiene change: `ellipsoid`, `left` and
        // `right` are now `pub(crate)` and reachable from outside the crate
        // only through these read accessors, which must return the field values.
        let pj = unit_latlong();
        assert_eq!(pj.left(), IoUnits::Radians);
        assert_eq!(pj.right(), IoUnits::Radians);
        let ell = pj.ellipsoid();
        let named = Ellipsoid::named("WGS84").unwrap();
        assert_eq!(ell.a, named.a);
        assert_eq!(ell.es, named.es);
    }

    #[test]
    fn latlong_round_trip() {
        let pj = unit_latlong();
        let lam = 0.2;
        let phi = 0.5;
        let fwd = pj.forward(Coord::new(lam, phi, 0.0, 0.0)).unwrap();
        let inv = pj.inverse(fwd).unwrap();
        let r = inv.v();
        assert!((r[0] - lam).abs() < 1e-12);
        assert!((r[1] - phi).abs() < 1e-12);
    }

    #[test]
    fn non_finite_forward_yields_error_coord() {
        let pj = unit_latlong();
        let out = pj.forward(Coord::new(f64::NAN, 0.0, 0.0, 0.0)).unwrap();
        // `Coord::is_error()` flags any non-finite component, so it stays
        // `true` here regardless of *why* `v[0]` is non-finite — but that
        // alone doesn't distinguish "the HUGE_VAL error sentinel" from "a
        // genuine NaN that flowed through the transform". Pin the latter
        // explicitly: see `nan_forward_propagates_through_not_folded_to_huge_val`.
        assert!(out.is_error());
    }

    /// Regression test: `forward_impl`'s non-finite guard is ported from
    /// PROJ's `fwd_prepare` `if (HUGE_VAL == coo.v[0] || ...)`, an *exact*
    /// equality check against the `+infinity` error sentinel — not a general
    /// `!is_finite()` test. A NaN `lam` must flow through the transform
    /// ("NaN in -> NaN out") rather than being folded into `Coord::error()`
    /// (all four slots `+infinity`); only a genuine `+infinity` input should
    /// short-circuit to the error coordinate. Distinguishing the two matters
    /// because `is_error()` alone can't tell them apart (both have a
    /// non-finite component), and PROJ's own gie corpus has cases (e.g.
    /// `+proj=laea +lat_0=90 ... accept NaN NaN NaN NaN expect NaN NaN NaN
    /// NaN`) where a NaN input must produce a NaN *result*, not an error.
    #[test]
    fn nan_forward_propagates_through_not_folded_to_huge_val() {
        let pj = unit_latlong();
        let out = pj.forward(Coord::new(f64::NAN, 0.0, 0.0, 0.0)).unwrap();
        let v = out.v();
        assert!(v[0].is_nan(), "lam should be NaN, got {}", v[0]);
        assert_eq!(v[1], 0.0, "finite phi should pass through unchanged");
        assert_ne!(
            v[1],
            f64::INFINITY,
            "must not have collapsed to the Coord::error() sentinel"
        );

        // A genuine +infinity input, by contrast, must still hit the HUGE_VAL
        // guard and degrade to the error coordinate.
        let inf_out = pj
            .forward(Coord::new(f64::INFINITY, 0.0, 0.0, 0.0))
            .unwrap();
        assert!(inf_out.is_error());
        assert_eq!(inf_out.v(), [f64::INFINITY; 4]);
    }

    #[test]
    fn out_of_domain_latitude_rejected() {
        let pj = unit_latlong();
        let res = pj.forward(Coord::new(0.0, 2.0, 0.0, 0.0));
        assert!(res.is_err());
    }

    #[test]
    fn factors_exact_spherical_mercator_analytic() {
        use oxiproj_core::autodiff::Dual1;

        // Spherical Mercator Jacobian at lam=0.3, phi=0.5:
        //   x = k0*lam  => dx/dlam = k0, dx/dphi = 0
        //   y = k0*ln(tan(pi/4 + phi/2)) = k0*asinh(tan(phi))
        //      => dy/dlam = 0, dy/dphi = k0/cos(phi) = k0*sec(phi)
        let phi0 = 0.5_f64;
        let lam_d = Dual1::<2>::variable(0.3, 0);
        let phi_d = Dual1::<2>::variable(phi0, 1);

        let k0 = 1.0_f64;
        // x = k0 * lam
        let x = lam_d * Dual1::constant(k0);
        // y = k0 * asinh(tan(phi)) = k0 * ln(tan(phi) + sqrt(tan^2(phi) + 1))
        let t = phi_d.tan();
        let y = (t + (t * t + Dual1::constant(1.0)).sqrt()).ln() * Dual1::constant(k0);

        // Verify Jacobian entries
        assert!((x.d[0] - k0).abs() < 1e-12, "dx/dlam = k0");
        assert!(x.d[1].abs() < 1e-12, "dx/dphi = 0");
        assert!(y.d[0].abs() < 1e-12, "dy/dlam = 0");
        let sec_phi = 1.0 / phi0.cos();
        assert!(
            (y.d[1] - sec_phi * k0).abs() < 1e-9,
            "dy/dphi = sec(phi)*k0"
        );

        // FactorsExact: for spherical Mercator, h = k = sec(phi), omega = 0.
        // Mercator is conformal (h=k, omega=0) but NOT equal-area.
        // s = |det J| / (a^2 * cos phi) = sec(phi) / cos(phi) = sec^2(phi).
        let fe =
            oxiproj_core::FactorsExact::from_jacobian(x.d[0], x.d[1], y.d[0], y.d[1], phi0, 1.0);
        assert!(
            (fe.meridian_scale - sec_phi).abs() < 1e-9,
            "h={}",
            fe.meridian_scale
        );
        assert!(
            (fe.parallel_scale - sec_phi).abs() < 1e-9,
            "k={}",
            fe.parallel_scale
        );
        let expected_s = sec_phi * sec_phi;
        assert!(
            (fe.areal_scale - expected_s).abs() < 1e-9,
            "s={}",
            fe.areal_scale
        );
        assert!(
            fe.angular_distortion.abs() < 1e-9,
            "omega={}",
            fe.angular_distortion
        );
        assert!(fe.is_exact);
    }

    #[test]
    fn factors_exact_no_ad_proj_returns_error() {
        let pj = unit_latlong();
        let res = pj.factors_exact(Coord::new(0.1, 0.3, 0.0, 0.0));
        assert!(
            matches!(res, Err(oxiproj_core::ProjError::UnsupportedOperation)),
            "expected UnsupportedOperation, got {res:?}"
        );
    }

    #[test]
    fn non_finite_inverse_yields_error_coord() {
        // E2 finding 3: a lone non-finite inverse input must not bubble up as
        // an `Err` and fail a whole batch (matching PROJ's per-point graceful
        // degradation). It resolves as `Ok`, with the *coordinate* content
        // depending on why it's non-finite: a genuine `+infinity` (PROJ's
        // `HUGE_VAL` sentinel) degrades to `Coord::error()`, while NaN flows
        // through the transform instead — see
        // `nan_inverse_propagates_through_not_folded_to_huge_val` below for
        // that distinction pinned explicitly.
        let pj = unit_latlong();
        let out = pj.inverse(Coord::new(f64::NAN, 0.0, 0.0, 0.0)).unwrap();
        assert!(out.is_error());
        let out2 = pj
            .inverse(Coord::new(0.0, f64::INFINITY, 0.0, 0.0))
            .unwrap();
        assert!(out2.is_error());
    }

    /// Regression test: `inverse_impl`'s non-finite guard, like
    /// `forward_impl`'s, is an *exact* equality check against the `HUGE_VAL`
    /// (+infinity) error sentinel (ported from PROJ's `inv_prepare`), not a
    /// general `!is_finite()` test. See
    /// `nan_forward_propagates_through_not_folded_to_huge_val` for the
    /// rationale.
    #[test]
    fn nan_inverse_propagates_through_not_folded_to_huge_val() {
        let pj = unit_latlong();
        let out = pj.inverse(Coord::new(f64::NAN, 0.0, 0.0, 0.0)).unwrap();
        let v = out.v();
        assert!(v[0].is_nan(), "lam should be NaN, got {}", v[0]);
        assert_eq!(v[1], 0.0, "finite phi should pass through unchanged");
        assert_ne!(
            v[1],
            f64::INFINITY,
            "must not have collapsed to the Coord::error() sentinel"
        );
    }

    #[test]
    fn factors_rejects_non_angular_left() {
        // E2 finding 4: factors on a non-angular boundary op (Cartesian /
        // Whatever) must report UnsupportedOperation, not silently reinterpret
        // the coordinate as degrees.
        let mut pj = unit_latlong();
        pj.left = IoUnits::Cartesian;
        assert!(matches!(
            pj.factors(Coord::new(0.1, 0.3, 0.0, 0.0)),
            Err(ProjError::UnsupportedOperation)
        ));
        assert!(matches!(
            pj.factors_exact(Coord::new(0.1, 0.3, 0.0, 0.0)),
            Err(ProjError::UnsupportedOperation)
        ));
        pj.left = IoUnits::Whatever;
        assert!(matches!(
            pj.factors(Coord::new(0.1, 0.3, 0.0, 0.0)),
            Err(ProjError::UnsupportedOperation)
        ));
    }

    #[test]
    fn factors_applies_geoc_latitude_adjustment() {
        // E2 finding 4: with `+geoc` the latitude is converted geocentric ->
        // geographic before the derivative (mirroring fwd_prepare / PROJ
        // pj_factors), so the parallel scale k of the identity op is evaluated
        // at the adjusted latitude.
        let phi_in = 0.5_f64;
        let lam_in = 0.2_f64;

        let mut geo = unit_latlong();
        geo.geoc = true;
        let f_geoc = geo.factors(Coord::new(lam_in, phi_in, 0.0, 0.0)).unwrap();

        let plain = unit_latlong(); // geoc = false
        let f_plain = plain.factors(Coord::new(lam_in, phi_in, 0.0, 0.0)).unwrap();

        // For the identity op, k = sqrt(1 - es*sin^2 phi)/cos phi at the
        // (possibly adjusted) latitude. Compute the expected geoc value directly.
        // `+geoc` converts geocentric -> geographic (PJ_INV, `rone_es`), matching
        // PROJ `pj_factors`/`fwd_prepare`.
        let ell = &geo.ellipsoid;
        let phi_adj = (ell.rone_es * phi_in.tan()).atan();
        let expected_k = (1.0 - ell.es * phi_adj.sin().powi(2)).sqrt() / phi_adj.cos();
        assert!(
            (f_geoc.k - expected_k).abs() < 1e-6,
            "geoc k: got {}, want {}",
            f_geoc.k,
            expected_k
        );
        // The adjustment moves the latitude, so geoc and non-geoc k differ.
        assert!(
            (f_geoc.k - f_plain.k).abs() > 1e-4,
            "geoc adjustment must change k: geoc={}, plain={}",
            f_geoc.k,
            f_plain.k
        );
    }
}