oxiproj-transformations 0.1.1

#![forbid(unsafe_code)]
//! Vertical grid shift (`vgridshift`) — port of PROJ `src/transformations/vgridshift.cpp`.
//! Reads a GTX binary grid and applies bilinear-interpolated vertical (z) shifts in meters.

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

#[derive(Debug)]
struct VGridShift {
    grid: GridSet,
}

impl Operation for VGridShift {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let lon_deg = v[0].to_degrees();
        let lat_deg = v[1].to_degrees();
        let shifts = sample_grid(&self.grid, lat_deg, lon_deg).ok_or(ProjError::OutsideGrid)?;
        Ok(Coord::new(v[0], v[1], v[2] + shifts[0], v[3]))
    }

    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        // For vertical-only shift the horizontal position is unchanged in both CRS,
        // so we can directly look up the shift and subtract it.
        let v = c.v();
        let lon_deg = v[0].to_degrees();
        let lat_deg = v[1].to_degrees();
        let shifts = sample_grid(&self.grid, lat_deg, lon_deg).ok_or(ProjError::OutsideGrid)?;
        Ok(Coord::new(v[0], v[1], v[2] - shifts[0], v[3]))
    }

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

/// Construct a `vgridshift` transform from parsed parameters.
///
/// Requires `+grids=<name>` and a populated [`GridRegistry`](crate::GridRegistry).
/// The named grid is parsed as GTX (single-band, big-endian).
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let grid_name = p.params.get_str("grids").ok_or(ProjError::MissingArg)?;
    let registry = p.registry.ok_or(ProjError::FileNotFound)?;
    let data = registry
        .get_grid(grid_name)
        .ok_or(ProjError::FileNotFound)?;
    let grid = read_gtx(data, grid_name)?;
    if grid.bands.is_empty() {
        return Err(ProjError::FileNotFound);
    }
    Ok(TransBuild::new(
        Box::new(VGridShift { grid }),
        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())
        }
    }

    struct TestParams {
        grids: String,
    }
    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, _: &str) -> Option<f64> {
            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"
        }
    }

    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(),
        };
        (reg, ell, tp)
    }

    #[test]
    fn test_vgridshift_forward_adds_shift() {
        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");
    }

    #[test]
    fn test_vgridshift_inverse_subtracts_shift() {
        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, 10.0, 0.0);
        let out = tb.operation.inverse_4d(input).unwrap();
        assert!(
            (out.v()[2] - 0.0).abs() < 1e-6,
            "z should be 0 after inverse"
        );
    }

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