oxiproj-transformations 0.1.2

#![forbid(unsafe_code)]
//! TIN-based coordinate shift (`tinshift`) — port of PROJ `src/transformations/tinshift.cpp`.
//! Reads a PROJ `triangulation_file` JSON grid and interpolates coordinate shifts using
//! barycentric interpolation within triangles of an irregular triangulated network.
//!
//! Both directions are supported, matching PROJ's `Evaluator::forward` /
//! `Evaluator::inverse`: the forward transform searches the triangulation in
//! source space and interpolates the target coordinates, while the inverse
//! searches in target space and interpolates the source coordinates. Coordinates
//! pass through in the grid's native units (PROJ's `PJ_IO_UNITS_WHATEVER`), so no
//! radian/degree scaling is applied here.

use crate::{TransBuild, TransParams};
use oxiproj_core::{Coord, IoUnits, Operation, ProjError, ProjResult};
use oxiproj_grids::{read_tinshift, TinShiftIndex};

#[derive(Debug)]
struct TinShiftTransform {
    index: TinShiftIndex,
}

impl Operation for TinShiftTransform {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let (x, y, z) = self
            .index
            .forward(v[0], v[1], v[2])
            .ok_or(ProjError::OutsideGrid)?;
        Ok(Coord::new(x, y, z, v[3]))
    }

    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        let (x, y, z) = self
            .index
            .inverse(v[0], v[1], v[2])
            .ok_or(ProjError::OutsideGrid)?;
        Ok(Coord::new(x, y, z, v[3]))
    }

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

/// Construct a `tinshift` transform from parsed parameters.
///
/// Requires a grid name (`+grids=<name>`, or PROJ's `+file=<name>`) resolvable
/// through the [`GridRegistry`](crate::GridRegistry) in `p.registry`. The named
/// grid is parsed as a PROJ `triangulation_file` JSON document.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let grid_name = p
        .params
        .get_str("grids")
        .or_else(|| p.params.get_str("file"))
        .ok_or(ProjError::MissingArg)?;
    let registry = p.registry.ok_or(ProjError::FileNotFound)?;
    let data = registry
        .get_grid(grid_name)
        .ok_or(ProjError::FileNotFound)?;
    let index = read_tinshift(data)?;
    Ok(TransBuild::new(
        Box::new(TinShiftTransform { index }),
        IoUnits::Whatever,
        IoUnits::Whatever,
    ))
}

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

    /// Minimal horizontal triangulation matching PROJ's
    /// `test/unit/test_tinshift.cpp` `getMinValidContent()`.
    /// Source vertices (0,0),(0,1),(1,1) map to target
    /// (101,101),(100,101),(100,100); single triangle [0,1,2].
    fn make_horizontal_json() -> Vec<u8> {
        let s = r#"{
  "file_type": "triangulation_file",
  "format_version": "1.0",
  "input_crs": "EPSG:2393",
  "output_crs": "EPSG:3067",
  "transformed_components": [ "horizontal" ],
  "vertices_columns": [ "source_x", "source_y", "target_x", "target_y" ],
  "triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ],
  "vertices": [ [0, 0, 101, 101], [0, 1, 100, 101], [1, 1, 100, 100] ],
  "triangles": [ [0, 1, 2] ]
}"#;
        s.as_bytes().to_vec()
    }

    // Copied from PROJ 9.8.0 data/tests/tinshift_simplified_n60_n2000.json
    // (MIT/X11; Copyright (c) 2020, Even Rouault). Vertical-only transform.
    fn make_vertical_json() -> Vec<u8> {
        let s = r#"{
  "file_type": "triangulation_file",
  "format_version": "1.0",
  "input_crs": "EPSG:2393+5717",
  "output_crs": "EPSG:2393+5941",
  "transformed_components": [ "vertical" ],
  "vertices_columns": [ "source_x", "source_y", "source_z", "target_z" ],
  "triangles_columns": [ "idx_vertex1", "idx_vertex2", "idx_vertex3" ],
  "vertices": [ [3188607.0, 6688748.0, 23.123, 23.4133],
                [3184981.0, 6725255.0, 8.044, 8.34499],
                [3220912.0, 6699508.0, 1.724, 2.0101] ],
  "triangles": [ [0, 1, 2] ]
}"#;
        s.as_bytes().to_vec()
    }

    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_transform(data: Vec<u8>) -> TransBuild {
        let mut map = std::collections::HashMap::new();
        map.insert("test.json".to_string(), data);
        let reg = InMemReg(map);
        let tp = TestParams {
            grids: "test.json".to_string(),
        };
        let ell = oxiproj_core::Ellipsoid::named("WGS84").expect("WGS84 ellipsoid");
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        new(&params).expect("build tinshift")
    }

    #[test]
    fn declares_whatever_units() {
        // PROJ's tinshift is PJ_IO_UNITS_WHATEVER on both sides: coordinates
        // pass through raw in the grid's native units.
        let tb = make_transform(make_horizontal_json());
        assert_eq!(tb.left, IoUnits::Whatever);
        assert_eq!(tb.right, IoUnits::Whatever);
    }

    #[test]
    fn forward_matches_proj() {
        // PROJ unit test: forward(0.5, 0.75) -> (100.25, 100.5); z passes through.
        let tb = make_transform(make_horizontal_json());
        let input = Coord::new(0.5, 0.75, 1000.0, 0.0);
        let v = tb.operation.forward_4d(input).expect("inside mesh").v();
        assert!((v[0] - 100.25).abs() < 1e-9, "x: got {}", v[0]);
        assert!((v[1] - 100.5).abs() < 1e-9, "y: got {}", v[1]);
        assert_eq!(v[2], 1000.0, "z passes through");
    }

    #[test]
    fn inverse_is_available_and_round_trips() {
        // Regression: inverse was previously unimplemented (NoInverseOp).
        let tb = make_transform(make_horizontal_json());
        assert!(tb.operation.has_inverse());

        let input = Coord::new(0.5, 0.75, 1000.0, 0.0);
        let fwd = tb.operation.forward_4d(input).expect("forward");
        let back = tb.operation.inverse_4d(fwd).expect("inverse");
        let v = back.v();
        assert!((v[0] - 0.5).abs() < 1e-9, "x round-trip: got {}", v[0]);
        assert!((v[1] - 0.75).abs() < 1e-9, "y round-trip: got {}", v[1]);
        assert_eq!(v[2], 1000.0);
    }

    #[test]
    fn vertical_forward_and_inverse() {
        // Verified against PROJ cct: (3210000, 6700000, 10.0) -> z 10.28863673.
        let tb = make_transform(make_vertical_json());
        let input = Coord::new(3_210_000.0, 6_700_000.0, 10.0, 0.0);
        let fwd = tb.operation.forward_4d(input).expect("forward");
        let v = fwd.v();
        assert_eq!(v[0], 3_210_000.0);
        assert_eq!(v[1], 6_700_000.0);
        assert!((v[2] - 10.288_636_73).abs() < 1e-6, "z: got {}", v[2]);

        let back = tb.operation.inverse_4d(fwd).expect("inverse");
        assert!(
            (back.v()[2] - 10.0).abs() < 1e-6,
            "z round-trip: got {}",
            back.v()[2]
        );
    }

    #[test]
    fn outside_mesh_returns_err() {
        let tb = make_transform(make_horizontal_json());
        let far = Coord::new(-10.0, -10.0, 0.0, 0.0);
        assert_eq!(
            tb.operation.forward_4d(far).err(),
            Some(ProjError::OutsideGrid)
        );
        assert_eq!(
            tb.operation.inverse_4d(far).err(),
            Some(ProjError::OutsideGrid)
        );
    }

    #[test]
    fn missing_grids_param() {
        struct NoParams;
        impl crate::TransParamLookup for NoParams {
            fn get_str(&self, _: &str) -> Option<&str> {
                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, _: &str) -> bool {
                false
            }
        }
        let reg = InMemReg(std::collections::HashMap::new());
        let ell = oxiproj_core::Ellipsoid::named("WGS84").expect("WGS84 ellipsoid");
        let np = NoParams;
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &np,
            registry: Some(&reg),
        };
        assert_eq!(new(&params).err(), Some(ProjError::MissingArg));
    }
}