oxiproj-transformations 0.1.2

//! `unitconvert` conversion. Ported from src/conversions/unitconvert.cpp.
//!
//! The `unitconvert` pseudo-projection rescales coordinate components between
//! units. Horizontal components (slots 0 and 1) share a single linear/angular
//! factor derived from `+xy_in`/`+xy_out`, the vertical component (slot 2) uses
//! `+z_in`/`+z_out`, and the temporal component (slot 3) is converted through a
//! set of calendar-aware time unit conversions selected by `+t_in`/`+t_out`.

use crate::{TransBuild, TransParams};
use oxiproj_core::{
    find_angular_unit, find_linear_unit, Coord, IoUnits, Lp, Lpz, Operation, ProjError, ProjResult,
    Xy, Xyz,
};

/// Supported temporal units. Ported from src/conversions/unitconvert.cpp.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum TimeUnit {
    Mjd,
    DecimalYear,
    GpsWeek,
    YyyyMmDd,
}

/// Parse a temporal unit identifier into a [`TimeUnit`].
///
/// Ported from src/conversions/unitconvert.cpp.
fn parse_time_unit(id: &str) -> ProjResult<TimeUnit> {
    match id {
        "mjd" => Ok(TimeUnit::Mjd),
        "decimalyear" => Ok(TimeUnit::DecimalYear),
        "gps_week" => Ok(TimeUnit::GpsWeek),
        "yyyymmdd" => Ok(TimeUnit::YyyyMmDd),
        _ => Err(ProjError::IllegalArgValue),
    }
}

/// Return `true` for proleptic Gregorian leap years.
///
/// Ported from src/conversions/unitconvert.cpp.
fn is_leap_year(year: i64) -> bool {
    (year % 4 == 0 && year % 100 != 0) || year % 400 == 0
}

/// Number of days in `year`.
///
/// Ported from src/conversions/unitconvert.cpp.
fn days_in_year(year: i64) -> i64 {
    if is_leap_year(year) {
        366
    } else {
        365
    }
}

/// Number of days in `month` of `year` (1-based month).
///
/// Ported from src/conversions/unitconvert.cpp.
fn days_in_month(year: i64, month: i64) -> i64 {
    const DAYS: [i64; 12] = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31];
    let m = month.clamp(1, 12);
    let idx = (m - 1) as usize;
    let mut d = DAYS[idx];
    if m == 2 && is_leap_year(year) {
        d += 1;
    }
    d
}

/// Day-of-year number for the given calendar date.
///
/// Ported from src/conversions/unitconvert.cpp.
fn daynumber_in_year(year: i64, month: i64, day: i64) -> i64 {
    let m = month.clamp(1, 12);
    let day = day.clamp(1, days_in_month(year, m));
    let mut daynum = day;
    let mut mm = 1;
    while mm < m {
        daynum += days_in_month(year, mm);
        mm += 1;
    }
    daynum
}

/// Convert a decimal year to Modified Julian Date.
///
/// Ported from src/conversions/unitconvert.cpp.
fn decimalyear_to_mjd(dy: f64) -> f64 {
    if !((-10000.0..=10000.0).contains(&dy)) {
        return 0.0;
    }
    let year = dy.floor() as i64;
    let frac = dy - year as f64;
    let mut mjd = ((year - 1859) * 365 + 14 + 31) as f64;
    mjd += frac * days_in_year(year) as f64;
    let mut y = year - 1;
    while y > 1858 {
        if is_leap_year(y) {
            mjd += 1.0;
        }
        y -= 1;
    }
    mjd
}

/// Convert a `YYYYMMDD` encoded date to Modified Julian Date.
///
/// Ported from src/conversions/unitconvert.cpp.
fn yyyymmdd_to_mjd(v: f64) -> f64 {
    let year = (v / 10000.0).floor() as i64;
    let month = ((v - year as f64 * 10000.0) / 100.0).floor() as i64;
    let day = (v - year as f64 * 10000.0 - month as f64 * 100.0).floor() as i64;
    let mut mjd = daynumber_in_year(year, month, day) as f64;
    let mut y = year - 1;
    while y > 1858 {
        mjd += days_in_year(y) as f64;
        y -= 1;
    }
    mjd + 13.0 + 31.0
}

/// Convert a Modified Julian Date to a decimal year.
///
/// Ported from src/conversions/unitconvert.cpp.
fn mjd_to_decimalyear(mjd: f64) -> f64 {
    let mut mjd_iter = 14.0 + 31.0;
    let mut year = 1859i64;
    while mjd >= mjd_iter {
        mjd_iter += days_in_year(year) as f64;
        year += 1;
    }
    year -= 1;
    mjd_iter -= days_in_year(year) as f64;
    year as f64 + (mjd - mjd_iter) / days_in_year(year) as f64
}

/// Convert a Modified Julian Date to a `YYYYMMDD` encoded date.
///
/// Ported from src/conversions/unitconvert.cpp.
fn mjd_to_yyyymmdd(mjd: f64) -> f64 {
    let date = mjd.round() as i64;
    let mut date_iter = 14 + 31;
    let mut year = 1859i64;
    while date >= date_iter {
        date_iter += days_in_year(year);
        year += 1;
    }
    year -= 1;
    date_iter -= days_in_year(year);
    let mut month = 1i64;
    while date_iter + days_in_month(year, month) <= date {
        date_iter += days_in_month(year, month);
        month += 1;
    }
    let day = date - date_iter + 1;
    year as f64 * 10000.0 + month as f64 * 100.0 + day as f64
}

impl TimeUnit {
    /// Convert a value expressed in `self` into Modified Julian Date.
    ///
    /// Ported from src/conversions/unitconvert.cpp.
    fn to_mjd(self, t: f64) -> f64 {
        match self {
            TimeUnit::Mjd => t,
            TimeUnit::DecimalYear => decimalyear_to_mjd(t),
            TimeUnit::GpsWeek => 44244.0 + t * 7.0,
            TimeUnit::YyyyMmDd => yyyymmdd_to_mjd(t),
        }
    }

    /// Convert a Modified Julian Date into a value expressed in `self`.
    ///
    /// Ported from src/conversions/unitconvert.cpp.
    fn convert_from_mjd(self, mjd: f64) -> f64 {
        match self {
            TimeUnit::Mjd => mjd,
            TimeUnit::DecimalYear => mjd_to_decimalyear(mjd),
            TimeUnit::GpsWeek => (mjd - 44244.0) / 7.0,
            TimeUnit::YyyyMmDd => mjd_to_yyyymmdd(mjd),
        }
    }
}

/// Whether a resolved unit identifier names a linear unit, an angular unit,
/// or is a bare numeric scale factor of unknown type.
///
/// Mirrors PROJ's `xy_in_is_linear`/`xy_out_is_linear` tri-state (`1` linear,
/// `0` angular, `-1` unknown) in `unitconvert.cpp`'s
/// `get_unit_conversion_factor`.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum UnitKind {
    Linear,
    Angular,
    /// A bare numeric factor (no named unit matched): PROJ leaves this
    /// "unknown" and never flags it as inconsistent against the other side.
    Unknown,
}

/// Resolve a unit identifier to a `(factor, optional name, kind)` triple.
///
/// Linear units are tried first, then angular units, then a bare numeric
/// scale factor. A zero, infinite or NaN numeric factor is rejected.
///
/// Ported from src/conversions/unitconvert.cpp.
fn resolve_factor(id: &str) -> ProjResult<(f64, Option<&'static str>, UnitKind)> {
    if let Some(u) = find_linear_unit(id) {
        return Ok((u.factor, Some(u.name), UnitKind::Linear));
    }
    if let Some(u) = find_angular_unit(id) {
        return Ok((u.factor, Some(u.name), UnitKind::Angular));
    }
    if let Ok(f) = id.parse::<f64>() {
        if f != 0.0 && f.is_finite() {
            return Ok((f, None, UnitKind::Unknown));
        }
    }
    Err(ProjError::IllegalArgValue)
}

/// The `unitconvert` conversion. Ported from src/conversions/unitconvert.cpp.
#[derive(Debug)]
struct UnitConvert {
    xy_factor: f64,
    z_factor: f64,
    t_in: Option<TimeUnit>,
    t_out: Option<TimeUnit>,
}

/// Ported from src/conversions/unitconvert.cpp.
impl Operation for UnitConvert {
    fn forward_2d(&self, lp: Lp) -> ProjResult<Xy> {
        Ok(Xy::new(lp.lam * self.xy_factor, lp.phi * self.xy_factor))
    }

    fn inverse_2d(&self, xy: Xy) -> ProjResult<Lp> {
        Ok(Lp::new(xy.x / self.xy_factor, xy.y / self.xy_factor))
    }

    fn forward_3d(&self, lpz: Lpz) -> ProjResult<Xyz> {
        Ok(Xyz::new(
            lpz.lam * self.xy_factor,
            lpz.phi * self.xy_factor,
            lpz.z * self.z_factor,
        ))
    }

    fn inverse_3d(&self, xyz: Xyz) -> ProjResult<Lpz> {
        Ok(Lpz::new(
            xyz.x / self.xy_factor,
            xyz.y / self.xy_factor,
            xyz.z / self.z_factor,
        ))
    }

    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let x = v[0] * self.xy_factor;
        let y = v[1] * self.xy_factor;
        let z = v[2] * self.z_factor;
        let mut t = v[3];
        if let Some(u) = self.t_in {
            t = u.to_mjd(t);
        }
        if let Some(u) = self.t_out {
            t = u.convert_from_mjd(t);
        }
        Ok(Coord::new(x, y, z, t))
    }

    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let x = v[0] / self.xy_factor;
        let y = v[1] / self.xy_factor;
        let z = v[2] / self.z_factor;
        let mut t = v[3];
        if let Some(u) = self.t_out {
            t = u.to_mjd(t);
        }
        if let Some(u) = self.t_in {
            t = u.convert_from_mjd(t);
        }
        Ok(Coord::new(x, y, z, t))
    }

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

/// Construct the `unitconvert` conversion. Ported from
/// src/conversions/unitconvert.cpp.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let pa = p.params;

    let mut left = IoUnits::Whatever;
    let mut right = IoUnits::Whatever;

    let mut xy_in_kind = UnitKind::Unknown;
    let mut xy_out_kind = UnitKind::Unknown;

    let mut xy_factor = 1.0;
    if let Some(s) = pa.get_str("xy_in") {
        let (f_in, name, kind) = resolve_factor(s)?;
        xy_factor = f_in;
        xy_in_kind = kind;
        if name == Some("Radian") {
            left = IoUnits::Radians;
        }
        if name == Some("Degree") {
            left = IoUnits::Degrees;
        }
    }
    if let Some(s) = pa.get_str("xy_out") {
        let (f_out, name, kind) = resolve_factor(s)?;
        xy_factor /= f_out;
        xy_out_kind = kind;
        if name == Some("Radian") {
            right = IoUnits::Radians;
        }
        if name == Some("Degree") {
            right = IoUnits::Degrees;
        }
    }
    // PROJ: "inconsistent unit type between xy_in and xy_out" — only flagged
    // when *both* sides resolved to a named (linear or angular) unit; a bare
    // numeric factor on either side is left untyped and never conflicts.
    if xy_in_kind != UnitKind::Unknown
        && xy_out_kind != UnitKind::Unknown
        && xy_in_kind != xy_out_kind
    {
        return Err(ProjError::IllegalArgValue);
    }

    let mut z_in_kind = UnitKind::Unknown;
    let mut z_out_kind = UnitKind::Unknown;

    let mut z_factor = 1.0;
    if let Some(s) = pa.get_str("z_in") {
        let (f_in, _, kind) = resolve_factor(s)?;
        z_factor = f_in;
        z_in_kind = kind;
    }
    if let Some(s) = pa.get_str("z_out") {
        let (f_out, _, kind) = resolve_factor(s)?;
        z_factor /= f_out;
        z_out_kind = kind;
    }
    // PROJ applies the identical "inconsistent unit type" check between
    // z_in and z_out.
    if z_in_kind != UnitKind::Unknown && z_out_kind != UnitKind::Unknown && z_in_kind != z_out_kind
    {
        return Err(ProjError::IllegalArgValue);
    }

    let t_in = match pa.get_str("t_in") {
        Some(s) => Some(parse_time_unit(s)?),
        None => None,
    };
    let t_out = match pa.get_str("t_out") {
        Some(s) => Some(parse_time_unit(s)?),
        None => None,
    };

    let op = UnitConvert {
        xy_factor,
        z_factor,
        t_in,
        t_out,
    };
    Ok(TransBuild::new(Box::new(op), left, right))
}

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

    #[derive(Default)]
    struct UcParams {
        xy_in: Option<String>,
        xy_out: Option<String>,
        z_in: Option<String>,
        z_out: Option<String>,
        t_in: Option<String>,
        t_out: Option<String>,
    }

    impl crate::TransParamLookup for UcParams {
        fn get_dms(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_f64(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_int(&self, _key: &str) -> Option<i64> {
            None
        }
        fn get_str(&self, key: &str) -> Option<&str> {
            match key {
                "xy_in" => self.xy_in.as_deref(),
                "xy_out" => self.xy_out.as_deref(),
                "z_in" => self.z_in.as_deref(),
                "z_out" => self.z_out.as_deref(),
                "t_in" => self.t_in.as_deref(),
                "t_out" => self.t_out.as_deref(),
                _ => None,
            }
        }
        fn get_bool(&self, _key: &str) -> bool {
            false
        }
        fn exists(&self, key: &str) -> bool {
            self.get_str(key).is_some()
        }
    }

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

    fn approx(a: f64, b: f64, tol: f64) -> bool {
        (a - b).abs() < tol
    }

    #[test]
    fn xy_meters_to_us_survey_feet() {
        let b = build(UcParams {
            xy_in: Some("m".into()),
            xy_out: Some("us-ft".into()),
            ..Default::default()
        });
        let op = b.operation;
        let fwd = op
            .forward_4d(Coord::new(1000000.0, 2000000.0, 0.0, 0.0))
            .unwrap();
        let v = fwd.v();
        assert!(approx(v[0], 3280833.333333333, 1e-6), "x = {}", v[0]);
        assert!(approx(v[1], 6561666.666666667, 1e-6), "y = {}", v[1]);

        let inv = op.inverse_4d(fwd).unwrap();
        let iv = inv.v();
        assert!(approx(iv[0], 1000000.0, 1e-6), "inv x = {}", iv[0]);
        assert!(approx(iv[1], 2000000.0, 1e-6), "inv y = {}", iv[1]);
    }

    #[test]
    fn z_meters_to_feet() {
        let b = build(UcParams {
            z_in: Some("m".into()),
            z_out: Some("ft".into()),
            ..Default::default()
        });
        let op = b.operation;
        let fwd = op.forward_4d(Coord::new(0.0, 0.0, 100.0, 0.0)).unwrap();
        let v = fwd.v();
        assert!(approx(v[2], 328.083989501, 1e-6), "z = {}", v[2]);

        let inv = op.inverse_4d(fwd).unwrap();
        assert!(approx(inv.v()[2], 100.0, 1e-6), "inv z = {}", inv.v()[2]);
    }

    #[test]
    fn t_decimalyear_to_mjd() {
        let b = build(UcParams {
            t_in: Some("decimalyear".into()),
            t_out: Some("mjd".into()),
            ..Default::default()
        });
        let op = b.operation;
        let fwd = op.forward_4d(Coord::new(0.0, 0.0, 0.0, 2017.5)).unwrap();
        let v = fwd.v();
        assert!(approx(v[3], 57936.5, 1e-3), "t = {}", v[3]);

        let inv = op.inverse_4d(fwd).unwrap();
        assert!(approx(inv.v()[3], 2017.5, 1e-6), "inv t = {}", inv.v()[3]);
    }

    #[test]
    fn combined_xy_z_t() {
        let b = build(UcParams {
            xy_in: Some("km".into()),
            xy_out: Some("m".into()),
            z_in: Some("m".into()),
            z_out: Some("ft".into()),
            t_in: Some("decimalyear".into()),
            t_out: Some("mjd".into()),
        });
        let op = b.operation;
        let fwd = op
            .forward_4d(Coord::new(1000.0, 2000.0, 100.0, 2017.5))
            .unwrap();
        let v = fwd.v();
        assert!(approx(v[0], 1000000.0, 1e-6), "x = {}", v[0]);
        assert!(approx(v[1], 2000000.0, 1e-6), "y = {}", v[1]);
        assert!(approx(v[2], 328.083989501, 1e-6), "z = {}", v[2]);
        assert!(approx(v[3], 57936.5, 1e-3), "t = {}", v[3]);
    }

    #[test]
    fn unknown_time_unit_errors() {
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = UcParams {
            t_in: Some("nope".into()),
            ..Default::default()
        };
        let tp = TransParams {
            ellipsoid: &ell,
            params: &params,
            registry: None,
        };
        assert_eq!(new(&tp).err(), Some(ProjError::IllegalArgValue));
    }

    #[test]
    fn unknown_xy_unit_errors() {
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = UcParams {
            xy_in: Some("frobnitz".into()),
            ..Default::default()
        };
        let tp = TransParams {
            ellipsoid: &ell,
            params: &params,
            registry: None,
        };
        assert_eq!(new(&tp).err(), Some(ProjError::IllegalArgValue));
    }

    /// PROJ (verified against `cs2cs +proj=unitconvert +xy_in=rad
    /// +xy_out=m`, error 1027 "inconsistent unit type between xy_in and
    /// xy_out") rejects mixing an angular `xy_in` with a linear `xy_out`
    /// rather than silently producing a meaningless scale factor.
    #[test]
    fn mixed_angular_linear_xy_units_errors() {
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = UcParams {
            xy_in: Some("rad".into()),
            xy_out: Some("m".into()),
            ..Default::default()
        };
        let tp = TransParams {
            ellipsoid: &ell,
            params: &params,
            registry: None,
        };
        assert_eq!(new(&tp).err(), Some(ProjError::IllegalArgValue));
    }

    /// Same check, reversed direction (linear in, angular out).
    #[test]
    fn mixed_linear_angular_xy_units_errors() {
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = UcParams {
            xy_in: Some("m".into()),
            xy_out: Some("deg".into()),
            ..Default::default()
        };
        let tp = TransParams {
            ellipsoid: &ell,
            params: &params,
            registry: None,
        };
        assert_eq!(new(&tp).err(), Some(ProjError::IllegalArgValue));
    }

    /// PROJ applies the identical "inconsistent unit type" check between
    /// z_in and z_out (verified against `cs2cs +proj=unitconvert +z_in=rad
    /// +z_out=m`, error 1027 "inconsistent unit type between z_in and
    /// z_out").
    #[test]
    fn mixed_angular_linear_z_units_errors() {
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = UcParams {
            z_in: Some("rad".into()),
            z_out: Some("m".into()),
            ..Default::default()
        };
        let tp = TransParams {
            ellipsoid: &ell,
            params: &params,
            registry: None,
        };
        assert_eq!(new(&tp).err(), Some(ProjError::IllegalArgValue));
    }

    /// Both angular is fine (same kind on both sides).
    #[test]
    fn both_angular_xy_units_ok() {
        let b = build(UcParams {
            xy_in: Some("rad".into()),
            xy_out: Some("deg".into()),
            ..Default::default()
        });
        let op = b.operation;
        let fwd = op
            .forward_4d(Coord::new(std::f64::consts::PI, 0.0, 0.0, 0.0))
            .unwrap();
        assert!(approx(fwd.v()[0], 180.0, 1e-9), "x = {}", fwd.v()[0]);
    }

    /// A bare numeric factor on one side must remain untyped/allowed against
    /// either a linear or an angular unit on the other side (PROJ only
    /// compares when *both* sides resolved to a named unit).
    #[test]
    fn bare_numeric_factor_against_linear_unit_ok() {
        assert!(resolve_factor("2.5").is_ok());
        let b = build(UcParams {
            xy_in: Some("2.0".into()),
            xy_out: Some("m".into()),
            ..Default::default()
        });
        let op = b.operation;
        let fwd = op.forward_4d(Coord::new(3.0, 0.0, 0.0, 0.0)).unwrap();
        assert!(approx(fwd.v()[0], 6.0, 1e-12), "x = {}", fwd.v()[0]);
    }

    /// A bare numeric factor on one side must also remain allowed against an
    /// angular unit on the other side.
    #[test]
    fn bare_numeric_factor_against_angular_unit_ok() {
        let b = build(UcParams {
            xy_in: Some("2.0".into()),
            xy_out: Some("rad".into()),
            ..Default::default()
        });
        assert!(b
            .operation
            .forward_4d(Coord::new(1.0, 0.0, 0.0, 0.0))
            .is_ok());
    }
}