oxiproj_engine/

create.rs

//! Top-level construction and transformation entry points.
//!
//! Ported from PROJ 9.8.0 `src/4D_api.cpp` (`proj_create`, `proj_trans`,
//! `proj_trans_array`) and `src/pipeline.cpp` (pipeline assembly).

use crate::context::Context;
use crate::params::{parse, ParamList};
use crate::pipeline::Pipeline;
use crate::pj::Pj;
use oxiproj_core::{Coord, Direction, Ellipsoid, IoUnits, ProjError, ProjResult};

/// The keys that, if present, mean a parameter set carries its own ellipsoid.
const ELLIPSOID_KEYS: [&str; 9] = ["R", "a", "ellps", "datum", "rf", "f", "es", "e", "b"];

fn has_ellipsoid_def(pl: &ParamList) -> bool {
    ELLIPSOID_KEYS.iter().any(|k| pl.exists(k))
}

// Pipeline boundary i/o units mirror PROJ 9.8.0 src/internal.cpp
// `pj_left`/`pj_right`, consumed by src/pipeline.cpp
// (`P->left = pj_left(front); P->right = pj_right(back)`). A step that runs
// inverted swaps its raw left/right when reported in forward sense, and the
// `Classic` sentinel is normalized to `Projected`.

/// Effective forward-input units of a single step, mirroring PROJ
/// `pj_left` (src/internal.cpp): an inverted step swaps its raw left/right,
/// and `Classic` is normalized to `Projected`.
fn effective_left(s: &Pj) -> IoUnits {
    let u = if s.inverted { s.right } else { s.left };
    match u {
        IoUnits::Classic => IoUnits::Projected,
        other => other,
    }
}

/// Effective forward-output units of a single step, mirroring PROJ
/// `pj_right` (src/internal.cpp): an inverted step swaps its raw left/right,
/// and `Classic` is normalized to `Projected`.
fn effective_right(s: &Pj) -> IoUnits {
    let u = if s.inverted { s.left } else { s.right };
    match u {
        IoUnits::Classic => IoUnits::Projected,
        other => other,
    }
}

/// Propagate `IoUnits::Whatever` boundary units of unit-agnostic steps
/// (e.g. `axisswap`, `push`, `pop`, `affine`) from their neighbours, mirroring
/// PROJ 9.8.0 `src/pipeline.cpp` (the two passes right after the steps are
/// built, before `P->left = pj_left(front); P->right = pj_right(back)`).
///
/// A step whose *effective* left and right are BOTH `Whatever` adopts the units
/// of an adjacent step so that angular pipelines whose last (or first) step is
/// unit-agnostic still report `Radians` at the pipeline boundary — which is what
/// drives the degree<->radian conversion in the CLI / gie harness
/// (`proj_angular_input`/`proj_angular_output`, i.e. `pj_left == RADIANS` /
/// `pj_right == RADIANS`). Without this, a `+proj=pipeline +step +proj=latlong
/// +step +proj=axisswap` pipeline reports `Whatever` output and its degrees are
/// never emitted (radians leak out); likewise a leading `push`/`pop` step would
/// leave the pipeline input `Whatever`, so degrees are fed to a following
/// `utm`/`merc` step verbatim as radians and blow up (out-of-domain / NaN).
///
/// PROJ reads neighbour units through `pj_left`/`pj_right` (the effective,
/// inversion-aware, `Classic`->`Projected`-normalized units — our
/// [`effective_left`]/[`effective_right`]) and writes the adopted value into the
/// step's *raw* `left` and `right` fields, setting BOTH to the same value so the
/// step is thereafter unit-consistent regardless of its `inverted` flag.
///
/// Pass 1 walks right-to-left: a `Whatever/Whatever` step takes its right
/// neighbour's `pj_left`, provided that neighbour is not itself fully
/// `Whatever`. Pass 2 walks left-to-right using the left neighbour's
/// `pj_right`. Two passes let a run of unit-agnostic steps fill in from whichever
/// side carries real units.
fn propagate_whatever_units(steps: &mut [Pj]) {
    let nsteps = steps.len();
    if nsteps < 2 {
        return;
    }

    // Pass 1: right-to-left (PROJ `for (i = nsteps - 2; i >= 0; --i)`).
    for i in (0..nsteps - 1).rev() {
        if effective_left(&steps[i]) == IoUnits::Whatever
            && effective_right(&steps[i]) == IoUnits::Whatever
        {
            let right_left = effective_left(&steps[i + 1]);
            let right_right = effective_right(&steps[i + 1]);
            if right_left != right_right || right_left != IoUnits::Whatever {
                steps[i].left = right_left;
                steps[i].right = right_left;
            }
        }
    }

    // Pass 2: left-to-right (PROJ `for (i = 1; i < nsteps; i++)`).
    for i in 1..nsteps {
        if effective_left(&steps[i]) == IoUnits::Whatever
            && effective_right(&steps[i]) == IoUnits::Whatever
        {
            let left_left = effective_left(&steps[i - 1]);
            let left_right = effective_right(&steps[i - 1]);
            if left_left != left_right || left_right != IoUnits::Whatever {
                steps[i].left = left_right;
                steps[i].right = left_right;
            }
        }
    }
}

/// Construct a [`Pj`] from a proj-string.
///
/// Handles both single operations and `+proj=pipeline` definitions. Pipelines
/// resolve a global ellipsoid (GRS80 by default) inherited by steps lacking an
/// explicit ellipsoid, and are driven with prepare/finalize bypassed so each
/// step does its own unit handling.
pub fn create(proj_string: &str) -> ProjResult<Pj> {
    let ctx = Context::new();
    let params = parse(proj_string);

    let is_pipeline = params
        .entries
        .iter()
        .any(|(k, v)| k == "proj" && v.as_deref() == Some("pipeline"));

    if is_pipeline {
        create_pipeline(&params, &ctx)
    } else {
        create_single(&params, &ctx)
    }
}

/// Construct a [`Pj`] from a proj-string using a supplied [`Context`].
///
/// Unlike [`create`], this allows the caller to register grid data in the
/// context before construction so that `+nadgrids=` and `+hgridshift` steps
/// can resolve their grids.
pub fn create_with_ctx(proj_string: &str, ctx: &Context) -> ProjResult<Pj> {
    let params = parse(proj_string);

    let is_pipeline = params
        .entries
        .iter()
        .any(|(k, v)| k == "proj" && v.as_deref() == Some("pipeline"));

    if is_pipeline {
        create_pipeline(&params, ctx)
    } else {
        create_single(&params, ctx)
    }
}

/// Construct a [`Pj`] that performs the **map projection only**, ignoring any
/// `+towgs84`/`+nadgrids`/`+geoidgrids` datum-shift directives.
///
/// This mirrors PROJ's legacy `proj` program (`src/apps/proj.cpp`), which drives
/// the coordinate through `pj_fwd`/`pj_inv` and runs the map projection with its
/// own unit/axis handling but **never** applies the datum shift — in contrast to
/// the modern `proj_create` API (exposed by [`create`]), which synthesizes a
/// cs2cs-style `cart`/`helmert`/`cart` pipeline whenever `+towgs84` (or a
/// `+datum` carrying one) is present.
///
/// For a `+proj=pipeline` spec, or any spec without datum-shift parameters, this
/// is identical to [`create`]; only the towgs84/nadgrids/geoidgrids injection is
/// dropped. Prime-meridian, false easting/northing, unit and `+axis` handling
/// are all preserved.
pub fn create_projection_only(proj_string: &str) -> ProjResult<Pj> {
    let ctx = Context::new();
    let params = parse(proj_string);

    let is_pipeline = params
        .entries
        .iter()
        .any(|(k, v)| k == "proj" && v.as_deref() == Some("pipeline"));

    if is_pipeline {
        create_pipeline(&params, &ctx)
    } else {
        create_single_projection_only(&params, &ctx)
    }
}

/// Projection-only single-operation build: like [`create_single`] but never
/// synthesizes the `+towgs84`/`+nadgrids`/`+geoidgrids` datum-shift steps.
/// The `+axis` re-ordering (a unit/axis concern, applied by PROJ's
/// `fwd_finalize` regardless of datum shift) is still honored.
fn create_single_projection_only(params: &ParamList, ctx: &Context) -> ProjResult<Pj> {
    // `+lon_wrap` is parsed and validated by `registry::build_single_op`
    // (which every path below eventually reaches), matching PROJ's
    // `is_long_wrap_set`/`long_wrap_center` (`src/init.cpp`); no rejection
    // needed here.

    // Classic `+axis=` re-orders / sign-flips output axes, mirroring PROJ's
    // fwd_finalize (`enu` is the identity and needs no axisswap step). The
    // `axisswap` conversion consumes its own `+axis=` (named-axis form), so
    // `+axis=` is never the classic directive for it.
    let name_is_axisswap = params
        .entries
        .iter()
        .any(|(k, v)| k == "proj" && v.as_deref() == Some("axisswap"));
    let axis = if name_is_axisswap {
        None
    } else {
        params
            .get_str("axis")
            .map(|s| s.trim().to_string())
            .filter(|s| !s.is_empty() && s != "enu")
    };

    // No axis synthesis: the projection alone, with the towgs84/nadgrids/
    // geoidgrids directives simply ignored (they never reach `build_single_op`).
    if axis.is_none() {
        return create_single_core(params, ctx);
    }

    let src_ell = crate::setup::setup_ellipsoid(params)?;
    let mut pipeline_str = String::from("+proj=pipeline");
    if params
        .entries
        .iter()
        .any(|(k, v)| (k == "inv" || k == "inverted") && v.is_none())
    {
        pipeline_str.push_str(" +inv");
    }
    pipeline_str.push_str(" +step ");
    pipeline_str.push_str(&build_projection_step(params, &src_ell));
    if let Some(ax) = &axis {
        pipeline_str.push_str(" +step +proj=axisswap +axis=");
        pipeline_str.push_str(ax);
    }

    let pipeline_params = parse(&pipeline_str);
    create_pipeline(&pipeline_params, ctx)
}

fn extract_towgs84_from_params(params: &ParamList) -> Option<String> {
    // 1. Check for explicit +towgs84=... param
    if let Some(s) = params.get_str("towgs84") {
        if !s.is_empty() {
            return Some(s.to_string());
        }
    }
    // 2. Check +datum=<name> and look up in datum table
    if let Some(datum_name) = params.get_str("datum") {
        if !datum_name.is_empty() {
            if let Some(d) = oxiproj_core::find_datum(datum_name) {
                if let Some(rest) = d.defn.strip_prefix("towgs84=") {
                    return Some(rest.to_string());
                }
            }
        }
    }
    None
}

/// Pad a `+towgs84` value list to PROJ's canonical arity.
///
/// PROJ (`src/datum_set.cpp`) zero-fills the seven Helmert parameters: a list
/// with 4-6 values keeps its rotation/scale terms (zero-padded) rather than
/// being silently truncated to a 3-parameter shift, a list with more than 7
/// values is truncated to 7, and short lists (<=3) stay a 3-parameter shift
/// (padded to 3). The result is a comma-joined string that `+proj=helmert
/// +towgs84=` accepts directly.
fn normalize_towgs84(s: &str) -> String {
    let mut vals: Vec<String> = s.split(',').map(|v| v.trim().to_string()).collect();
    if vals.len() > 7 {
        vals.truncate(7);
    } else if vals.len() > 3 {
        while vals.len() < 7 {
            vals.push("0".to_string());
        }
    } else {
        while vals.len() < 3 {
            vals.push("0".to_string());
        }
    }
    vals.join(",")
}

/// Whether a comma-joined `+towgs84` value list is entirely zero (identity).
fn towgs84_is_identity(vals: &str) -> bool {
    vals.split(',')
        .all(|v| matches!(v.trim(), "0" | "0.0" | "-0" | "-0.0" | ""))
}

/// Whether an ellipsoid is (numerically) WGS84, matching the tolerance PROJ
/// uses in `cs2cs_emulation_setup` (`src/create.cpp`) to decide whether an
/// all-zero shift still needs a cartesian ellipsoid round-trip.
fn is_wgs84_ellipsoid(ell: &Ellipsoid) -> bool {
    (ell.a - 6378137.0).abs() < 1e-8 && (ell.es - 0.0066943799901413).abs() < 1e-15
}

/// Keys stripped from the projection step of an injected pipeline: the
/// datum-shift/axis/geoid directives (applied as their own steps), inversion
/// flags (applied at the pipeline level), and every ellipsoid-defining key
/// (re-emitted as explicit `+a`/`+es`).
///
/// `lon_wrap` is deliberately NOT stripped: PROJ applies the `+lon_wrap`
/// output wrap in `fwd_finalize`, before the `+axis` axisswap step
/// (`src/fwd.cpp`: `is_long_wrap_set` handling precedes `if (P->axisswap)`),
/// so it must stay on the projection step itself rather than being dropped
/// or hoisted onto the pipeline as a whole.
fn strip_from_projection_step(key: &str) -> bool {
    matches!(
        key,
        "towgs84" | "nadgrids" | "geoidgrids" | "axis" | "inv" | "inverted"
    ) || ELLIPSOID_KEYS.contains(&key)
}

/// Build the projection step string for an injected pipeline: the original
/// parameters minus the stripped keys, with the resolved source ellipsoid
/// re-emitted as explicit `+a`/`+es` so the projection runs on exactly the
/// ellipsoid PROJ uses for `P->fwd` (`P->a`/`P->es`), including datum-derived
/// ones.
fn build_projection_step(params: &ParamList, src_ell: &Ellipsoid) -> String {
    let mut out: String = params
        .entries
        .iter()
        .filter(|(k, _)| !strip_from_projection_step(k))
        .map(|(k, v)| match v {
            Some(val) => format!("+{}={}", k, val),
            None => format!("+{}", k),
        })
        .collect::<Vec<_>>()
        .join(" ");
    out.push_str(&format!(" +a={} +es={}", src_ell.a, src_ell.es));
    out
}

fn extract_nadgrids_from_params(params: &ParamList) -> Option<String> {
    if let Some(s) = params.get_str("nadgrids") {
        if !s.is_empty() {
            return Some(s.to_string());
        }
    }
    None
}

/// Extract a non-empty classic `+geoidgrids=` value.
///
/// PROJ (`cs2cs_emulation_setup`, `src/create.cpp`) turns `+geoidgrids=<grids>`
/// into a `proj=vgridshift grids=<grids>` helper stored in `P->vgridshift` and
/// applied by `fwd_prepare`/`inv_finalize` to move between geometric
/// (ellipsoidal) and orthometric heights.
fn extract_geoidgrids_from_params(params: &ParamList) -> Option<String> {
    if let Some(s) = params.get_str("geoidgrids") {
        if !s.is_empty() {
            return Some(s.to_string());
        }
    }
    None
}

fn create_single(params: &ParamList, ctx: &Context) -> ProjResult<Pj> {
    // `+lon_wrap` is parsed and validated by `registry::build_single_op` on
    // whichever step's `+proj=` operation it lands on (either the plain
    // single operation below, or the re-emitted projection step of a
    // synthesized datum-shift/axis pipeline — see
    // `strip_from_projection_step`, which deliberately keeps `lon_wrap` in
    // that step's parameters so its own `fwd_finalize` wraps the projection
    // output before any injected `+step +proj=axisswap` runs, mirroring
    // PROJ's `fwd_finalize` order (`is_long_wrap_set` wrap, then
    // `P->axisswap`; `src/fwd.cpp`).
    let name = params
        .entries
        .iter()
        .find(|(k, _)| k == "proj")
        .and_then(|(_, v)| v.as_deref())
        .ok_or(ProjError::MissingArg)?;

    // Resolve the source ellipsoid once (handles +ellps, +datum, +a/+rf, ...).
    let src_ell = crate::setup::setup_ellipsoid(params)?;

    // Build the operation (with cs2cs datum-shift / geoidgrids / axis emulation)
    // WITHOUT consuming a top-level `+inv`, then apply the inversion here.
    let mut pj = build_single_emulated(name, params, src_ell, ctx)?;
    if params
        .entries
        .iter()
        .any(|(k, v)| (k == "inv" || k == "inverted") && v.is_none())
    {
        pj.inverted = true;
    }
    Ok(pj)
}

/// Whether `name` refers to a map projection (or the longlat/latlong
/// geographic pass-through), as opposed to a transformation/conversion.
///
/// The `oxiproj-projections` crate builds the former and returns
/// [`ProjError::InvalidOp`] for the latter (see `registry::build_single_op`,
/// which routes on exactly this distinction). A projection with otherwise-bad
/// parameters still fails with a *non*-`InvalidOp` error, so it is correctly
/// classified as a projection here. Used to gate the cs2cs datum-shift
/// emulation, which PROJ applies only to geodetic-CRS (projection/longlat)
/// operations.
fn op_is_projection(name: &str, params: &ParamList, ellipsoid: &Ellipsoid) -> bool {
    let phi0 = params.get_dms("lat_0").unwrap_or(0.0);
    let k0 = params
        .get_f64("k_0")
        .or_else(|| params.get_f64("k"))
        .unwrap_or(1.0);
    let view = crate::params::ParamView(params);
    let pp = oxiproj_projections::ProjParams {
        ellipsoid,
        phi0,
        k0,
        params: &view,
    };
    !matches!(
        oxiproj_projections::build(name, &pp),
        Err(ProjError::InvalidOp)
    )
}

/// Build a single operation applying PROJ's cs2cs datum-shift / geoidgrids /
/// axis emulation, WITHOUT consuming a top-level `+inv` (the caller applies
/// inversion by setting `.inverted`).
///
/// Used both for a standalone single op ([`create_single`]) and for each
/// pipeline step ([`create_pipeline`]), so that a step carrying
/// `+towgs84` / `+datum` / `+nadgrids` / `+geoidgrids` / `+axis` gets the same
/// synthesized `cart`/`helmert`/`vgridshift`/`axisswap` sub-pipeline PROJ builds
/// per-PJ in `pj_init` (`P->helmert`, `P->cart`, `P->vgridshift`, `P->axisswap`
/// applied by `fwd_prepare`/`fwd_finalize`). Without this, a pipeline step's
/// `+towgs84` was silently ignored (the shift only ran for standalone ops).
fn build_single_emulated(
    name: &str,
    params: &ParamList,
    src_ell: Ellipsoid,
    ctx: &Context,
) -> ProjResult<Pj> {
    let towgs84 = extract_towgs84_from_params(params);
    let nadgrids = extract_nadgrids_from_params(params);
    // Classic `+geoidgrids=` attaches a vertical grid shift (geometric <->
    // orthometric height) directly to a projection, e.g.
    // `+proj=merc +geoidgrids=egm96_15.gtx`.
    let geoidgrids = extract_geoidgrids_from_params(params);
    // Classic `+axis=` re-orders / sign-flips output axes. PROJ's
    // cs2cs_emulation_setup (`src/create.cpp`) skips the axisswap synthesis
    // when the order is already "enu"; so do we.
    //
    // The `axisswap` conversion consumes its OWN `+axis=` parameter (the
    // named-axis form `+axis=neu|nue|swd`); for it, `+axis=` is NOT the classic
    // cs2cs output-reordering directive and must be left on the operation
    // rather than synthesized into an extra (and invalid, order-less) axisswap
    // step. So never treat `+axis=` as the classic directive for `axisswap`.
    let axis = if name == "axisswap" {
        None
    } else {
        params
            .get_str("axis")
            .map(|s| s.trim().to_string())
            .filter(|s| !s.is_empty() && s != "enu")
    };

    // The cs2cs datum-shift / geoidgrids / classic-axis emulation only applies
    // to map projections and the longlat/latlong geographic pass-through — the
    // geodetic-CRS operations PROJ runs it for. A transformation/conversion
    // (helmert, cart, geocent, unitconvert, push/pop, ...) consumes these
    // parameters itself (helmert's own `+towgs84`) or ignores them, and must
    // NOT be wrapped: in particular `+proj=helmert +towgs84=...` would otherwise
    // synthesize its own inner `+proj=helmert +towgs84=...` step and recurse
    // without bound. So build the plain operation for anything without a
    // datum-shift directive, or that is not a projection.
    let needs_emulation =
        towgs84.is_some() || nadgrids.is_some() || geoidgrids.is_some() || axis.is_some();
    if !needs_emulation || !op_is_projection(name, params, &src_ell) {
        return crate::registry::build_single_op(name, params, src_ell, ctx);
    }

    // Prepare-side (pre-projection) steps, in PROJ fwd_prepare order.
    //
    // Precedence: when BOTH `+nadgrids` and `+towgs84` are present, PROJ uses
    // the grid shift and ignores the Helmert. `cs2cs_emulation_setup`
    // (`src/create.cpp`) builds `P->hgridshift` first, then gates the Helmert
    // on it: `p = P->hgridshift ? nullptr : pj_param_exists(params, "towgs84")`
    // — so nadgrids wins. Check nadgrids first to match.
    let mut prepare_steps: Vec<String> = Vec::new();
    if let Some(grid) = &nadgrids {
        // fwd_prepare applies hgridshift INVERSE on the forward (source ->
        // WGS84) path; PROJ: `proj_trans(P->hgridshift, PJ_INV, coo)`.
        prepare_steps.push(format!("+proj=hgridshift +grids={grid} +inv"));
    } else if let Some(tw) = &towgs84 {
        let vals = normalize_towgs84(tw);
        if towgs84_is_identity(&vals) {
            // Identity shift: PROJ (do_cart) still performs the cartesian
            // round-trip when the ellipsoid differs from WGS84, else nothing.
            if !is_wgs84_ellipsoid(&src_ell) {
                prepare_steps.push("+proj=cart +ellps=WGS84".to_string());
                prepare_steps.push(format!(
                    "+proj=cart +inv +a={} +es={}",
                    src_ell.a, src_ell.es
                ));
            }
        } else {
            // fwd_prepare: cart(WGS84) FWD -> helmert INV -> cart(local) INV.
            prepare_steps.push("+proj=cart +ellps=WGS84".to_string());
            prepare_steps.push(format!(
                "+proj=helmert +towgs84={vals} +convention=position_vector +exact +inv"
            ));
            prepare_steps.push(format!(
                "+proj=cart +inv +a={} +es={}",
                src_ell.a, src_ell.es
            ));
        }
    }

    // `+geoidgrids=` -> vgridshift, applied FORWARD on the forward path *after*
    // the horizontal datum shift, matching PROJ fwd_prepare order
    // (`if (P->vgridshift) coo = proj_trans(P->vgridshift, PJ_FWD, coo)` runs
    // after hgridshift/helmert; `src/fwd.cpp`). Running the pipeline in reverse
    // then applies it INVERSE before the horizontal shift, matching
    // inv_finalize (`proj_trans(P->vgridshift, PJ_INV, coo)` before
    // hgridshift/helmert; `src/inv.cpp`). No `+inv` flag: forward is FWD.
    if let Some(grids) = &geoidgrids {
        prepare_steps.push(format!("+proj=vgridshift +grids={grids}"));
    }

    // If the shift collapsed to nothing and there is no axis synthesis, build
    // the plain single operation (avoids a needless one-step pipeline). No
    // inversion is applied here — the caller does that.
    if prepare_steps.is_empty() && axis.is_none() {
        return crate::registry::build_single_op(name, params, src_ell, ctx);
    }

    // Assemble: [prepare...] projection [axisswap]. `+axis` is applied last,
    // mirroring PROJ's fwd_finalize (axisswap after the projection output).
    // The synthesized pipeline is NOT inverted here; the caller applies any
    // top-level `+inv` by flipping the returned wrapper's `.inverted`.
    let mut pipeline_str = String::from("+proj=pipeline");
    for step in &prepare_steps {
        pipeline_str.push_str(" +step ");
        pipeline_str.push_str(step);
    }
    pipeline_str.push_str(" +step ");
    pipeline_str.push_str(&build_projection_step(params, &src_ell));
    if let Some(ax) = &axis {
        pipeline_str.push_str(" +step +proj=axisswap +axis=");
        pipeline_str.push_str(ax);
    }

    let pipeline_params = parse(&pipeline_str);
    create_pipeline(&pipeline_params, ctx)
}

fn create_single_core(params: &ParamList, ctx: &Context) -> ProjResult<Pj> {
    let name = params
        .entries
        .iter()
        .find(|(k, _)| k == "proj")
        .and_then(|(_, v)| v.as_deref())
        .ok_or(ProjError::MissingArg)?;

    let ellipsoid = crate::setup::setup_ellipsoid(params)?;
    let mut pj = crate::registry::build_single_op(name, params, ellipsoid, ctx)?;

    if params
        .entries
        .iter()
        .any(|(k, v)| (k == "inv" || k == "inverted") && v.is_none())
    {
        pj.inverted = true;
    }
    Ok(pj)
}

fn create_pipeline(params: &ParamList, ctx: &Context) -> ProjResult<Pj> {
    // Split entries into a global section (before the first `step`) and one
    // vec per subsequent step. The `step` markers are not retained.
    //
    // Structural validation mirrors PROJ 9.8.0 `pj_create_pipeline`
    // (`src/pipeline.cpp`, the first `for` loop): a malformed pipeline shape is
    // rejected at creation time with the malformed-pipeline error class.
    let mut global: Vec<(String, Option<String>)> = Vec::new();
    let mut steps_entries: Vec<Vec<(String, Option<String>)>> = Vec::new();
    let mut seen_pipeline = false;
    let mut seen_step = false;
    for (k, v) in &params.entries {
        if k == "step" && v.is_none() {
            if !seen_pipeline {
                // `+step` before `+proj=pipeline` (PROJ: "+step before
                // +proj=pipeline").
                return Err(ProjError::InvalidOp);
            }
            seen_step = true;
            steps_entries.push(Vec::new());
            continue;
        }
        if k == "proj" && v.as_deref() == Some("pipeline") {
            if seen_pipeline {
                // A second `+proj=pipeline` — either a duplicate in the globals
                // or a nested pipeline as a step (PROJ: "Nesting only allowed
                // when child pipelines are wrapped in '+init's").
                return Err(ProjError::InvalidOp);
            }
            seen_pipeline = true;
            // Drop the pipeline marker from globals; it is never looked up.
            continue;
        }
        if !seen_step {
            // Global (pre-first-step) section. PROJ forbids `proj=`/`o_proj=`
            // operators here ("proj= operator before first step not allowed").
            if k == "proj" || k == "o_proj" {
                return Err(ProjError::InvalidOp);
            }
            global.push((k.clone(), v.clone()));
        } else if let Some(last) = steps_entries.last_mut() {
            last.push((k.clone(), v.clone()));
        }
    }

    if !seen_pipeline {
        // No `+proj=pipeline` marker at all (PROJ: "no pipeline def").
        return Err(ProjError::InvalidOp);
    }

    let global_pl = ParamList { entries: global };

    // Global ellipsoid: explicit definition, else GRS80.
    //
    // GRS80 default per src/pipeline.cpp set_ellipsoid (a = 6378137,
    // f = 1/298.257222101).
    let global_ellipsoid = if has_ellipsoid_def(&global_pl) {
        crate::setup::setup_ellipsoid(&global_pl)?
    } else {
        Ellipsoid::from_a_rf(6378137.0, 298.257222101)?
    };

    // Global parameters inherited by every step. PROJ's pipeline constructor
    // (`src/pipeline.cpp`) appends all global args to each step's argument list
    // *after* the step-specific args, so step-local keys win first-match
    // lookups. We exclude the ellipsoid keys (already propagated via
    // `global_ellipsoid`, which also carries the GRS80 pipeline default),
    // inversion flags (handled by `top_inverted` / per-step `inv`), the
    // per-step `omit_*` flags, and the pipeline `proj` marker.
    let appended_globals: Vec<(String, Option<String>)> = global_pl
        .entries
        .iter()
        .filter(|(k, _)| {
            k != "inv"
                && k != "inverted"
                && k != "step"
                && k != "proj"
                && k != "omit_fwd"
                && k != "omit_inv"
                && !ELLIPSOID_KEYS.contains(&k.as_str())
        })
        .cloned()
        .collect();

    let mut steps: Vec<Pj> = Vec::new();
    for step_vec in steps_entries {
        let step_pl = ParamList { entries: step_vec };
        let name = step_pl.get_str("proj").ok_or(ProjError::MissingArg)?;
        // get_str returns "" for a bare `proj`; that is still missing.
        if name.is_empty() {
            return Err(ProjError::MissingArg);
        }
        let inverted = step_pl
            .entries
            .iter()
            .any(|(k, v)| k == "inv" && v.is_none());
        let omit_fwd = step_pl
            .entries
            .iter()
            .any(|(k, v)| k == "omit_fwd" && v.is_none());
        let omit_inv = step_pl
            .entries
            .iter()
            .any(|(k, v)| k == "omit_inv" && v.is_none());
        let step_ellipsoid = if has_ellipsoid_def(&step_pl) {
            crate::setup::setup_ellipsoid(&step_pl)?
        } else {
            global_ellipsoid
        };
        // Append the inherited globals after the step-local entries so
        // step-local keys take precedence on first-match lookups.
        let mut augmented = step_pl.entries.clone();
        augmented.extend(appended_globals.iter().cloned());
        let augmented_pl = ParamList { entries: augmented };
        // Build the step WITH cs2cs datum-shift / geoidgrids / axis emulation
        // (via `build_single_emulated`), so a step carrying `+towgs84` /
        // `+datum` / `+nadgrids` / `+geoidgrids` / `+axis` gets the same
        // synthesized sub-pipeline PROJ applies per-PJ — matching PROJ's
        // per-step `pj_init` (`P->helmert`/`P->cart`/`P->vgridshift`). The step
        // inversion / omit flags are applied to the (possibly wrapper) result.
        let mut step_pj = build_single_emulated(name, &augmented_pl, step_ellipsoid, ctx)?;
        step_pj.inverted = inverted;
        step_pj.omit_fwd = omit_fwd;
        step_pj.omit_inv = omit_inv;
        steps.push(step_pj);
    }

    if steps.is_empty() {
        return Err(ProjError::MissingArg);
    }

    // Apply pipeline optimization: eliminate redundant/cancelling steps.
    // Run after the empty-pipeline guard (a truly empty pipeline string is an
    // error regardless) but before computing boundary i/o units (those must
    // reflect the post-optimization boundary steps).
    let mut steps = crate::pipeline_opt::optimize_pipeline(steps);

    if steps.is_empty() {
        // All steps were optimized away (e.g. every step was a complete no-op).
        return Err(ProjError::MissingArg);
    }

    // Fill in `Whatever` boundary units of unit-agnostic steps from their
    // neighbours (PROJ `src/pipeline.cpp`), so a pipeline whose boundary step is
    // unit-agnostic (axisswap / push / pop / affine) still reports the correct
    // angular boundary units. Must run before the `left`/`right` computation
    // below so front/back units reflect the propagated values.
    propagate_whatever_units(&mut steps);

    let top_inverted = global_pl
        .entries
        .iter()
        .any(|(k, v)| k == "inv" && v.is_none());

    // --- Malformed-pipeline unit-consistency validation ---
    //
    // PROJ `pj_create_pipeline` (`src/pipeline.cpp`) checks, AFTER the
    // Whatever-propagation, that each step's forward output units match the next
    // step's forward input units, failing with the malformed-pipeline error
    // class if they don't. A `Whatever` boundary on either side is compatible
    // (skipped). Using the EFFECTIVE (inversion-aware, Classic->Projected
    // normalized) units matches PROJ's `pj_right(step[i])` / `pj_left(step[i+1])`.
    // This rejects e.g. `merc` straight into `merc` (projected metres into
    // angular radians) or `merc` into `helmert` (projected metres into
    // cartesian).
    for pair in steps.windows(2) {
        let curr_right = effective_right(&pair[0]);
        let next_left = effective_left(&pair[1]);
        if curr_right == IoUnits::Whatever || next_left == IoUnits::Whatever {
            continue;
        }
        if curr_right != next_left {
            return Err(ProjError::InvalidOp);
        }
    }

    // --- Forward-path availability validation (reject un-runnable steps) ---
    //
    // PROJ `pj_create_pipeline` requires a valid FORWARD path: every step that
    // is not omitted in the forward direction must provide the operation for
    // the direction it will actually run. A step whose net direction is inverse
    // (its own `+inv`, XOR the pipeline's global `+inv`) but whose operation has
    // no inverse (e.g. forward-only `urm5`, or a singular `affine`) is rejected
    // at creation with the no-inverse error — matching PROJ erroring rather than
    // producing NaN at transform time. The `omit` role swaps under a global
    // `+inv` (the forward pass then runs the inverse traversal), so a step
    // omitted in the forward direction is the `omit_inv` one when `top_inverted`.
    for step in &steps {
        let omitted_in_fwd = if top_inverted {
            step.omit_inv
        } else {
            step.omit_fwd
        };
        if omitted_in_fwd {
            continue;
        }
        let net_inverted = top_inverted ^ step.inverted;
        if net_inverted && !step.inner_has_inverse() {
            return Err(ProjError::NoInverseOp);
        }
    }

    // Pipeline i/o units mirror PROJ src/pipeline.cpp
    // (`P->left = pj_left(front); P->right = pj_right(back)`): take the
    // EFFECTIVE units of the boundary steps so an inverted boundary step
    // (e.g. `+proj=cart +inv`, whose effective output is geographic radians)
    // propagates correctly to the pipeline's forward-sense units.
    //
    // When the whole pipeline carries a global `+inv` (`top_inverted`), its
    // forward pass runs the steps' inverse traversal (last step first, each
    // inverted). The forward-sense boundary units are therefore SWAPPED
    // relative to the non-inverted pipeline: the forward input is what the last
    // step's forward would output (`pj_right(back)`), and the forward output is
    // what the first step's forward would input (`pj_left(front)`). Without this
    // swap a single-step `+proj=pipeline +inv +step +proj=urm5 +inv` reports its
    // boundary units backwards, so the degree<->radian conversion at the
    // boundary is skipped and the (radian-domain) step is fed raw degrees.
    let (left, right) = if top_inverted {
        (
            steps
                .last()
                .map(effective_right)
                .unwrap_or(IoUnits::Whatever),
            steps
                .first()
                .map(effective_left)
                .unwrap_or(IoUnits::Whatever),
        )
    } else {
        (
            steps
                .first()
                .map(effective_left)
                .unwrap_or(IoUnits::Whatever),
            steps
                .last()
                .map(effective_right)
                .unwrap_or(IoUnits::Whatever),
        )
    };

    Ok(Pj {
        operation: Box::new(Pipeline { steps }),
        ellipsoid: global_ellipsoid,
        // A pipeline wrapper does not compute map-scale factors on its own
        // (its boundary units are not a single projection's lon/lat), so the
        // factors metric simply carries the global ellipsoid's `es`.
        factors_es: global_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,
        // The pipeline wrapper itself never carries `+lon_wrap`: a global
        // `+lon_wrap` is inherited by every step through `appended_globals`
        // (like `+lon_0`) and applied by the step whose own forward pass
        // produces the pipeline's angular output, mirroring PROJ's per-step
        // `pj_init_ctx` inheritance (`src/pipeline.cpp`). Setting it here too
        // would double-wrap (a no-op since `adjlon` is idempotent, but not
        // the single source of truth the field's contract implies).
        lon_wrap_center: None,
        is_latlong: false,
        left,
        right,
        inverted: top_inverted,
        bypass_prepare_finalize: true,
        omit_fwd: false,
        omit_inv: false,
        ad_proj: None,
        op_name: String::new(),
        // The pipeline wrapper itself is not an `axisswap` step.
        axisswap_order: None,
    })
}

/// Transform a single coordinate in the given direction.
///
/// The `inverted` flag is honored inside [`Pj::forward`]/[`Pj::inverse`], so the
/// direction is not swapped here.
pub fn trans(pj: &Pj, dir: Direction, c: Coord) -> ProjResult<Coord> {
    match dir {
        Direction::Fwd => pj.forward(c),
        Direction::Inv => pj.inverse(c),
        Direction::Ident => Ok(c),
    }
}

/// Transform a slice of coordinates in place.
///
/// A failed coordinate becomes [`Coord::error`]; the loop does not abort early.
pub fn trans_array(pj: &Pj, dir: Direction, coords: &mut [Coord]) -> ProjResult<()> {
    for c in coords.iter_mut() {
        match trans(pj, dir, *c) {
            Ok(r) => *c = r,
            Err(_) => *c = Coord::error(),
        }
    }
    Ok(())
}

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

    // audit confirmed[11]: end-to-end coverage (through the public `create`
    // entry point, not just `setup_ellipsoid` in isolation) that a degenerate
    // ellipsoid is rejected for both a plain single operation and a pipeline
    // step. Cross-checked live against Homebrew PROJ 9.7.0: `projinfo
    // "+proj=longlat +a=0 +b=0 +no_defs +type=crs"` fails with error 1027
    // (`pj_init_ctx: Must specify ellipsoid or sphere`).
    #[test]
    fn create_rejects_zero_semi_major_axis() {
        assert!(create("+proj=longlat +a=0 +b=0 +ellps=WGS84").is_err());
        // `+a=0` must win over an explicit `+ellps` baseline (setup.rs step
        // 3a: "+a overrides the size"), so the trailing `+ellps=WGS84` above
        // must not mask the rejection; assert the same without it too.
        assert!(create("+proj=longlat +a=0 +b=0").is_err());
    }

    #[test]
    fn create_rejects_negative_semi_major_axis() {
        assert!(create("+proj=longlat +a=-6378137 +b=-6356752").is_err());
    }

    #[test]
    fn create_rejects_negative_flattening() {
        assert!(create("+proj=merc +a=6378137 +rf=-298.257223563").is_err());
    }

    #[test]
    fn create_pipeline_rejects_degenerate_step_ellipsoid() {
        // A single bad step must fail the whole pipeline build, matching
        // PROJ's `pj_create_pipeline` (a step's `pj_create_argv` failure
        // aborts the whole pipeline construction).
        assert!(create(
            "+proj=pipeline +step +proj=merc +ellps=WGS84 +step +proj=longlat +a=0 +b=0"
        )
        .is_err());
    }

    #[test]
    fn merc_forward_and_round_trip() {
        let pj = create("+proj=merc +ellps=WGS84").unwrap();
        let input = Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0);
        let fwd = trans(&pj, Direction::Fwd, input).unwrap();
        let f = fwd.v();
        assert!((f[0] - 1335833.8895192828).abs() < 1e-6, "x got {}", f[0]);
        assert!(
            (f[1] - 7_326_837.715_045_549).abs() < 1e-6,
            "y got {}",
            f[1]
        );
        let inv = trans(&pj, Direction::Inv, fwd).unwrap();
        let i = inv.v();
        assert!((i[0] - 12.0 * DEG_TO_RAD).abs() < 1e-9, "lam got {}", i[0]);
        assert!((i[1] - 55.0 * DEG_TO_RAD).abs() < 1e-9, "phi got {}", i[1]);
    }

    #[test]
    fn utm_known_values() {
        let pj = create("+proj=utm +zone=32 +ellps=WGS84").unwrap();
        let at_origin = trans(
            &pj,
            Direction::Fwd,
            Coord::new(9.0 * DEG_TO_RAD, 0.0, 0.0, 0.0),
        )
        .unwrap();
        let o = at_origin.v();
        assert!((o[0] - 500000.0).abs() < 1e-6, "x got {}", o[0]);
        assert!((o[1] - 0.0).abs() < 1e-6, "y got {}", o[1]);

        let p = trans(
            &pj,
            Direction::Fwd,
            Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let pv = p.v();
        assert!(
            (pv[0] - 691_875.632_137_542).abs() < 1e-6,
            "x got {}",
            pv[0]
        );
        assert!(
            (pv[1] - 6_098_907.825_129_169).abs() < 1e-6,
            "y got {}",
            pv[1]
        );
    }

    #[test]
    fn etmerc_central_meridian() {
        let pj = create("+proj=etmerc +lon_0=9 +ellps=WGS84").unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(9.0 * DEG_TO_RAD, 50.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let o = out.v();
        assert!(o[0].abs() < 1e-6, "x got {}", o[0]);
        assert!(
            (o[1] - 5_540_847.041_684_148).abs() < 1e-6,
            "y got {}",
            o[1]
        );
    }

    #[test]
    fn pipeline_utm_round_trip() {
        let pj = create(
            "+proj=pipeline +step +proj=utm +zone=32 +ellps=WGS84 +step +proj=utm +zone=32 +ellps=WGS84 +inv",
        )
        .unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let o = out.v();
        assert!((o[0] - 12.0 * DEG_TO_RAD).abs() < 1e-9, "lam got {}", o[0]);
        assert!((o[1] - 55.0 * DEG_TO_RAD).abs() < 1e-9, "phi got {}", o[1]);
    }

    #[test]
    fn trans_array_maps_in_place() {
        let pj = create("+proj=merc +ellps=WGS84").unwrap();
        let mut coords = [
            Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0),
            Coord::new(0.0, 0.0, 0.0, 0.0),
        ];
        trans_array(&pj, Direction::Fwd, &mut coords).unwrap();
        assert!((coords[0].v()[0] - 1335833.8895192828).abs() < 1e-6);
        assert!(coords[1].v()[0].abs() < 1e-6);
    }

    /// Build a minimal NTv2 file in memory with a uniform shift for test use.
    ///
    /// Covers lon [-1°, 0°] and lat [0°, 1°] — a 2×2 grid whose every node
    /// carries the same `lat_sec`/`lon_sec` shift (arcseconds).
    fn make_uniform_shift_ntv2(lat_sec: f32, lon_sec: f32) -> Vec<u8> {
        let mut buf = Vec::new();

        // File overview header (11 records × 16 bytes = 176 bytes)
        // Record 0: NUM_OREC
        buf.extend_from_slice(b"NUM_OREC");
        buf.extend_from_slice(&11i32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        // Record 1: NUM_SREC
        buf.extend_from_slice(b"NUM_SREC");
        buf.extend_from_slice(&11i32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        // Record 2: NUM_FILE (1 subfile)
        buf.extend_from_slice(b"NUM_FILE");
        buf.extend_from_slice(&1u32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        // Record 3: GS_TYPE = "SECONDS "
        buf.extend_from_slice(b"GS_TYPE ");
        buf.extend_from_slice(b"SECONDS ");
        // Records 4-10: unused filler (7 × 16 = 112 bytes)
        buf.extend_from_slice(&[0u8; 112]);

        // Sub-grid header (11 records × 16 bytes = 176 bytes)
        // Record 0: SUB_NAME
        buf.extend_from_slice(b"SUB_NAME");
        buf.extend_from_slice(b"TESTGRID");
        // Record 1: PARENT
        buf.extend_from_slice(b"PARENT  ");
        buf.extend_from_slice(b"NONE    ");
        // Record 2: CREATED
        buf.extend_from_slice(b"CREATED ");
        buf.extend_from_slice(b"20240101");
        // Record 3: UPDATED
        buf.extend_from_slice(b"UPDATED ");
        buf.extend_from_slice(b"20240101");
        // Record 4: S_LAT = 0 arcsec (equator)
        buf.extend_from_slice(b"S_LAT   ");
        buf.extend_from_slice(&0.0f64.to_le_bytes());
        // Record 5: N_LAT = 3600 arcsec = 1 degree
        buf.extend_from_slice(b"N_LAT   ");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        // Record 6: E_LONG = 0 arcsec (east boundary, positive-westward = 0°W = 0°E)
        buf.extend_from_slice(b"E_LONG  ");
        buf.extend_from_slice(&0.0f64.to_le_bytes());
        // Record 7: W_LONG = 3600 arcsec (west boundary, positive-westward = 1°W = -1°E)
        buf.extend_from_slice(b"W_LONG  ");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        // Record 8: LAT_INC = 3600 arcsec (1 degree → 2 rows)
        buf.extend_from_slice(b"LAT_INC ");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        // Record 9: LONG_INC = 3600 arcsec (1 degree → 2 cols)
        buf.extend_from_slice(b"LONG_INC");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        // Record 10: GS_COUNT = 4 (2×2 grid)
        buf.extend_from_slice(b"GS_COUNT");
        buf.extend_from_slice(&4i32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);

        // Data: 4 cells × 4 f32 = 64 bytes (uniform shift everywhere)
        for _ in 0..4 {
            buf.extend_from_slice(&lat_sec.to_le_bytes()); // lat shift
            buf.extend_from_slice(&lon_sec.to_le_bytes()); // lon shift
            buf.extend_from_slice(&0.0f32.to_le_bytes()); // lat accuracy
            buf.extend_from_slice(&0.0f32.to_le_bytes()); // lon accuracy
        }

        buf
    }

    /// Zero-shift NTv2 grid: `hgridshift` is an identity over the covered region.
    fn make_zero_shift_ntv2() -> Vec<u8> {
        make_uniform_shift_ntv2(0.0, 0.0)
    }

    #[test]
    fn nadgrids_pipeline_constructed() {
        // Verify that +nadgrids= on a single proj-string builds a pipeline
        // containing an hgridshift step, and that the grid is used for the shift.
        // We use a zero-shift grid so the round-trip is identity.
        let mut ctx = crate::context::Context::new();
        ctx.register_grid("test.gsb", make_zero_shift_ntv2());

        // merc with nadgrids — the point (lon=-0.5°, lat=0.5°) is inside the grid
        let pj = create_with_ctx("+proj=merc +ellps=WGS84 +nadgrids=test.gsb", &ctx).unwrap();

        // With zero shift, result should match plain merc
        let plain = create("+proj=merc +ellps=WGS84").unwrap();
        let coord = Coord::new(-0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);
        let shifted = trans(&pj, Direction::Fwd, coord).unwrap();
        let expected = trans(&plain, Direction::Fwd, coord).unwrap();
        let sv = shifted.v();
        let ev = expected.v();
        assert!(
            (sv[0] - ev[0]).abs() < 1e-3,
            "x: got {}, expected {}",
            sv[0],
            ev[0]
        );
        assert!(
            (sv[1] - ev[1]).abs() < 1e-3,
            "y: got {}, expected {}",
            sv[1],
            ev[1]
        );
    }

    #[test]
    fn pipeline_omit_fwd_skips_middle_step() {
        // Pipeline: noop -> axisswap(+omit_fwd) -> noop
        // Forward: middle step omitted => x,y unchanged
        // Inverse: middle step applied => x,y swapped
        let pj = create(
            "+proj=pipeline \
             +step +proj=noop \
             +step +proj=axisswap +order=2,1 +omit_fwd \
             +step +proj=noop",
        )
        .unwrap();

        let input = Coord::new(1.0, 2.0, 3.0, 0.0);

        // Forward: middle step omitted, so output equals input (noop+noop)
        let fwd = trans(&pj, Direction::Fwd, input).unwrap();
        let fv = fwd.v();
        assert!(
            (fv[0] - 1.0).abs() < 1e-12,
            "forward x should be 1.0, got {}",
            fv[0]
        );
        assert!(
            (fv[1] - 2.0).abs() < 1e-12,
            "forward y should be 2.0, got {}",
            fv[1]
        );

        // Inverse: middle step included (axisswap), so x and y are swapped
        let inv = trans(&pj, Direction::Inv, input).unwrap();
        let iv = inv.v();
        assert!(
            (iv[0] - 2.0).abs() < 1e-12,
            "inverse x should be 2.0 (swapped), got {}",
            iv[0]
        );
        assert!(
            (iv[1] - 1.0).abs() < 1e-12,
            "inverse y should be 1.0 (swapped), got {}",
            iv[1]
        );
    }

    #[test]
    fn optimizer_removes_double_omit_step() {
        // A step with both +omit_fwd and +omit_inv is never executed — the
        // pipeline optimizer removes it. The remaining noop steps let coordinates
        // pass through unchanged.
        let pj = create(
            "+proj=pipeline \
             +step +proj=noop \
             +step +proj=noop +omit_fwd +omit_inv \
             +step +proj=noop",
        )
        .expect("pipeline with double-omit middle step");
        let input = Coord::new(1.0, 2.0, 3.0, 0.0);
        let fwd = trans(&pj, Direction::Fwd, input).expect("forward through optimized pipeline");
        let fv = fwd.v();
        assert!(
            (fv[0] - 1.0).abs() < 1e-12,
            "x should pass through unchanged, got {}",
            fv[0]
        );
        assert!(
            (fv[1] - 2.0).abs() < 1e-12,
            "y should pass through unchanged, got {}",
            fv[1]
        );
    }

    #[test]
    fn optimizer_cancels_axisswap_pair() {
        // Two consecutive axisswap(2,1) steps cancel each other (axisswap is
        // self-inverse). After optimization the pipeline reduces to noop only,
        // and the input coordinate passes through unchanged.
        let pj = create(
            "+proj=pipeline \
             +step +proj=noop \
             +step +proj=axisswap +order=2,1 \
             +step +proj=axisswap +order=2,1",
        )
        .expect("pipeline with noop + cancelling axisswap pair");
        let input = Coord::new(3.0, 7.0, 0.0, 0.0);
        let fwd = trans(&pj, Direction::Fwd, input).expect("forward through optimized pipeline");
        let fv = fwd.v();
        assert!(
            (fv[0] - 3.0).abs() < 1e-12,
            "x should be unchanged after axisswap cancellation, got {}",
            fv[0]
        );
        assert!(
            (fv[1] - 7.0).abs() < 1e-12,
            "y should be unchanged after axisswap cancellation, got {}",
            fv[1]
        );
    }

    #[test]
    fn optimizer_keeps_non_cancelling_three_cycle_pair() {
        // Regression for the axisswap-cancellation bug: `+order=2,3,1` twice is
        // a 3-cycle squared (net permutation 3,1,2), NOT the identity. The
        // optimizer must NOT remove the pair. Verified against PROJ `cct`:
        //   echo "1 2 3 0" | cct +proj=pipeline \
        //     +step +proj=axisswap +order=2,3,1 \
        //     +step +proj=axisswap +order=2,3,1   ->  3 1 2 0
        let pj = create(
            "+proj=pipeline \
             +step +proj=axisswap +order=2,3,1 \
             +step +proj=axisswap +order=2,3,1",
        )
        .expect("two 3-cycle axisswap steps");
        let out = trans(&pj, Direction::Fwd, Coord::new(1.0, 2.0, 3.0, 4.0))
            .unwrap()
            .v();
        assert_eq!(
            out,
            [3.0, 1.0, 2.0, 4.0],
            "3-cycle applied twice must be the net permutation 3,1,2, not identity"
        );
    }

    #[test]
    fn optimizer_keeps_swap_then_negate_pair() {
        // `+order=2,1` then `+order=1,-2` is not the identity — keep both.
        // Verified against PROJ `cct`:
        //   echo "1 2 3 0" | cct +proj=pipeline \
        //     +step +proj=axisswap +order=2,1 \
        //     +step +proj=axisswap +order=1,-2   ->  2 -1 3 0
        let pj = create(
            "+proj=pipeline \
             +step +proj=axisswap +order=2,1 \
             +step +proj=axisswap +order=1,-2",
        )
        .expect("swap then negate axisswap steps");
        let out = trans(&pj, Direction::Fwd, Coord::new(1.0, 2.0, 3.0, 4.0))
            .unwrap()
            .v();
        assert_eq!(
            out,
            [2.0, -1.0, 3.0, 4.0],
            "swap-then-negate must not be optimized to identity"
        );
    }

    #[test]
    fn optimizer_keeps_rotation_twice_pair() {
        // `+order=-2,1` twice is a 90° rotation squared (= 180° rotation),
        // not the identity. Verified against PROJ `cct`:
        //   echo "1 2 3 0" | cct +proj=pipeline \
        //     +step +proj=axisswap +order=-2,1 \
        //     +step +proj=axisswap +order=-2,1   ->  -1 -2 3 0
        let pj = create(
            "+proj=pipeline \
             +step +proj=axisswap +order=-2,1 \
             +step +proj=axisswap +order=-2,1",
        )
        .expect("two rotation axisswap steps");
        let out = trans(&pj, Direction::Fwd, Coord::new(1.0, 2.0, 3.0, 4.0))
            .unwrap()
            .v();
        assert_eq!(
            out,
            [-1.0, -2.0, 3.0, 4.0],
            "90° rotation applied twice must be a 180° rotation, not identity"
        );
    }

    #[test]
    fn optimizer_cancels_true_inverse_axisswap_pair() {
        // A permutation followed by its own inverse (`+order=2,3,1` then the
        // same `+inv`) is the identity and IS safely cancelled. A leading `noop`
        // keeps the pipeline non-empty after the pair is removed; the output
        // must equal the input.
        let pj = create(
            "+proj=pipeline \
             +step +proj=noop \
             +step +proj=axisswap +order=2,3,1 \
             +step +proj=axisswap +order=2,3,1 +inv",
        )
        .expect("noop + axisswap and its inverse");
        let out = trans(&pj, Direction::Fwd, Coord::new(1.0, 2.0, 3.0, 4.0))
            .unwrap()
            .v();
        assert_eq!(
            out,
            [1.0, 2.0, 3.0, 4.0],
            "a permutation and its inverse compose to the identity"
        );
    }

    #[test]
    fn optimizer_all_cancelled_returns_error() {
        // Every step has both +omit_fwd and +omit_inv — all steps are removed by
        // the optimizer, leaving zero steps. This must return an error.
        let result = create(
            "+proj=pipeline \
             +step +proj=noop +omit_fwd +omit_inv \
             +step +proj=noop +omit_fwd +omit_inv",
        );
        assert!(
            result.is_err(),
            "all-cancelled pipeline must return an error"
        );
    }

    // ---------------------------------------------------------------------
    // E1 findings: single-string datum shift / nadgrids / axis direction,
    // pipeline global inheritance, towgs84 padding. Expected values are
    // derived from Homebrew PROJ 9.x `cs2cs`/`proj` (see test comments).
    // ---------------------------------------------------------------------

    fn assert_xy(got: Coord, ex: [f64; 2], tol: f64, label: &str) {
        let g = got.v();
        assert!(
            (g[0] - ex[0]).abs() < tol,
            "{label} x: got {}, want {}",
            g[0],
            ex[0]
        );
        assert!(
            (g[1] - ex[1]).abs() < tol,
            "{label} y: got {}, want {}",
            g[1],
            ex[1]
        );
    }

    #[test]
    fn towgs84_3param_matches_cs2cs() {
        // cs2cs +proj=longlat +datum=WGS84 +to +proj=utm +zone=32 +ellps=intl
        //       +towgs84=-87,-98,-121   at 9 0  ->  500083.152337543 120.954458759
        let pj = create("+proj=utm +zone=32 +ellps=intl +towgs84=-87,-98,-121").unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(9.0 * DEG_TO_RAD, 0.0, 0.0, 0.0),
        )
        .unwrap();
        assert_xy(
            out,
            [500083.152337543, 120.954458759],
            1e-4,
            "utm intl towgs84",
        );
    }

    #[test]
    fn towgs84_7param_matches_cs2cs() {
        // merc intl +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993 at 174 -41
        // -> 19370336.340988029 -4984632.904443800
        let pj =
            create("+proj=merc +ellps=intl +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993")
                .unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(174.0 * DEG_TO_RAD, -41.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        assert_xy(
            out,
            [19370336.34098803, -4984632.9044438],
            1e-2,
            "merc intl 7-param towgs84",
        );
    }

    #[test]
    fn towgs84_datum_resolves_ellipsoid() {
        // +datum=carthage must resolve to ellps=clrk80ign for BOTH the cart
        // round-trip and the projection step, matching the explicit form.
        let via_datum = create("+proj=merc +datum=carthage").unwrap();
        let via_explicit = create("+proj=merc +ellps=clrk80ign +towgs84=-263.0,6.0,431.0").unwrap();
        let p = Coord::new(10.0 * DEG_TO_RAD, 34.0 * DEG_TO_RAD, 0.0, 0.0);
        let a = trans(&via_datum, Direction::Fwd, p).unwrap();
        let b = trans(&via_explicit, Direction::Fwd, p).unwrap();
        // cs2cs ground truth for the explicit form.
        assert_xy(
            a,
            [1113152.342921798, 4004375.48761513],
            1e-2,
            "merc datum=carthage",
        );
        let (av, bv) = (a.v(), b.v());
        assert!(
            (av[0] - bv[0]).abs() < 1e-6 && (av[1] - bv[1]).abs() < 1e-6,
            "datum vs explicit ellps mismatch: {av:?} vs {bv:?}"
        );
    }

    #[test]
    fn towgs84_4param_zero_padded_not_dropped() {
        // finding 4: a 4-6 element list keeps its rotation term (padded to 7),
        // never silently truncated to a 3-parameter shift.
        assert_eq!(normalize_towgs84("1,2,3,4"), "1,2,3,4,0,0,0");
        assert_eq!(normalize_towgs84("1,2"), "1,2,0");
        assert_eq!(normalize_towgs84("1,2,3,4,5,6,7,8"), "1,2,3,4,5,6,7");
        // End-to-end: +towgs84=1,2,3,4 equals the explicit +towgs84=1,2,3,4,0,0,0.
        let four = create("+proj=merc +ellps=intl +towgs84=1,2,3,4").unwrap();
        let out = trans(
            &four,
            Direction::Fwd,
            Coord::new(10.0 * DEG_TO_RAD, 40.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        // cs2cs ground truth for +towgs84=1,2,3,4,0,0,0.
        assert_xy(
            out,
            [1113337.903887998, 4838632.996508759],
            1e-2,
            "merc intl 4-param towgs84",
        );
    }

    #[test]
    fn nadgrids_applies_hgridshift_inverse() {
        // finding 2: on the forward path PROJ applies hgridshift INVERSE.
        let mut ctx = crate::context::Context::new();
        ctx.register_grid("shift.gsb", make_uniform_shift_ntv2(5.0, -8.0));
        let p = Coord::new(-0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);

        let engine = create_with_ctx("+proj=merc +ellps=WGS84 +nadgrids=shift.gsb", &ctx).unwrap();
        let inv_ref = create_with_ctx(
            "+proj=pipeline +step +proj=hgridshift +grids=shift.gsb +inv \
             +step +proj=merc +ellps=WGS84",
            &ctx,
        )
        .unwrap();
        let fwd_ref = create_with_ctx(
            "+proj=pipeline +step +proj=hgridshift +grids=shift.gsb \
             +step +proj=merc +ellps=WGS84",
            &ctx,
        )
        .unwrap();

        let e = trans(&engine, Direction::Fwd, p).unwrap().v();
        let i = trans(&inv_ref, Direction::Fwd, p).unwrap().v();
        let f = trans(&fwd_ref, Direction::Fwd, p).unwrap().v();

        assert!(
            (e[0] - i[0]).abs() < 1e-9 && (e[1] - i[1]).abs() < 1e-9,
            "engine must apply hgridshift +inv: engine={e:?} inv_ref={i:?}"
        );
        // With a nonzero shift the forward direction is measurably different,
        // proving the injected step really is inverse.
        assert!(
            (e[0] - f[0]).abs() > 1e-3 || (e[1] - f[1]).abs() > 1e-3,
            "engine must NOT match forward hgridshift: engine={e:?} fwd_ref={f:?}"
        );
    }

    #[test]
    fn nadgrids_projection_uses_source_ellipsoid() {
        // The injected projection step keeps the source ellipsoid (intl),
        // not the pipeline GRS80/WGS84 default.
        let mut ctx = crate::context::Context::new();
        ctx.register_grid("shift.gsb", make_uniform_shift_ntv2(5.0, -8.0));
        let p = Coord::new(-0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);
        let engine = create_with_ctx("+proj=merc +ellps=intl +nadgrids=shift.gsb", &ctx).unwrap();
        let reference = create_with_ctx(
            "+proj=pipeline +step +proj=hgridshift +grids=shift.gsb +inv \
             +step +proj=merc +ellps=intl",
            &ctx,
        )
        .unwrap();
        let e = trans(&engine, Direction::Fwd, p).unwrap().v();
        let r = trans(&reference, Direction::Fwd, p).unwrap().v();
        assert!(
            (e[0] - r[0]).abs() < 1e-9 && (e[1] - r[1]).abs() < 1e-9,
            "nadgrids projection must use intl: engine={e:?} reference={r:?}"
        );
    }

    #[test]
    fn nadgrids_takes_precedence_over_towgs84() {
        // When BOTH +nadgrids and +towgs84 are present, PROJ uses the grid
        // shift and ignores the Helmert: cs2cs_emulation_setup builds
        // P->hgridshift first, then `p = P->hgridshift ? nullptr : towgs84`
        // (src/create.cpp). The engine must match the nadgrids-only result and
        // NOT the towgs84-only result.
        let mut ctx = crate::context::Context::new();
        ctx.register_grid("shift.gsb", make_uniform_shift_ntv2(5.0, -8.0));
        let p = Coord::new(-0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);

        let both = create_with_ctx(
            "+proj=merc +ellps=WGS84 +towgs84=100,200,300 +nadgrids=shift.gsb",
            &ctx,
        )
        .unwrap();
        let grid_only =
            create_with_ctx("+proj=merc +ellps=WGS84 +nadgrids=shift.gsb", &ctx).unwrap();
        let helmert_only =
            create_with_ctx("+proj=merc +ellps=WGS84 +towgs84=100,200,300", &ctx).unwrap();

        let b = trans(&both, Direction::Fwd, p).unwrap().v();
        let g = trans(&grid_only, Direction::Fwd, p).unwrap().v();
        let h = trans(&helmert_only, Direction::Fwd, p).unwrap().v();

        assert!(
            (b[0] - g[0]).abs() < 1e-9 && (b[1] - g[1]).abs() < 1e-9,
            "nadgrids must win over towgs84: both={b:?} grid_only={g:?}"
        );
        assert!(
            (b[0] - h[0]).abs() > 1e-3 || (b[1] - h[1]).abs() > 1e-3,
            "with both present the towgs84 path must NOT be used: both={b:?} helmert_only={h:?}"
        );
    }

    #[test]
    fn pipeline_global_param_inherited_by_steps() {
        // finding 3: a global +lon_0 is inherited by each step (PROJ pipeline.cpp).
        // Verified against PROJ cct: GLOBAL == STEP-LOCAL, both != NONE.
        let p = Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0);
        let global = create("+proj=pipeline +lon_0=10 +step +proj=merc +ellps=WGS84").unwrap();
        let local = create("+proj=pipeline +step +proj=merc +lon_0=10 +ellps=WGS84").unwrap();
        let none = create("+proj=pipeline +step +proj=merc +ellps=WGS84").unwrap();
        let g = trans(&global, Direction::Fwd, p).unwrap().v();
        let l = trans(&local, Direction::Fwd, p).unwrap().v();
        let n = trans(&none, Direction::Fwd, p).unwrap().v();
        assert!(
            (g[0] - l[0]).abs() < 1e-9 && (g[1] - l[1]).abs() < 1e-9,
            "global lon_0 must equal step-local lon_0: {g:?} vs {l:?}"
        );
        assert!(
            (g[0] - n[0]).abs() > 1.0,
            "global lon_0 must change the result vs none: {g:?} vs {n:?}"
        );
    }

    #[test]
    fn axis_wsu_negates_output() {
        // finding 6: +axis=wsu negates x and y (verified against cs2cs).
        let pj = create("+proj=merc +axis=wsu +ellps=WGS84").unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        assert_xy(
            out,
            [-1335833.8895192828, -7_326_837.715_045_549],
            1e-6,
            "merc axis=wsu",
        );
    }

    #[test]
    fn axis_neu_swaps_output() {
        // finding 6: +axis=neu swaps x and y.
        let pj = create("+proj=merc +axis=neu +ellps=WGS84").unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        assert_xy(
            out,
            [7_326_837.715_045_549, 1335833.8895192828],
            1e-6,
            "merc axis=neu",
        );
    }

    #[test]
    fn axis_round_trips() {
        let pj = create("+proj=merc +axis=wsu +ellps=WGS84").unwrap();
        let input = Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 0.0, 0.0);
        let fwd = trans(&pj, Direction::Fwd, input).unwrap();
        let inv = trans(&pj, Direction::Inv, fwd).unwrap().v();
        assert!(
            (inv[0] - 12.0 * DEG_TO_RAD).abs() < 1e-9 && (inv[1] - 55.0 * DEG_TO_RAD).abs() < 1e-9,
            "axis round trip: {inv:?}"
        );
    }

    #[test]
    fn lon_wrap_wraps_output_into_center_band() {
        // audit gap 3: `+lon_wrap=<center>` must wrap angular OUTPUT into
        // [center-180, center+180), matching PROJ `fwd_finalize`'s
        // `is_long_wrap_set` branch. Reference values cross-checked live
        // against Homebrew PROJ 9.7.0 `cct` on an explicit
        // deg->rad->longlat(+lon_wrap=180)->rad->deg pipeline:
        //   -170 10  -> 190 10  (wrapped: -170 + 360 = 190)
        //     10 10  ->  10 10  (already inside [0, 360))
        //    190 10  -> 190 10  (already inside [0, 360))
        let pj = create("+proj=longlat +lon_wrap=180 +ellps=WGS84").unwrap();

        let cases = [
            (-170.0, 10.0, 190.0, 10.0),
            (10.0, 10.0, 10.0, 10.0),
            (190.0, 10.0, 190.0, 10.0),
        ];
        for (lam_in, phi_in, lam_want, phi_want) in cases {
            let out = trans(
                &pj,
                Direction::Fwd,
                Coord::new(lam_in * DEG_TO_RAD, phi_in * DEG_TO_RAD, 0.0, 0.0),
            )
            .unwrap();
            let v = out.v();
            assert!(
                (v[0] - lam_want * DEG_TO_RAD).abs() < 1e-9,
                "lam_in={lam_in}: got {} deg, want {lam_want} deg",
                v[0] * oxiproj_core::RAD_TO_DEG
            );
            assert!(
                (v[1] - phi_want * DEG_TO_RAD).abs() < 1e-9,
                "lam_in={lam_in}: phi got {} deg, want {phi_want} deg",
                v[1] * oxiproj_core::RAD_TO_DEG
            );
        }
    }

    #[test]
    fn lon_wrap_oracle_sweep() {
        // audit gap 3: broader sweep, including boundary/wraparound inputs,
        // cross-checked live against Homebrew PROJ 9.7.0 `cct` on
        // `+proj=pipeline +lon_wrap=180 +step +proj=longlat +ellps=WGS84`
        // (`-z 0 -t 0` pinned; `+proj=longlat` rejects HUGE_VAL z/t inside a
        // pipeline step). -180 -> 180 and 360 -> 0 exercise the raw-input
        // `adjlon` normalization in `fwd_prepare` running BEFORE
        // `fwd_finalize`'s `+lon_wrap` wrap ever sees the value (both are
        // pre-existing `forward_impl` prepare-stage behavior, untouched by
        // this fix) composing with the new finalize-stage wrap.
        let pj = create("+proj=pipeline +lon_wrap=180 +step +proj=longlat +ellps=WGS84").unwrap();
        let cases: [(f64, f64); 11] = [
            (-170.0, 190.0),
            (-180.0, 180.0),
            (0.0, 0.0),
            (10.0, 10.0),
            (170.0, 170.0),
            (179.9999, 179.9999),
            (180.0, 180.0),
            (180.0001, 180.0001),
            (190.0, 190.0),
            (359.9999, 359.9999),
            (360.0, 0.0),
        ];
        for (lam_in, lam_want) in cases {
            let out = trans(
                &pj,
                Direction::Fwd,
                Coord::new(lam_in * DEG_TO_RAD, 5.0 * DEG_TO_RAD, 0.0, 0.0),
            )
            .unwrap();
            let got_deg = out.v()[0] * oxiproj_core::RAD_TO_DEG;
            assert!(
                (got_deg - lam_want).abs() < 1e-6,
                "lam_in={lam_in}: got {got_deg} deg, want {lam_want} deg"
            );
        }
    }

    #[test]
    fn lon_wrap_absent_leaves_longitude_unwrapped() {
        // Baseline: without `+lon_wrap`, -170 stays -170 (already inside the
        // default (-180, 180] band adjlon normalizes to), unlike the
        // +lon_wrap=180 case above which relocates it to 190.
        let pj = create("+proj=longlat +ellps=WGS84").unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(-170.0 * DEG_TO_RAD, 10.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let v = out.v();
        assert!(
            (v[0] - (-170.0 * DEG_TO_RAD)).abs() < 1e-9,
            "lam got {}",
            v[0]
        );
    }

    #[test]
    fn lon_wrap_precedes_axisswap_in_synthesized_pipeline() {
        // audit gap 3: when `+axis` forces the `+lon_wrap`/axisswap
        // synthesized-pipeline path (`create_single` / `strip_from_projection_step`),
        // the wrap must run BEFORE the axisswap step, matching PROJ
        // `fwd_finalize`'s `is_long_wrap_set` wrap preceding `P->axisswap`
        // (`src/fwd.cpp`). Reference cross-checked live against Homebrew PROJ
        // 9.7.0 `cct` on the equivalent deg->rad->longlat(+lon_wrap=180
        // +axis=wsu)->rad->deg pipeline: -170 10 -> -190 -10 (wrap to 190,
        // then axisswap wsu negates both components).
        let pj = create("+proj=longlat +lon_wrap=180 +axis=wsu +ellps=WGS84").unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(-170.0 * DEG_TO_RAD, 10.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let v = out.v();
        assert!(
            (v[0] - (-190.0 * DEG_TO_RAD)).abs() < 1e-9,
            "x got {} deg",
            v[0] * oxiproj_core::RAD_TO_DEG
        );
        assert!(
            (v[1] - (-10.0 * DEG_TO_RAD)).abs() < 1e-9,
            "y got {} deg",
            v[1] * oxiproj_core::RAD_TO_DEG
        );
    }

    #[test]
    fn lon_wrap_pipeline_global_wraps_single_step_output() {
        // A global `+lon_wrap` on a `+proj=pipeline` is inherited by every
        // step (like `+lon_0`; see `create.rs`'s `appended_globals`), and
        // applied on each step's own true forward pass. Cross-checked live
        // against Homebrew PROJ 9.7.0 `cct` on the identical raw pipeline
        // string (-170 10 -> 190 10, z/t pinned since `+proj=longlat`
        // rejects HUGE_VAL z/t inside a pipeline step).
        let pj = create("+proj=pipeline +lon_wrap=180 +step +proj=longlat +ellps=WGS84").unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(-170.0 * DEG_TO_RAD, 10.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let v = out.v();
        assert!(
            (v[0] - (190.0 * DEG_TO_RAD)).abs() < 1e-9,
            "x got {} deg",
            v[0] * oxiproj_core::RAD_TO_DEG
        );
    }

    #[test]
    fn lon_wrap_pipeline_global_undone_by_inverted_round_trip_step() {
        // A two-step round-trip pipeline (plain forward, then the same
        // operation `+inv`) does NOT keep the first step's `+lon_wrap`
        // relocation: the second (inverted) step runs its own `inv_finalize`
        // equivalent, which in PROJ unconditionally re-applies the plain
        // `+over`-gated `adjlon` (`src/inv.cpp`) regardless of `+lon_wrap` —
        // PROJ's `inv_finalize` has no `is_long_wrap_set` branch at all, only
        // `fwd_finalize` does. So the relocation set up by step 0's forward
        // pass is normalized straight back to the standard band by step 1's
        // inverse pass. Cross-checked live against Homebrew PROJ 9.7.0 `cct`
        // on the identical raw pipeline string: -170 10 -> -170 10 (NOT 190).
        let pj = create(
            "+proj=pipeline +lon_wrap=180 +step +proj=longlat +ellps=WGS84 +step +proj=longlat +ellps=WGS84 +inv",
        )
        .unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(-170.0 * DEG_TO_RAD, 10.0 * DEG_TO_RAD, 0.0, 0.0),
        )
        .unwrap();
        let v = out.v();
        assert!(
            (v[0] - (-170.0 * DEG_TO_RAD)).abs() < 1e-9,
            "x got {} deg",
            v[0] * oxiproj_core::RAD_TO_DEG
        );
    }

    #[test]
    fn lon_wrap_rejects_excessive_center() {
        // PROJ `src/init.cpp`: `if (!(fabs(long_wrap_center) < 10 * M_TWOPI))`
        // rejects an excessive center with error 1027 ("Invalid value for
        // lon_wrap"). Cross-checked live: Homebrew PROJ 9.7.0 `projinfo
        // "+proj=longlat +lon_wrap=999999 +ellps=WGS84 +type=crs"` fails with
        // exactly this message.
        assert!(create("+proj=longlat +lon_wrap=999999 +ellps=WGS84").is_err());
    }

    #[test]
    fn lon_wrap_accepts_bare_flag_as_zero_center() {
        // A bare `+lon_wrap` (no `=value`) must not error; PROJ's numeric
        // param parse falls back to 0.0 for a missing operand.
        assert!(create("+proj=longlat +lon_wrap +ellps=WGS84").is_ok());
    }

    /// Build a 2×2 GTX grid over lat/lon [0°,1°] with a uniform vertical shift
    /// (metres). Big-endian, matching `oxiproj_grids::read_gtx`.
    fn make_uniform_gtx(shift_m: f32) -> Vec<u8> {
        let mut buf = Vec::new();
        buf.extend_from_slice(&0.0f64.to_be_bytes()); // south_lat (deg)
        buf.extend_from_slice(&0.0f64.to_be_bytes()); // west_lon (deg)
        buf.extend_from_slice(&1.0f64.to_be_bytes()); // lat_inc (deg)
        buf.extend_from_slice(&1.0f64.to_be_bytes()); // lon_inc (deg)
        buf.extend_from_slice(&2i32.to_be_bytes()); // rows
        buf.extend_from_slice(&2i32.to_be_bytes()); // cols
        for _ in 0..4 {
            buf.extend_from_slice(&shift_m.to_be_bytes());
        }
        buf
    }

    #[test]
    fn geoidgrids_injects_vgridshift() {
        // E2 finding 1: classic `+geoidgrids=` must no longer be silently
        // ignored. PROJ (cs2cs_emulation_setup) turns it into a vgridshift
        // helper applied FORWARD on the forward path (fwd_prepare), moving
        // between geometric and orthometric height. PROJ's vgridshift forward
        // SUBTRACTS the grid value (default forward_multiplier = -1.0), so a
        // uniform +10 m geoid grid REMOVES 10 m from the height before the
        // projection runs (verified: `cct +proj=vgridshift` yields z-10).
        let mut ctx = crate::context::Context::new();
        ctx.register_grid("geoid.gtx", make_uniform_gtx(10.0));

        // Point (lon=0.5°, lat=0.5°) is inside the grid; carry a nonzero height.
        let p = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 5.0, 0.0);

        let engine =
            create_with_ctx("+proj=merc +ellps=WGS84 +geoidgrids=geoid.gtx", &ctx).unwrap();
        // Equivalent explicit pipeline: vgridshift (FWD) then the projection.
        let fwd_ref = create_with_ctx(
            "+proj=pipeline +step +proj=vgridshift +grids=geoid.gtx \
             +step +proj=merc +ellps=WGS84",
            &ctx,
        )
        .unwrap();
        let plain = create("+proj=merc +ellps=WGS84").unwrap();

        let e = trans(&engine, Direction::Fwd, p).unwrap().v();
        let r = trans(&fwd_ref, Direction::Fwd, p).unwrap().v();
        let n = trans(&plain, Direction::Fwd, p).unwrap().v();

        // Correct wiring & direction: engine equals the explicit vgridshift+merc.
        assert!(
            (e[0] - r[0]).abs() < 1e-9 && (e[1] - r[1]).abs() < 1e-9 && (e[2] - r[2]).abs() < 1e-6,
            "geoidgrids must equal explicit vgridshift pipeline: engine={e:?} ref={r:?}"
        );
        // Not silently ignored: height differs from plain merc by the -10 m
        // shift (PROJ vgridshift forward subtracts the grid value).
        assert!(
            (e[2] - (n[2] - 10.0)).abs() < 1e-6,
            "geoidgrids must subtract the +10 m geoid shift from height: engine z={}, plain z={}",
            e[2],
            n[2]
        );
        // Horizontal output is unaffected by the vertical shift.
        assert!(
            (e[0] - n[0]).abs() < 1e-6 && (e[1] - n[1]).abs() < 1e-6,
            "geoidgrids must not move x/y: engine={e:?} plain={n:?}"
        );
    }

    #[test]
    fn geoidgrids_round_trips() {
        // Forward then inverse restores the original height: the injected
        // vgridshift runs INVERSE on the inverse path (inv_finalize).
        let mut ctx = crate::context::Context::new();
        ctx.register_grid("geoid.gtx", make_uniform_gtx(10.0));
        let engine =
            create_with_ctx("+proj=merc +ellps=WGS84 +geoidgrids=geoid.gtx", &ctx).unwrap();
        let p = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 5.0, 0.0);
        let fwd = trans(&engine, Direction::Fwd, p).unwrap();
        let inv = trans(&engine, Direction::Inv, fwd).unwrap().v();
        assert!(
            (inv[0] - 0.5 * DEG_TO_RAD).abs() < 1e-9
                && (inv[1] - 0.5 * DEG_TO_RAD).abs() < 1e-9
                && (inv[2] - 5.0).abs() < 1e-6,
            "geoidgrids round trip: {inv:?}"
        );
    }

    #[test]
    fn hgridshift_resolves_grid_from_proj_data_via_public_api() {
        // E4 finding: grid loading must work end-to-end from a local PROJ_DATA
        // directory through the public engine API (`create`), not only via the
        // in-memory registry. Write a synthetic NTv2 grid to a unique temp dir,
        // point PROJ_DATA at it, and run `+proj=hgridshift`. `cargo nextest`
        // runs every test in its own process, so mutating PROJ_DATA here cannot
        // race other tests; the prior value is still saved and restored.
        use std::io::Write;
        let ns = std::time::SystemTime::now()
            .duration_since(std::time::UNIX_EPOCH)
            .map(|d| d.as_nanos())
            .unwrap_or(0);
        let dir = std::env::temp_dir().join(format!("oxiproj_e4_projdata_{ns}"));
        std::fs::create_dir_all(&dir).expect("create temp PROJ_DATA dir");
        // +1 deg latitude shift, 0 deg longitude shift over lon[-1,0], lat[0,1].
        let grid_bytes = make_uniform_shift_ntv2(3600.0, 0.0);
        {
            let mut f = std::fs::File::create(dir.join("e4_shift.gsb")).expect("create grid file");
            f.write_all(&grid_bytes).expect("write grid bytes");
        }

        let saved = std::env::var_os("PROJ_DATA");
        std::env::set_var("PROJ_DATA", &dir);

        // Point inside the grid (lon=-0.5 deg, lat=0.5 deg).
        let p = Coord::new(-0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);
        let disk_out = create("+proj=hgridshift +grids=e4_shift.gsb")
            .ok()
            .and_then(|pj| trans(&pj, Direction::Fwd, p).ok().map(|c| c.v()));

        // Restore env and remove temp files before any assertion can unwind.
        match saved {
            Some(v) => std::env::set_var("PROJ_DATA", v),
            None => std::env::remove_var("PROJ_DATA"),
        }
        let _ = std::fs::remove_dir_all(&dir);

        let out = disk_out
            .expect("hgridshift must build and transform using a grid resolved from PROJ_DATA");
        // The +1 deg latitude shift was applied; longitude is unchanged.
        assert!(
            (out[0] - (-0.5 * DEG_TO_RAD)).abs() < 1e-9,
            "lon should be unchanged, got {}",
            out[0]
        );
        assert!(
            (out[1] - 1.5 * DEG_TO_RAD).abs() < 1e-9,
            "lat should be shifted +1 deg, got {}",
            out[1]
        );

        // Cross-check: the identical grid registered in memory yields the same
        // result, tying the disk-resolution path to the already-verified
        // in-memory registration path.
        let mut ctx = crate::context::Context::new();
        ctx.register_grid("e4_shift.gsb", grid_bytes);
        let mem = create_with_ctx("+proj=hgridshift +grids=e4_shift.gsb", &ctx)
            .expect("in-memory hgridshift");
        let mem_out = trans(&mem, Direction::Fwd, p)
            .expect("in-memory transform")
            .v();
        assert!(
            (out[0] - mem_out[0]).abs() < 1e-12 && (out[1] - mem_out[1]).abs() < 1e-12,
            "disk-resolved result must match in-memory-registered result: disk={out:?} mem={mem_out:?}"
        );
    }

    #[test]
    fn moll_factors_use_sphere_metric() {
        // Defect B: Mollweide is sphere-only (PROJ `moll.cpp` sets `P->es = 0`),
        // so its map-scale factors must be computed on the sphere even when an
        // ellipsoid is supplied. Ground truth from PROJ `proj -S +proj=moll
        // +ellps=WGS84` at (0,0):  <1.11072 0.900316 1 ...>  — areal scale s = 1
        // (Mollweide is equal-area), meridional h = 1.11072.
        let pj = create("+proj=moll +ellps=WGS84").unwrap();
        let f = pj.factors(Coord::new(0.0, 0.0, 0.0, 0.0)).unwrap();
        assert!(
            (f.s - 1.0).abs() < 1e-6,
            "moll areal scale must be 1.0 (sphere metric), got {}",
            f.s
        );
        assert!(
            (f.h - 1.110_720_734_5).abs() < 1e-4,
            "moll meridional scale must match PROJ 1.11072, got {}",
            f.h
        );
        // Exact-AD factors path must agree with the numeric one (both es=0).
        let fe = pj.factors_exact(Coord::new(0.0, 0.0, 0.0, 0.0)).unwrap();
        assert!(
            (fe.areal_scale - 1.0).abs() < 1e-6,
            "moll exact areal scale must be 1.0, got {}",
            fe.areal_scale
        );
    }

    #[test]
    fn sphere_only_projections_report_unit_areal_scale() {
        // sinu, hammer, eck4 are all sphere-only equal-area projections; given an
        // ellipsoid PROJ still reports areal scale 1.0 (verified vs `proj -S`).
        for name in ["sinu", "hammer", "eck4"] {
            let pj = create(&format!("+proj={name} +ellps=WGS84")).unwrap();
            let f = pj
                .factors(Coord::new(20.0 * DEG_TO_RAD, 40.0 * DEG_TO_RAD, 0.0, 0.0))
                .unwrap();
            assert!(
                (f.s - 1.0).abs() < 1e-6,
                "{name} areal scale must be 1.0 (sphere metric), got {}",
                f.s
            );
        }
    }

    #[test]
    fn ellipsoidal_projection_factors_keep_eccentricity() {
        // Genuinely ellipsoidal projections must NOT be forced onto the sphere
        // metric. Ellipsoidal Mercator at (20,40): PROJ `proj -S` reports
        // h = k = 1.30360 (ellipsoidal), distinct from the spherical
        // sec(40°) = 1.305407. Confirms `factors_es` stays = ellipsoid.es here.
        let pj = create("+proj=merc +ellps=WGS84").unwrap();
        let f = pj
            .factors(Coord::new(20.0 * DEG_TO_RAD, 40.0 * DEG_TO_RAD, 0.0, 0.0))
            .unwrap();
        assert!(
            (f.h - 1.303_600_689_3).abs() < 1e-4,
            "merc meridional scale must match ellipsoidal PROJ 1.30360, got {}",
            f.h
        );
        let spherical = 1.0 / 40.0_f64.to_radians().cos();
        assert!(
            (f.h - spherical).abs() > 1e-3,
            "merc factors must be ellipsoidal, not spherical {spherical}: got {}",
            f.h
        );
    }

    #[test]
    fn geocent_produces_ecef_via_engine() {
        // E2 finding 2: `+proj=geocent` must return real ECEF, not the input
        // radians. Ground truth from Homebrew PROJ 9.x:
        //   echo "12 55 100" | cct +proj=geocent +ellps=WGS84
        //   -> 3586525.7610  762339.5841  5201465.4384
        let pj = create("+proj=geocent +ellps=WGS84").unwrap();
        let out = trans(
            &pj,
            Direction::Fwd,
            Coord::new(12.0 * DEG_TO_RAD, 55.0 * DEG_TO_RAD, 100.0, 0.0),
        )
        .unwrap()
        .v();
        assert!((out[0] - 3586525.761017917).abs() < 1e-3, "X = {}", out[0]);
        assert!((out[1] - 762339.584102928).abs() < 1e-3, "Y = {}", out[1]);
        assert!((out[2] - 5201465.438406702).abs() < 1e-3, "Z = {}", out[2]);
        // Inverse restores the geodetic input.
        let back = trans(&pj, Direction::Inv, Coord::new(out[0], out[1], out[2], 0.0))
            .unwrap()
            .v();
        assert!(
            (back[0] - 12.0 * DEG_TO_RAD).abs() < 1e-9,
            "lam = {}",
            back[0]
        );
        assert!(
            (back[1] - 55.0 * DEG_TO_RAD).abs() < 1e-9,
            "phi = {}",
            back[1]
        );
        assert!((back[2] - 100.0).abs() < 1e-3, "h = {}", back[2]);
    }
}