oxiproj-transformations 0.1.2

#![forbid(unsafe_code)]
//! Hermetic regression for PROJ's multi-grid per-point selection on the
//! `+proj=hgridshift` path (`pj_hgrid_apply` first-containing-grid-wins scan).
//!
//! Two synthetic NTv2 grids with overlapping-but-different extents and *distinct*
//! shift values are listed together in `+grids=A,B`. A point that falls in A only
//! must get A's shift; a point in B only must get B's shift; a point in the
//! overlap must get the shift of whichever grid is listed *first* (the
//! first-containing-grid-wins rule, and its inverse when the order is swapped).
//! No `PROJ_DATA`, no network, no Homebrew — the grids are built in memory.

use oxiproj_core::{Coord, Ellipsoid, ProjError, DEG_TO_RAD};
use oxiproj_transformations::{build, GridRegistry, TransParamLookup, TransParams};

// ---------------------------------------------------------------------------
// Synthetic NTv2 grid builder (uniform shift over a configurable extent).
// ---------------------------------------------------------------------------

/// Build a minimal little-endian NTv2 (.gsb) grid with a single sub-grid holding
/// a uniform `(lat_shift, lon_shift)` in arc-seconds over the geographic box
/// `lon ∈ [west_lon_deg, east_lon_deg]`, `lat ∈ [south_lat_deg, north_lat_deg]`.
/// NTv2 longitude is positive-WEST (`read_ntv2` re-normalizes to east-positive).
fn build_ntv2(
    south_lat_deg: f64,
    north_lat_deg: f64,
    west_lon_deg: f64,
    east_lon_deg: f64,
    lat_shift_arcsec: f32,
    lon_shift_arcsec: f32,
) -> Vec<u8> {
    let s_lat = south_lat_deg * 3600.0;
    let n_lat = north_lat_deg * 3600.0;
    let e_long = -east_lon_deg * 3600.0;
    let w_long = -west_lon_deg * 3600.0;
    let n_rows = 2usize;
    let n_cols = 2usize;
    let lat_inc = (n_lat - s_lat) / (n_rows - 1) as f64;
    let long_inc = (w_long - e_long) / (n_cols - 1) as f64;
    let gs_count = (n_rows * n_cols) as i32;

    let mut buf = Vec::new();
    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);

    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(&s_lat.to_le_bytes());
    buf.extend_from_slice(b"N_LAT   ");
    buf.extend_from_slice(&n_lat.to_le_bytes());
    buf.extend_from_slice(b"E_LONG  ");
    buf.extend_from_slice(&e_long.to_le_bytes());
    buf.extend_from_slice(b"W_LONG  ");
    buf.extend_from_slice(&w_long.to_le_bytes());
    buf.extend_from_slice(b"LAT_INC ");
    buf.extend_from_slice(&lat_inc.to_le_bytes());
    buf.extend_from_slice(b"LONG_INC");
    buf.extend_from_slice(&long_inc.to_le_bytes());
    buf.extend_from_slice(b"GS_COUNT");
    buf.extend_from_slice(&gs_count.to_le_bytes());
    buf.extend_from_slice(&[0u8; 4]);
    assert_eq!(buf.len(), 352);

    for _ in 0..gs_count {
        buf.extend_from_slice(&lat_shift_arcsec.to_le_bytes());
        buf.extend_from_slice(&lon_shift_arcsec.to_le_bytes());
        buf.extend_from_slice(&0.0f32.to_le_bytes());
        buf.extend_from_slice(&0.0f32.to_le_bytes());
    }
    buf
}

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

struct GridsParam(String);
impl TransParamLookup for GridsParam {
    fn get_str(&self, key: &str) -> Option<&str> {
        (key == "grids").then_some(self.0.as_str())
    }
    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"
    }
}

// Grid A: lat ∈ [0,4], lon ∈ [-4,0], +1° latitude shift (3600 arc-second).
// Grid B: lat ∈ [2,6], lon ∈ [-6,-2], +2° latitude shift (7200 arc-second).
// Overlap: lat ∈ [2,4], lon ∈ [-4,-2].
fn grid_a() -> Vec<u8> {
    build_ntv2(0.0, 4.0, -4.0, 0.0, 3600.0, 0.0)
}
fn grid_b() -> Vec<u8> {
    build_ntv2(2.0, 6.0, -6.0, -2.0, 7200.0, 0.0)
}

fn registry() -> InMemReg {
    let mut map = std::collections::HashMap::new();
    map.insert("grid_a.gsb".to_string(), grid_a());
    map.insert("grid_b.gsb".to_string(), grid_b());
    InMemReg(map)
}

/// Forward the geographic point `(lon_deg, lat_deg)` through `+grids=<spec>` and
/// return the output latitude in degrees.
fn forward_lat_deg(grids_spec: &str, lon_deg: f64, lat_deg: f64) -> Result<f64, ProjError> {
    let reg = registry();
    let ell = Ellipsoid::named("WGS84").expect("WGS84");
    let tp = GridsParam(grids_spec.to_string());
    let params = TransParams {
        ellipsoid: &ell,
        params: &tp,
        registry: Some(&reg),
    };
    let tb = build("hgridshift", &params)?;
    let input = Coord::new(lon_deg * DEG_TO_RAD, lat_deg * DEG_TO_RAD, 0.0, 0.0);
    let out = tb.operation.forward_4d(input)?;
    Ok(out.v()[1].to_degrees())
}

#[test]
fn point_in_grid_a_only_gets_a_shift() {
    // (lon=-1, lat=1): inside A (lat[0,4], lon[-4,0]), outside B (lat[2,6]).
    let lat = forward_lat_deg("grid_a.gsb,grid_b.gsb", -1.0, 1.0).expect("forward");
    assert!(
        (lat - 2.0).abs() < 1e-9,
        "A-only point must get A's +1°: got {lat}"
    );
}

#[test]
fn point_in_grid_b_only_gets_b_shift() {
    // (lon=-5, lat=5): inside B (lat[2,6], lon[-6,-2]), outside A (lon[-4,0]).
    // A is listed first but does not contain the point, so B is selected.
    let lat = forward_lat_deg("grid_a.gsb,grid_b.gsb", -5.0, 5.0).expect("forward");
    assert!(
        (lat - 7.0).abs() < 1e-9,
        "B-only point must get B's +2°: got {lat}"
    );
}

#[test]
fn overlap_point_takes_first_listed_grid() {
    // (lon=-3, lat=3): inside BOTH. With A listed first, A wins (+1°).
    let lat_ab = forward_lat_deg("grid_a.gsb,grid_b.gsb", -3.0, 3.0).expect("forward A,B");
    assert!(
        (lat_ab - 4.0).abs() < 1e-9,
        "overlap with A first must get A's +1°: got {lat_ab}"
    );
    // Swapping the order selects B (+2°) for the same point, proving selection
    // is by list order among containing grids, not a fixed grid.
    let lat_ba = forward_lat_deg("grid_b.gsb,grid_a.gsb", -3.0, 3.0).expect("forward B,A");
    assert!(
        (lat_ba - 5.0).abs() < 1e-9,
        "overlap with B first must get B's +2°: got {lat_ba}"
    );
}

#[test]
fn point_outside_all_grids_errors() {
    // (lon=20, lat=20): outside both A and B.
    let err = forward_lat_deg("grid_a.gsb,grid_b.gsb", 20.0, 20.0)
        .expect_err("outside all grids must error");
    assert_eq!(err, ProjError::OutsideGrid);
}