oxiproj-transformations 0.1.2

#![forbid(unsafe_code)]
//! Vertical grid shift (`vgridshift`) — port of PROJ `src/transformations/vgridshift.cpp`.
//!
//! Reads a GTX binary grid **or** a GeoTIFF vertical/geoid grid (band role
//! `geoid_undulation` / `vertical_offset`, in metres) and applies
//! bilinear-interpolated vertical (z) shifts. Honors `+multiplier` (with the
//! legacy VERTCON `.gtx`→`.tif` millimetre hack), the `+t_epoch`/`+t_final`
//! time-restriction gate, and PROJ's comma-separated `+grids=` fallback list
//! (including the per-entry `@` optional prefix).

use crate::{TransBuild, TransParams};
use oxiproj_core::{Coord, IoUnits, Operation, ProjError, ProjResult};
use oxiproj_grids::{read_geotiff_hierarchy, read_gtx, sample_grid, GridSet};

use super::gridshift::{
    is_tiff, parse_time_gate, resolve_grid_list, time_gate_applies, vertical_gridset,
};

#[derive(Debug)]
struct VGridShift {
    /// Fallback list of vertical grids (single band each, metres); the first
    /// grid containing the point wins, matching PROJ's `pj_vgrid_value`.
    grids: Vec<GridSet>,
    /// `+multiplier`: the grid value is scaled by this before being added to z
    /// in the forward direction (PROJ default `-1.0`, i.e. z -= shift).
    multiplier: f64,
    /// `+t_epoch` (0 disables gating).
    t_epoch: f64,
    /// `+t_final` (0 disables gating).
    t_final: f64,
}

impl VGridShift {
    /// Interpolated grid value at `(lat_deg, lon_deg)`, scaled by `+multiplier`,
    /// searching the fallback list. Mirrors PROJ's `pj_vgrid_value`.
    fn value(&self, lat_deg: f64, lon_deg: f64) -> ProjResult<f64> {
        for gs in &self.grids {
            if let Some(shifts) = sample_grid(gs, lat_deg, lon_deg) {
                return Ok(shifts[0] * self.multiplier);
            }
        }
        Err(ProjError::OutsideGrid)
    }
}

impl Operation for VGridShift {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        // PROJ 9.8.0 `vgridshift.cpp`: `point.xyz.z += pj_vgrid_value(...,
        // forward_multiplier)`, where the value already includes the multiplier
        // (default `-1.0`, so forward SUBTRACTS the interpolated grid value).
        let v = c.v();
        // Empty grid list (all-optional grids absent): pass through unchanged.
        if self.grids.is_empty() {
            return Ok(c);
        }
        // Time-restriction gate: outside the bracket the point is unchanged.
        if !time_gate_applies(self.t_epoch, self.t_final, v[3]) {
            return Ok(c);
        }
        let value = self.value(v[1].to_degrees(), v[0].to_degrees())?;
        Ok(Coord::new(v[0], v[1], v[2] + value, v[3]))
    }

    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        // Mirror of forward: PROJ's `pj_vgridshift_reverse_3d` does
        // `point.xyz.z -= pj_vgrid_value(..., forward_multiplier)`. For a
        // vertical-only shift the horizontal position is unchanged.
        let v = c.v();
        if self.grids.is_empty() {
            return Ok(c);
        }
        if !time_gate_applies(self.t_epoch, self.t_final, v[3]) {
            return Ok(c);
        }
        let value = self.value(v[1].to_degrees(), v[0].to_degrees())?;
        Ok(Coord::new(v[0], v[1], v[2] - value, v[3]))
    }

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

/// Names of the legacy VERTCON GTX grids that PROJ re-distributes as GeoTIFF;
/// the `.gtx` files stored millimetres, the `.tif` files metres, so a
/// `+multiplier=0.001` (mm→m) must be reset to 1.0 when the resolved grid is a
/// GeoTIFF (PROJ `deal_with_vertcon_gtx_hack`).
const VERTCON_GTX: [&str; 3] = ["vertconw.gtx", "vertconc.gtx", "vertcone.gtx"];

/// Construct a `vgridshift` transform from parsed parameters.
///
/// Requires `+grids=<name>[,<name>...]` (comma-separated, `@`-optional entries).
/// Each grid is resolved registry-first, then from the local
/// `PROJ_DATA`/`PROJ_LIB` directories and the unified disk cache, then (with the
/// `network` feature) the PROJ CDN — see `super::gridshift::resolve_grid_bytes`.
/// GeoTIFF grids (TIFF magic) are decoded and reduced to their
/// geoid/vertical-offset band; everything else is parsed as GTX (single-band,
/// big-endian, metres). Honors `+multiplier`, `+t_epoch`/`+t_final`.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let grid_name = p.params.get_str("grids").ok_or(ProjError::MissingArg)?;

    // PROJ: historical default forward_multiplier = -1.0, overridable by
    // +multiplier.
    let mut multiplier = p.params.get_f64("multiplier").unwrap_or(-1.0);
    let (t_epoch, t_final) = parse_time_gate(p);

    let resolved = resolve_grid_list(p.registry, grid_name)?;

    let mut grids = Vec::with_capacity(resolved.len());
    let mut first_is_tiff = false;
    for (idx, (name, bytes)) in resolved.iter().enumerate() {
        if is_tiff(bytes) {
            if idx == 0 {
                first_is_tiff = true;
            }
            // Multi-IFD GeoTIFF vertical grids can carry a parent/child sub-grid
            // hierarchy just like horizontal ones; read every georeferenced IFD
            // ordered child-first so the first-containing-grid-wins scan picks
            // the finest sub-grid covering the point (PROJ `GTiffVGridShiftSet`).
            for gs in read_geotiff_hierarchy(bytes, name)? {
                let gs = vertical_gridset(gs)?;
                if gs.bands.is_empty() {
                    return Err(ProjError::FileNotFound);
                }
                grids.push(gs);
            }
        } else {
            let gs = read_gtx(bytes, name)?;
            if gs.bands.is_empty() {
                return Err(ProjError::FileNotFound);
            }
            grids.push(gs);
        }
    }

    // VERTCON .gtx→.tif millimetre hack: a +multiplier=0.001 supplied for the
    // legacy mm-stored VERTCON grids must be reset to 1.0 when the actual grid
    // resolved to a (metre-stored) GeoTIFF.
    if multiplier == 0.001 && first_is_tiff && VERTCON_GTX.contains(&grid_name) {
        multiplier = 1.0;
    }

    Ok(TransBuild::new(
        Box::new(VGridShift {
            grids,
            multiplier,
            t_epoch,
            t_final,
        }),
        IoUnits::Radians,
        IoUnits::Radians,
    ))
}

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

    fn build_gtx_uniform(shift: f32) -> Vec<u8> {
        // 2×2 GTX grid covering lat 0-1 deg, lon 0-1 deg
        let mut buf = Vec::new();
        buf.extend_from_slice(&0.0f64.to_be_bytes()); // south_lat
        buf.extend_from_slice(&0.0f64.to_be_bytes()); // west_lon
        buf.extend_from_slice(&1.0f64.to_be_bytes()); // lat_inc
        buf.extend_from_slice(&1.0f64.to_be_bytes()); // lon_inc
        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.to_be_bytes());
        }
        buf
    }

    struct InMemReg(std::collections::HashMap<String, Vec<u8>>);
    impl crate::GridRegistry for InMemReg {
        fn get_grid(&self, name: &str) -> Option<&[u8]> {
            self.0.get(name).map(|v| v.as_slice())
        }
    }

    #[derive(Default)]
    struct TestParams {
        grids: String,
        multiplier: Option<f64>,
        t_epoch: Option<f64>,
        t_final: Option<f64>,
    }
    impl crate::TransParamLookup for TestParams {
        fn get_str(&self, key: &str) -> Option<&str> {
            if key == "grids" {
                Some(&self.grids)
            } else {
                None
            }
        }
        fn get_f64(&self, key: &str) -> Option<f64> {
            match key {
                "multiplier" => self.multiplier,
                "t_epoch" => self.t_epoch,
                "t_final" => self.t_final,
                _ => None,
            }
        }
        fn get_dms(&self, _: &str) -> Option<f64> {
            None
        }
        fn get_int(&self, _: &str) -> Option<i64> {
            None
        }
        fn get_bool(&self, _: &str) -> bool {
            false
        }
        fn exists(&self, key: &str) -> bool {
            key == "grids"
                || (key == "multiplier" && self.multiplier.is_some())
                || (key == "t_epoch" && self.t_epoch.is_some())
                || (key == "t_final" && self.t_final.is_some())
        }
    }

    fn make_params(data: Vec<u8>) -> (InMemReg, oxiproj_core::Ellipsoid, TestParams) {
        let mut map = std::collections::HashMap::new();
        map.insert("test.gtx".to_string(), data);
        let reg = InMemReg(map);
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let tp = TestParams {
            grids: "test.gtx".to_string(),
            ..Default::default()
        };
        (reg, ell, tp)
    }

    #[test]
    fn test_vgridshift_forward_subtracts_shift() {
        // End-to-end parity with PROJ 9.7.0 `cct` on a uniform +10 m GTX grid:
        //   echo "0.5 0.5 0" | cct +proj=vgridshift +grids=geoid10.gtx
        //   ->  0.5  0.5  -10.0000
        // PROJ's forward SUBTRACTS the grid value (forward_multiplier = -1.0).
        let (reg, ell, tp) = make_params(build_gtx_uniform(10.0));
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&p).unwrap();
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);
        let out = tb.operation.forward_4d(input).unwrap();
        assert!(
            (out.v()[2] - (-10.0)).abs() < 1e-6,
            "z should be -10 (PROJ cct forward), got {}",
            out.v()[2]
        );
    }

    #[test]
    fn test_vgridshift_inverse_adds_shift() {
        // End-to-end parity with PROJ 9.7.0 `cct -I` on the same grid:
        //   echo "0.5 0.5 0" | cct -I +proj=vgridshift +grids=geoid10.gtx
        //   ->  0.5  0.5  +10.0000
        // PROJ's inverse ADDS the grid value.
        let (reg, ell, tp) = make_params(build_gtx_uniform(10.0));
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&p).unwrap();
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);
        let out = tb.operation.inverse_4d(input).unwrap();
        assert!(
            (out.v()[2] - 10.0).abs() < 1e-6,
            "z should be +10 (PROJ cct inverse), got {}",
            out.v()[2]
        );
    }

    #[test]
    fn test_vgridshift_outside_returns_err() {
        let (reg, ell, tp) = make_params(build_gtx_uniform(5.0));
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&p).unwrap();
        let far = Coord::new(45.0 * DEG_TO_RAD, 45.0 * DEG_TO_RAD, 0.0, 0.0);
        assert_eq!(
            tb.operation.forward_4d(far).err(),
            Some(ProjError::OutsideGrid)
        );
    }

    #[test]
    fn test_vgridshift_round_trip() {
        let (reg, ell, tp) = make_params(build_gtx_uniform(7.5));
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&p).unwrap();
        let input = Coord::new(0.3 * DEG_TO_RAD, 0.4 * DEG_TO_RAD, 100.0, 2020.0);
        let fwd = tb.operation.forward_4d(input).unwrap();
        let inv = tb.operation.inverse_4d(fwd).unwrap();
        let vi = inv.v();
        let vi0 = input.v();
        assert!((vi[2] - vi0[2]).abs() < 1e-9, "z round-trip error");
    }

    #[test]
    fn multiplier_scales_and_signs_the_shift() {
        // +multiplier=1.0: forward ADDS the grid value (opposite of the default
        // -1.0), matching PROJ `pj_vgrid_value` (value * multiplier).
        let (reg, ell, mut tp) = make_params(build_gtx_uniform(10.0));
        tp.multiplier = Some(1.0);
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&p).unwrap();
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);
        let out = tb.operation.forward_4d(input).unwrap();
        assert!(
            (out.v()[2] - 10.0).abs() < 1e-6,
            "z should be +10 with multiplier=+1, got {}",
            out.v()[2]
        );
        // A non-unit multiplier scales the shift (vertcon mm→m style, e.g. 0.5).
        let (reg2, ell2, mut tp2) = make_params(build_gtx_uniform(10.0));
        tp2.multiplier = Some(0.5);
        let p2 = crate::TransParams {
            ellipsoid: &ell2,
            params: &tp2,
            registry: Some(&reg2),
        };
        let tb2 = new(&p2).unwrap();
        let out2 = tb2.operation.forward_4d(input).unwrap();
        assert!(
            (out2.v()[2] - 5.0).abs() < 1e-6,
            "z should be +5 with multiplier=0.5, got {}",
            out2.v()[2]
        );
    }

    #[test]
    fn time_gate_skips_points_outside_bracket() {
        // +t_epoch=2000, +t_final=2010: only points with t < 2000 are shifted.
        let (reg, ell, mut tp) = make_params(build_gtx_uniform(10.0));
        tp.t_epoch = Some(2000.0);
        tp.t_final = Some(2010.0);
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&p).unwrap();
        // t = 1995 < t_epoch and t_final > t_epoch => shift applied (z -= 10).
        let inside = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 1995.0);
        let out_in = tb.operation.forward_4d(inside).unwrap();
        assert!(
            (out_in.v()[2] - (-10.0)).abs() < 1e-6,
            "in-bracket z should be -10, got {}",
            out_in.v()[2]
        );
        // t = 2005 >= t_epoch => point passes through unchanged.
        let outside = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 2005.0);
        let out_out = tb.operation.forward_4d(outside).unwrap();
        assert!(
            (out_out.v()[2] - 0.0).abs() < 1e-12,
            "out-of-bracket z unchanged, got {}",
            out_out.v()[2]
        );
    }

    #[test]
    fn comma_split_optional_prefix_falls_back() {
        // `+grids=@missing.gtx,test.gtx`: the leading `@` marks the first entry
        // optional, so resolution silently skips it and uses the second grid.
        let (reg, ell, mut tp) = make_params(build_gtx_uniform(10.0));
        tp.grids = "@missing.gtx,test.gtx".to_string();
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&p).unwrap();
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 0.0, 0.0);
        let out = tb.operation.forward_4d(input).unwrap();
        assert!(
            (out.v()[2] - (-10.0)).abs() < 1e-6,
            "fell back to second grid, z should be -10, got {}",
            out.v()[2]
        );
    }

    #[test]
    fn required_missing_grid_errors() {
        // A non-`@` (required) missing grid aborts the whole list.
        let (reg, ell, mut tp) = make_params(build_gtx_uniform(10.0));
        tp.grids = "missing.gtx,test.gtx".to_string();
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        assert_eq!(new(&p).err(), Some(ProjError::FileNotFound));
    }

    #[test]
    fn all_optional_absent_is_pass_through() {
        // Every entry optional and absent => empty grid list => identity.
        let (reg, ell, mut tp) = make_params(build_gtx_uniform(10.0));
        tp.grids = "@nope1.gtx,@nope2.gtx".to_string();
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&p).unwrap();
        let input = Coord::new(0.5 * DEG_TO_RAD, 0.5 * DEG_TO_RAD, 3.0, 0.0);
        let out = tb.operation.forward_4d(input).unwrap();
        assert!((out.v()[2] - 3.0).abs() < 1e-12, "pass-through z unchanged");
    }
}