oxiproj-engine 0.1.0

//! 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,
    }
}

/// 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)
    }
}

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
}

fn build_datum_shift_pipeline(
    params: &ParamList,
    towgs84_str: &str,
    ctx: &Context,
) -> ProjResult<Pj> {
    // Determine source ellipsoid name
    let src_ellps = params
        .get_str("ellps")
        .filter(|s| !s.is_empty())
        .map(|s| s.to_string())
        .or_else(|| {
            params
                .get_str("datum")
                .filter(|s| !s.is_empty())
                .and_then(|d| oxiproj_core::find_datum(d))
                .map(|d| d.ellipse_id.to_string())
        })
        .unwrap_or_else(|| "WGS84".to_string());

    // Parse towgs84 values to detect all-zero (identity)
    let vals: Vec<&str> = towgs84_str.split(',').collect();

    // All-zero means identity: skip datum shift injection
    let all_zero = vals
        .iter()
        .all(|v| matches!(v.trim(), "0" | "0.0" | "-0" | "-0.0"));
    if all_zero {
        return create_single_core(params, ctx);
    }

    // Build helmert params string
    let helmert_params = if vals.len() >= 7 {
        format!(
            "+proj=helmert +x={} +y={} +z={} +rx={} +ry={} +rz={} +s={} +convention=position_vector",
            vals[0].trim(),
            vals[1].trim(),
            vals[2].trim(),
            vals[3].trim(),
            vals[4].trim(),
            vals[5].trim(),
            vals[6].trim()
        )
    } else {
        format!(
            "+proj=helmert +x={} +y={} +z={}",
            vals.first().map(|v| v.trim()).unwrap_or("0"),
            vals.get(1).map(|v| v.trim()).unwrap_or("0"),
            vals.get(2).map(|v| v.trim()).unwrap_or("0")
        )
    };

    // Build the original projection step params (exclude towgs84, datum)
    let other_params: String = params
        .entries
        .iter()
        .filter(|(k, _)| k != "towgs84" && k != "datum")
        .map(|(k, v)| match v {
            Some(val) => format!("+{}={}", k, val),
            None => format!("+{}", k),
        })
        .collect::<Vec<_>>()
        .join(" ");

    // Synthetic pipeline:
    // step 1: cart (geodetic -> ECEF) using src ellipsoid
    // step 2: helmert (datum shift)
    // step 3: cart +inv (ECEF -> geodetic) using WGS84
    // step 4: the original projection
    let pipeline_str = format!(
        "+proj=pipeline \
         +step +proj=cart +ellps={src_ellps} \
         +step {helmert_params} \
         +step +proj=cart +inv +ellps=WGS84 \
         +step {other_params}"
    );

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

fn create_single(params: &ParamList, ctx: &Context) -> ProjResult<Pj> {
    // Inject datum-shift pipeline when towgs84 params are present
    if let Some(towgs84_str) = extract_towgs84_from_params(params) {
        return build_datum_shift_pipeline(params, &towgs84_str, ctx);
    }
    create_single_core(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.
    let mut global: Vec<(String, Option<String>)> = Vec::new();
    let mut steps_entries: Vec<Vec<(String, Option<String>)>> = Vec::new();
    let mut seen_step = false;
    for (k, v) in &params.entries {
        if k == "step" && v.is_none() {
            seen_step = true;
            steps_entries.push(Vec::new());
            continue;
        }
        if !seen_step {
            // Drop the pipeline marker from globals; it is never looked up.
            if k == "proj" && v.as_deref() == Some("pipeline") {
                continue;
            }
            global.push((k.clone(), v.clone()));
        } else if let Some(last) = steps_entries.last_mut() {
            last.push((k.clone(), v.clone()));
        }
    }

    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)?
    };

    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 step_ellipsoid = if has_ellipsoid_def(&step_pl) {
            crate::setup::setup_ellipsoid(&step_pl)?
        } else {
            global_ellipsoid
        };
        let mut step_pj = crate::registry::build_single_op(name, &step_pl, step_ellipsoid, ctx)?;
        step_pj.inverted = inverted;
        steps.push(step_pj);
    }

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

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

    // 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.
    let left = steps
        .first()
        .map(effective_left)
        .unwrap_or(IoUnits::Whatever);
    let right = steps
        .last()
        .map(effective_right)
        .unwrap_or(IoUnits::Whatever);

    Ok(Pj {
        operation: Box::new(Pipeline { steps }),
        ellipsoid: global_ellipsoid,
        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,
        is_latlong: false,
        left,
        right,
        inverted: top_inverted,
        bypass_prepare_finalize: true,
    })
}

/// 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;

    #[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);
    }
}