oxiproj-transformations 0.1.2

#![forbid(unsafe_code)]
//! Horizontal grid shift (`hgridshift`) — port of PROJ `src/transformations/hgridshift.cpp`.
//!
//! Reads an NTv2 binary grid (lat/lon shifts in arc-seconds) **or** a GeoTIFF
//! horizontal-offset grid (semantics-aware role/unit/positive per GDAL_METADATA)
//! and applies bilinear-interpolated lat/lon shifts to input coordinates in
//! radians. Honors the `+t_epoch`/`+t_final` time-restriction gate and PROJ's
//! comma-separated `+grids=` fallback list (per-entry `@` optional prefix).

use crate::{TransBuild, TransParams};
use oxiproj_core::{Coord, IoUnits, Operation, ProjError, ProjResult};
use oxiproj_grids::{read_geotiff_hierarchy, read_ntv2, sample_grid, GridSet};
use std::f64::consts::PI;

use super::gridshift::{
    geotiff_horizontal_offset, is_tiff, parse_time_gate, resolve_grid_list, resolve_hshift_plan,
    time_gate_applies, HShiftPlan,
};

/// Arc-seconds to radians: π / (180 × 3600) = π / 648000
const ARCSEC_TO_RAD: f64 = PI / 648_000.0;
const MAX_ITER: usize = 10;
const ITER_TOL: f64 = 1e-12;

#[derive(Debug)]
struct HGridShift {
    /// Fallback list of horizontal grids; the first grid containing the point
    /// wins (PROJ `pj_hgrid_apply`).
    grids: Vec<GridSet>,
    /// Semantics-aware GeoTIFF recipe; `Some` for GeoTIFF grids, `None` for the
    /// classic NTv2 arc-second convention.
    hplan: Option<HShiftPlan>,
    /// `+t_epoch` (0 disables gating).
    t_epoch: f64,
    /// `+t_final` (0 disables gating).
    t_final: f64,
}

impl HGridShift {
    /// Longitude/latitude offset (radians) of the interpolated `shifts`, per the
    /// grid convention (NTv2 arc-second, or the GeoTIFF [`HShiftPlan`]).
    fn offset(&self, shifts: &[f64]) -> (f64, f64) {
        match &self.hplan {
            Some(plan) => geotiff_horizontal_offset(plan, shifts),
            // NTv2: shifts[0]=lat shift arcsec, shifts[1]=lon shift arcsec
            // (read_ntv2 already normalizes longitude to east-positive).
            None => (shifts[1] * ARCSEC_TO_RAD, shifts[0] * ARCSEC_TO_RAD),
        }
    }
}

impl Operation for HGridShift {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        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 lon_deg = v[0].to_degrees();
        let lat_deg = v[1].to_degrees();
        for gs in &self.grids {
            if let Some(shifts) = sample_grid(gs, lat_deg, lon_deg) {
                let (dlon, dlat) = self.offset(&shifts);
                return Ok(Coord::new(v[0] + dlon, v[1] + dlat, v[2], v[3]));
            }
        }
        Err(ProjError::OutsideGrid)
    }

    fn inverse_4d(&self, c: Coord) -> ProjResult<Coord> {
        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 mut lon_r = v[0];
        let mut lat_r = v[1];
        for gs in &self.grids {
            // Check if the point falls in this grid (approximate).
            if sample_grid(gs, lat_r.to_degrees(), lon_r.to_degrees()).is_none() {
                continue;
            }
            // Iterative fixed-point inversion.
            for _ in 0..MAX_ITER {
                let shifts = match sample_grid(gs, lat_r.to_degrees(), lon_r.to_degrees()) {
                    Some(s) => s,
                    None => break,
                };
                let (dlon, dlat) = self.offset(&shifts);
                let new_lon = v[0] - dlon;
                let new_lat = v[1] - dlat;
                let dlon_step = (new_lon - lon_r).abs();
                let dlat_step = (new_lat - lat_r).abs();
                lon_r = new_lon;
                lat_r = new_lat;
                if dlon_step < ITER_TOL && dlat_step < ITER_TOL {
                    break;
                }
            }
            return Ok(Coord::new(lon_r, lat_r, v[2], v[3]));
        }
        Err(ProjError::OutsideGrid)
    }

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

/// Construct an `hgridshift` 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, e.g. the NADCON5-derived `us_noaa_pahpgn.tif`) are
/// applied with role/arc-second/positive-aware semantics; NTv2 `.gsb` grids use
/// the classic arc-second convention. Honors `+t_epoch`/`+t_final`.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let grid_name = p.params.get_str("grids").ok_or(ProjError::MissingArg)?;
    let (t_epoch, t_final) = parse_time_gate(p);

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

    // Empty list (all-optional grids absent): a pass-through transform.
    if resolved.is_empty() {
        return Ok(TransBuild::new(
            Box::new(HGridShift {
                grids: Vec::new(),
                hplan: None,
                t_epoch,
                t_final,
            }),
            IoUnits::Radians,
            IoUnits::Radians,
        ));
    }

    let first_is_tiff = is_tiff(&resolved[0].1);
    let mut grids: Vec<GridSet> = Vec::new();
    let mut hplan: Option<HShiftPlan> = None;

    for (name, bytes) in &resolved {
        if first_is_tiff {
            // Multi-IFD GeoTIFF grids expose a parent/child sub-grid hierarchy
            // (e.g. de_geosn_NTv2_SN.tif: coarse RD83 parent + densified SN_*
            // children). `read_geotiff_hierarchy` returns every georeferenced
            // IFD ordered child-first, so the first-containing-grid-wins scan in
            // `forward_4d`/`inverse_4d` selects the finest sub-grid covering the
            // point, matching PROJ's `HorizontalShiftGrid::gridAt` recursion.
            for gs in read_geotiff_hierarchy(bytes, name)? {
                if hplan.is_none() {
                    hplan = Some(resolve_hshift_plan(&gs).ok_or(ProjError::UnsupportedOperation)?);
                }
                grids.push(gs);
            }
        } else {
            let more = read_ntv2(bytes, name)?;
            if more.is_empty() {
                return Err(ProjError::FileNotFound);
            }
            for gs in &more {
                if gs.bands.len() < 2 {
                    return Err(ProjError::FileNotFound);
                }
            }
            grids.extend(more);
        }
    }

    Ok(TransBuild::new(
        Box::new(HGridShift {
            grids,
            hplan,
            t_epoch,
            t_final,
        }),
        IoUnits::Radians,
        IoUnits::Radians,
    ))
}

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

    // Build a minimal NTv2 file with uniform lat_shift and lon_shift (arc-seconds),
    // covering ll=(0°,0°) to ur=(1°,-1°) (positive-west convention for NTv2).
    fn build_ntv2_with_shifts(lat_shift: f32, lon_shift: f32) -> Vec<u8> {
        let mut buf = Vec::new();
        // File header (11 × 16 = 176 bytes)
        buf.extend_from_slice(b"NUM_OREC");
        buf.extend_from_slice(&11i32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        buf.extend_from_slice(b"NUM_SREC");
        buf.extend_from_slice(&11i32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        buf.extend_from_slice(b"NUM_FILE");
        buf.extend_from_slice(&1u32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        buf.extend_from_slice(b"GS_TYPE ");
        buf.extend_from_slice(b"SECONDS ");
        buf.extend_from_slice(&[0u8; 112]);
        assert_eq!(buf.len(), 176);

        // Sub-grid header (11 × 16 = 176 bytes)
        buf.extend_from_slice(b"SUB_NAME");
        buf.extend_from_slice(b"TESTGRID");
        buf.extend_from_slice(b"PARENT  ");
        buf.extend_from_slice(b"NONE    ");
        buf.extend_from_slice(b"CREATED ");
        buf.extend_from_slice(b"20240101");
        buf.extend_from_slice(b"UPDATED ");
        buf.extend_from_slice(b"20240101");
        buf.extend_from_slice(b"S_LAT   ");
        buf.extend_from_slice(&0.0f64.to_le_bytes());
        buf.extend_from_slice(b"N_LAT   ");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        buf.extend_from_slice(b"E_LONG  ");
        buf.extend_from_slice(&0.0f64.to_le_bytes());
        buf.extend_from_slice(b"W_LONG  ");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        buf.extend_from_slice(b"LAT_INC ");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        buf.extend_from_slice(b"LONG_INC");
        buf.extend_from_slice(&3600.0f64.to_le_bytes());
        buf.extend_from_slice(b"GS_COUNT");
        buf.extend_from_slice(&4i32.to_le_bytes());
        buf.extend_from_slice(&[0u8; 4]);
        assert_eq!(buf.len(), 352);

        // Data: 4 cells (2×2), uniform shift values
        for _ in 0..4 {
            buf.extend_from_slice(&lat_shift.to_le_bytes());
            buf.extend_from_slice(&lon_shift.to_le_bytes());
            buf.extend_from_slice(&0.0f32.to_le_bytes()); // accuracy lat
            buf.extend_from_slice(&0.0f32.to_le_bytes()); // accuracy lon
        }
        buf
    }

    // Test-only registry
    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())
        }
    }

    // Minimal TransParamLookup impl for tests
    #[derive(Default)]
    struct TestParams {
        grids: String,
        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 {
                "t_epoch" => self.t_epoch,
                "t_final" => self.t_final,
                _ => None,
            }
        }
        fn get_dms(&self, _key: &str) -> Option<f64> {
            None
        }
        fn get_int(&self, _key: &str) -> Option<i64> {
            None
        }
        fn get_bool(&self, _key: &str) -> bool {
            false
        }
        fn exists(&self, key: &str) -> bool {
            key == "grids"
                || (key == "t_epoch" && self.t_epoch.is_some())
                || (key == "t_final" && self.t_final.is_some())
        }
    }

    #[test]
    fn test_hgridshift_forward_applies_arcsec_shift() {
        let data = build_ntv2_with_shifts(3600.0, 1800.0); // 1 deg lat, 0.5 deg lon
        let mut map = std::collections::HashMap::new();
        map.insert("test.gsb".to_string(), data);
        let reg = InMemReg(map);
        let tp = TestParams {
            grids: "test.gsb".to_string(),
            ..Default::default()
        };
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&params).unwrap();
        // Point at lon=-0.5°, lat=0.5° (inside the grid: lat 0-1°, lon -1-0°)
        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();
        let v = out.v();
        // Raw lon_shift=1800 arcsec is positive-WEST; read_ntv2 negates it to
        // east-positive (-0.5 deg), so the point moves WEST by 0.5 deg.
        // (Verified convention against PROJ cct +proj=hgridshift.)
        let expected_lon = (-0.5 - 0.5) * DEG_TO_RAD;
        let expected_lat = (0.5 + 1.0) * DEG_TO_RAD;
        assert!(
            (v[0] - expected_lon).abs() < 1e-9,
            "lon: got {}, expected {}",
            v[0],
            expected_lon
        );
        assert!(
            (v[1] - expected_lat).abs() < 1e-9,
            "lat: got {}, expected {}",
            v[1],
            expected_lat
        );
    }

    #[test]
    fn test_hgridshift_round_trip() {
        let data = build_ntv2_with_shifts(1800.0, 900.0);
        let mut map = std::collections::HashMap::new();
        map.insert("test.gsb".to_string(), data);
        let reg = InMemReg(map);
        let tp = TestParams {
            grids: "test.gsb".to_string(),
            ..Default::default()
        };
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&params).unwrap();
        // Point inside the grid: lat 0-1°, lon -1-0°
        let input = Coord::new(-0.7 * DEG_TO_RAD, 0.4 * DEG_TO_RAD, 5.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[0] - vi0[0]).abs() < 1e-9, "lon round-trip error");
        assert!((vi[1] - vi0[1]).abs() < 1e-9, "lat round-trip error");
    }

    #[test]
    fn test_hgridshift_outside_returns_err() {
        let data = build_ntv2_with_shifts(1.0, 1.0);
        let mut map = std::collections::HashMap::new();
        map.insert("test.gsb".to_string(), data);
        let reg = InMemReg(map);
        let tp = TestParams {
            grids: "test.gsb".to_string(),
            ..Default::default()
        };
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&params).unwrap();
        // Point far outside grid (grid covers lat 0-1, lon -1 to 0)
        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_hgridshift_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").unwrap();
        let np = NoParams;
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &np,
            registry: Some(&reg),
        };
        assert_eq!(new(&params).err(), Some(ProjError::MissingArg));
    }

    #[test]
    fn hgridshift_resolves_grid_from_disk_when_registry_empty() {
        // E4: with the in-memory registry empty, the grid must be resolved from
        // a local `PROJ_DATA` directory (registry -> local dirs -> CDN order).
        // `cargo nextest` runs every test in its own process, so mutating
        // PROJ_DATA here cannot race other tests; the prior value is saved and
        // restored regardless.
        use std::io::Write;
        let ns = std::time::SystemTime::now()
            .duration_since(std::time::UNIX_EPOCH)
            .map(|d| d.as_nanos())
            .unwrap_or(0);
        let dir = std::env::temp_dir().join(format!("oxiproj_hgridshift_disk_{ns}"));
        std::fs::create_dir_all(&dir).expect("create temp PROJ_DATA dir");
        // +1 deg latitude shift, 0 deg longitude shift, over lon[-1,0], lat[0,1].
        let grid = build_ntv2_with_shifts(3600.0, 0.0);
        {
            let mut f = std::fs::File::create(dir.join("disk_shift.gsb")).expect("create grid");
            f.write_all(&grid).expect("write grid");
        }

        let saved = std::env::var_os("PROJ_DATA");
        std::env::set_var("PROJ_DATA", &dir);

        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let tp = TestParams {
            grids: "disk_shift.gsb".to_string(),
            ..Default::default()
        };
        // Registry is present but empty: resolution must fall through to disk.
        let empty = InMemReg(std::collections::HashMap::new());
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&empty),
        };
        let built = new(&params);

        // Restore env and clean up before asserting so a failure never leaks.
        match saved {
            Some(v) => std::env::set_var("PROJ_DATA", v),
            None => std::env::remove_var("PROJ_DATA"),
        }
        let _ = std::fs::remove_dir_all(&dir);

        let tb = built.expect("hgridshift must resolve its grid from PROJ_DATA");
        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().v();
        assert!(
            (out[0] - (-0.5 * DEG_TO_RAD)).abs() < 1e-9,
            "lon unchanged, got {}",
            out[0]
        );
        assert!(
            (out[1] - 1.5 * DEG_TO_RAD).abs() < 1e-9,
            "lat shifted +1 deg, got {}",
            out[1]
        );
    }

    #[test]
    fn time_gate_skips_points_outside_bracket() {
        // +t_epoch=2000, +t_final=2010: only points with t < 2000 are shifted.
        let data = build_ntv2_with_shifts(3600.0, 0.0); // +1 deg lat, 0 lon
        let mut map = std::collections::HashMap::new();
        map.insert("test.gsb".to_string(), data);
        let reg = InMemReg(map);
        let tp = TestParams {
            grids: "test.gsb".to_string(),
            t_epoch: Some(2000.0),
            t_final: Some(2010.0),
        };
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&params).unwrap();
        // t = 1995 < t_epoch => shift applied (+1 deg lat).
        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().v();
        assert!(
            (out_in[1] - 1.5 * DEG_TO_RAD).abs() < 1e-9,
            "in-bracket lat shifted +1 deg, got {}",
            out_in[1].to_degrees()
        );
        // 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().v();
        assert!(
            (out_out[1] - 0.5 * DEG_TO_RAD).abs() < 1e-12,
            "out-of-bracket lat unchanged, got {}",
            out_out[1].to_degrees()
        );
    }

    #[test]
    fn comma_split_optional_prefix_falls_back() {
        // `+grids=@missing.gsb,test.gsb`: the optional first entry is skipped and
        // the second grid is used.
        let data = build_ntv2_with_shifts(3600.0, 0.0);
        let mut map = std::collections::HashMap::new();
        map.insert("test.gsb".to_string(), data);
        let reg = InMemReg(map);
        let tp = TestParams {
            grids: "@missing.gsb,test.gsb".to_string(),
            ..Default::default()
        };
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&params).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().v();
        assert!(
            (out[1] - 1.5 * DEG_TO_RAD).abs() < 1e-9,
            "fell back to second grid (+1 deg lat), got {}",
            out[1].to_degrees()
        );
    }

    #[test]
    fn required_missing_grid_errors() {
        // A non-`@` (required) missing grid aborts the whole list.
        let data = build_ntv2_with_shifts(3600.0, 0.0);
        let mut map = std::collections::HashMap::new();
        map.insert("test.gsb".to_string(), data);
        let reg = InMemReg(map);
        let tp = TestParams {
            grids: "missing.gsb,test.gsb".to_string(),
            ..Default::default()
        };
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let params = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        assert_eq!(new(&params).err(), Some(ProjError::FileNotFound));
    }
}