oxiproj_transformations/transformations/

gridshift.rs

#![forbid(unsafe_code)]
//! Generic grid shift (`gridshift`) — port of PROJ 9.8.0 generic gridshift.
//! Detects the grid format from raw bytes (NTv2, GTX, GeoTIFF) and applies
//! the appropriate horizontal, vertical, or 3-D shift. Acts as a dispatcher.

use crate::{TransBuild, TransParams};
use oxiproj_core::{Coord, IoUnits, Operation, ProjError, ProjResult};
use oxiproj_grids::{
    read_geotiff, read_geotiff_gdal_metadata, read_geotiff_hierarchy, read_gtx, read_ntv2,
    sample_grid, BandPositive, BandRole, BandUnit, GridBand, GridSet,
};
use std::f64::consts::PI;

const ARCSEC_TO_RAD: f64 = PI / 648_000.0;
const DEG_TO_RAD: f64 = PI / 180.0;
const MAX_ITER: usize = 10;
const ITER_TOL: f64 = 1e-12;
/// Relative tolerance for snapping a query that lands a hair outside the grid
/// back onto the edge cell, matching PROJ's `REL_TOLERANCE_HGRIDSHIFT`.
const REL_TOLERANCE_HGRIDSHIFT: f64 = 1e-5;

/// Interpolation kernel selected via `+interpolation=` on the generic
/// `gridshift`, mirroring PROJ `gridshift.cpp` (`bilinear` | `biquadratic`).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(crate) enum Interpolation {
    /// Standard 2×2 bilinear interpolation (PROJ default).
    Bilinear,
    /// 3×3-window quadratic interpolation (NOAA `qterp`), used e.g. for the
    /// NADCON5 GeoTIFF grids to reproduce the NCAT reference transformer.
    Biquadratic,
}

/// Return `true` if `data` begins with a TIFF magic signature (little-endian
/// `II\x2a\x00` or big-endian `MM\x00\x2a`), i.e. a GeoTIFF grid rather than an
/// NTv2 `.gsb` / GTX binary. Shared by `vgridshift`/`hgridshift`.
pub(crate) fn is_tiff(data: &[u8]) -> bool {
    data.starts_with(b"II\x2a\x00") || data.starts_with(b"MM\x00\x2a")
}

/// Split a `+grids=` value on commas and resolve each entry's raw bytes,
/// honoring PROJ's per-entry `@` "optional" prefix (`src/grids.cpp`
/// `pj_hgrid_init`/`pj_vgrid_init`/`pj_generic_grid_init`).
///
/// Each comma-separated entry may carry a leading `@`, meaning "may fail, skip
/// silently"; a non-`@` entry that cannot be resolved aborts the whole list with
/// its resolution error (PROJ's "could not find required grid(s)"). The returned
/// list preserves input order, giving the same first-containing-grid-wins
/// fallback semantics PROJ applies when sampling. An all-optional list whose
/// entries all fail resolves to an empty list (a pass-through transform), exactly
/// as PROJ leaves `Q->grids` empty in that case.
pub(crate) fn resolve_grid_list(
    registry: Option<&dyn crate::GridRegistry>,
    csv: &str,
) -> ProjResult<Vec<(String, Vec<u8>)>> {
    let mut out = Vec::new();
    for raw in csv.split(',') {
        let entry = raw.trim();
        if entry.is_empty() {
            continue;
        }
        let (name, can_fail) = match entry.strip_prefix('@') {
            Some(rest) => (rest.trim(), true),
            None => (entry, false),
        };
        if name.is_empty() {
            continue;
        }
        match resolve_grid_bytes(registry, name) {
            Ok(bytes) => out.push((name.to_string(), bytes)),
            Err(e) => {
                if !can_fail {
                    return Err(e);
                }
            }
        }
    }
    Ok(out)
}

/// Parse the optional `+t_epoch` / `+t_final` time-restriction parameters shared
/// by `hgridshift` and `vgridshift` (PROJ `hgridshift.cpp` / `vgridshift.cpp`).
///
/// `+t_final=now` is resolved to the current decimal year exactly as PROJ does
/// (`1900 + tm_year + tm_yday/365`). Absent parameters yield `0.0`, which
/// disables gating (the shift is always applied).
pub(crate) fn parse_time_gate(p: &TransParams) -> (f64, f64) {
    let t_epoch = p.params.get_f64("t_epoch").unwrap_or(0.0);
    let t_final = match p.params.get_f64("t_final") {
        Some(v) if v != 0.0 => v,
        _ => {
            // A number wasn't passed (or it was 0); honor the "now" keyword.
            if p.params.get_str("t_final") == Some("now") {
                decimal_year_now()
            } else {
                0.0
            }
        }
    };
    (t_epoch, t_final)
}

/// Return `true` when a grid shift should be applied to a point observed at
/// decimal year `t`, per PROJ's `pj_*gridshift_forward_4d` gating: apply
/// unconditionally when either bound is 0, otherwise only when
/// `t < t_epoch && t_final > t_epoch`.
pub(crate) fn time_gate_applies(t_epoch: f64, t_final: f64, t: f64) -> bool {
    if t_final == 0.0 || t_epoch == 0.0 {
        return true;
    }
    t < t_epoch && t_final > t_epoch
}

/// Current time as a decimal year, mirroring PROJ's
/// `1900.0 + date->tm_year + date->tm_yday / 365.0` (UTC here; the sub-day
/// difference from local time is negligible for a year fraction). Used only for
/// `+t_final=now`.
fn decimal_year_now() -> f64 {
    let secs = std::time::SystemTime::now()
        .duration_since(std::time::UNIX_EPOCH)
        .map(|d| d.as_secs() as i64)
        .unwrap_or(0);
    let days = secs.div_euclid(86_400);
    let (year, month, day) = civil_from_days(days);
    let doy = day_of_year(year, month, day);
    year as f64 + doy as f64 / 365.0
}

/// Convert a count of days since 1970-01-01 to a `(year, month, day)` civil date
/// (Howard Hinnant's `civil_from_days`, valid across the full proleptic
/// Gregorian range).
fn civil_from_days(days: i64) -> (i64, u32, u32) {
    let z = days + 719_468;
    let era = if z >= 0 { z } else { z - 146_096 } / 146_097;
    let doe = z - era * 146_097;
    let yoe = (doe - doe / 1460 + doe / 36_524 - doe / 146_096) / 365;
    let y = yoe + era * 400;
    let doy = doe - (365 * yoe + yoe / 4 - yoe / 100);
    let mp = (5 * doy + 2) / 153;
    let d = (doy - (153 * mp + 2) / 5 + 1) as u32;
    let m = if mp < 10 { mp + 3 } else { mp - 9 } as u32;
    let year = if m <= 2 { y + 1 } else { y };
    (year, m, d)
}

/// Zero-based day of year (matching C's `tm_yday`) for a civil date.
fn day_of_year(year: i64, month: u32, day: u32) -> u32 {
    const CUM: [u32; 12] = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334];
    let leap = (year % 4 == 0 && year % 100 != 0) || (year % 400 == 0);
    let mut doy = CUM[(month.clamp(1, 12) - 1) as usize] + day.saturating_sub(1);
    if leap && month > 2 {
        doy += 1;
    }
    doy
}

/// Pick the index of the vertical-correction band in a GeoTIFF vertical grid,
/// mirroring PROJ's `GTiffVGridShiftSet::open` band selection: prefer a band
/// whose GDAL_METADATA `DESCRIPTION` marks it as a geoid/vertical offset,
/// otherwise fall back to band 0.
pub(crate) fn resolve_vertical_band(gs: &GridSet) -> Option<usize> {
    if gs.bands.is_empty() {
        return None;
    }
    for (i, band) in gs.bands.iter().enumerate() {
        if matches!(
            band.semantics.role,
            BandRole::GeoidUndulation | BandRole::VerticalOffset
        ) {
            return Some(i);
        }
    }
    // No described vertical channel: default to the first band (single-band
    // geoid grids such as us_nga_egm96_15.tif carry no per-band DESCRIPTION).
    Some(0)
}

/// Reduce a multi-band GeoTIFF vertical grid to a single-band [`GridSet`]
/// holding only the resolved geoid/vertical-offset band, so the shared
/// `ShiftKind::Vertical` sampler (which reads band 0) applies the correct
/// channel in metres.
pub(crate) fn vertical_gridset(mut gs: GridSet) -> ProjResult<GridSet> {
    let idx = resolve_vertical_band(&gs).ok_or(ProjError::FileNotFound)?;
    if idx != 0 {
        gs.bands.swap(0, idx);
    }
    gs.bands.truncate(1);
    Ok(gs)
}

/// Resolved recipe for applying a 2-band GeoTIFF horizontal-offset grid,
/// mirroring PROJ's `GTiffHGridShiftSet::open` (`src/grids.cpp`): which band is
/// the latitude / longitude offset, the arc-second/degree/radian → radian
/// conversion factor, and whether positive longitude values point east.
#[derive(Debug, Clone, Copy)]
pub(crate) struct HShiftPlan {
    idx_lat: usize,
    idx_long: usize,
    conv_factor_to_rad: f64,
    positive_east: bool,
}

/// Convert a declared [`BandUnit`] to its radian conversion factor, or `None`
/// for a non-angular unit (which cannot be a horizontal offset).
fn unit_conv_factor(unit: BandUnit) -> Option<f64> {
    match unit {
        // PROJ's default when UNITTYPE is absent is arc-second (see the
        // `convFactorToRadian = ARC_SECOND_TO_RADIAN` initializer).
        BandUnit::ArcSecond | BandUnit::Unknown => Some(ARCSEC_TO_RAD),
        BandUnit::Degree => Some(DEG_TO_RAD),
        BandUnit::Radian => Some(1.0),
        BandUnit::Metre => None,
    }
}

/// Resolve the [`HShiftPlan`] for a GeoTIFF grid set, replicating PROJ's band
/// identification and its NTv2-inspired defaults (band 0 = latitude offset,
/// band 1 = longitude offset, arc-second, positive east) when the
/// GDAL_METADATA descriptions are absent.
///
/// Returns `None` when the grid cannot be interpreted as a horizontal-offset
/// grid (fewer than 2 bands; only one of lat/long described; a described but
/// non-horizontal grid; or a non-angular unit).
pub(crate) fn resolve_hshift_plan(gs: &GridSet) -> Option<HShiftPlan> {
    let nbands = gs.bands.len();
    if nbands < 2 {
        return None;
    }

    // PROJ defaults, inspired from NTv2.
    let mut idx_lat = 0usize;
    let mut idx_long = 1usize;
    let mut found_lat = false;
    let mut found_long = false;
    let mut found_any_desc = false;

    for (i, band) in gs.bands.iter().enumerate() {
        match band.semantics.role {
            BandRole::LatitudeOffset => {
                idx_lat = i;
                found_lat = true;
                found_any_desc = true;
            }
            BandRole::LongitudeOffset => {
                idx_long = i;
                found_long = true;
                found_any_desc = true;
            }
            BandRole::Unknown => {}
            // A described but non-horizontal role (geoid/vertical) still counts
            // as "descriptions present".
            BandRole::GeoidUndulation | BandRole::VerticalOffset => {
                found_any_desc = true;
            }
        }
    }

    // Descriptions present but neither offset channel found: not horizontal.
    if found_any_desc && !found_lat && !found_long {
        return None;
    }
    // Exactly one of the two channels described: PROJ treats this as an error.
    if found_lat != found_long {
        return None;
    }
    if idx_lat >= nbands || idx_long >= nbands {
        return None;
    }

    // Unit (PROJ requires lat and long to share it); default arc-second.
    let conv_factor_to_rad = unit_conv_factor(gs.bands[idx_lat].semantics.unit)?;
    // positive_value on the longitude band; default east.
    let positive_east = gs.bands[idx_long].semantics.positive != BandPositive::West;

    Some(HShiftPlan {
        idx_lat,
        idx_long,
        conv_factor_to_rad,
        positive_east,
    })
}

/// Apply an [`HShiftPlan`] to interpolated band samples, returning the
/// longitude and latitude offsets in radians (`(dlon, dlat)`), matching PROJ's
/// `GTiffHGrid::valueAt` (arc-second → radian, then negate longitude when
/// positive values point west).
pub(crate) fn geotiff_horizontal_offset(plan: &HShiftPlan, shifts: &[f64]) -> (f64, f64) {
    let raw_lat = shifts.get(plan.idx_lat).copied().unwrap_or(0.0);
    let raw_long = shifts.get(plan.idx_long).copied().unwrap_or(0.0);
    let dlat = raw_lat * plan.conv_factor_to_rad;
    let mut dlon = raw_long * plan.conv_factor_to_rad;
    if !plan.positive_east {
        dlon = -dlon;
    }
    (dlon, dlat)
}

/// Build a semantics-aware horizontal GeoTIFF shift operation from an already
/// decoded [`GridSet`]. Used by the `gridshift` unit tests to exercise the
/// horizontal-offset application path directly from an in-memory [`GridSet`].
/// Returns an error if the grid cannot be interpreted as a horizontal-offset
/// grid.
#[cfg(test)]
pub(crate) fn build_horizontal_geotiff(gs: GridSet) -> ProjResult<TransBuild> {
    let plan = resolve_hshift_plan(&gs).ok_or(ProjError::UnsupportedOperation)?;
    Ok(TransBuild::new(
        Box::new(GridShift {
            grids: vec![gs],
            kind: ShiftKind::HorizontalGeoTiff,
            hplan: Some(plan),
            idx_z: None,
            interpolation: Interpolation::Bilinear,
            no_z_transform: false,
        }),
        IoUnits::Radians,
        IoUnits::Radians,
    ))
}

/// Port of the NOAA `qterp` Fortran routine (a parabola through three equally
/// spaced samples). `x` is the fractional position in `[0, 2]`; `f0`, `f1`, `f2`
/// are the samples at 0, 1, 2. Matches PROJ `gridshift.cpp`'s
/// `quadraticInterpol` lambda exactly.
fn quadratic_interpol(x: f64, f0: f64, f1: f64, f2: f64) -> f64 {
    let df0 = f1 - f0;
    let df1 = f2 - f1;
    let d2f0 = df1 - df0;
    f0 + x * df0 + 0.5 * x * (x - 1.0) * d2f0
}

/// Biquadratic (3×3-window quadratic) interpolation of every band of `gs` at the
/// geographic position `(lat_deg, lon_deg)`, returning one value per band in the
/// bands' stored units (arc-second / metre — the same raw values `sample_grid`
/// returns), or `None` when the point is outside the grid or hits nodata.
///
/// This is a faithful port of PROJ `gridshift.cpp`'s `grid_interpolate`
/// biquadratic branch (its `qterp` window-shift + `quadraticInterpol`), operating
/// in this crate's north-first raster convention. Grids narrower/shorter than
/// 3 cells fall back to bilinear, exactly as PROJ does
/// (`grid->width() < 3 || grid->height() < 3`).
fn sample_grid_biquadratic(gs: &GridSet, lat_deg: f64, lon_deg: f64) -> Option<Vec<f64>> {
    if gs.bands.is_empty() {
        return None;
    }
    let mut results = Vec::with_capacity(gs.bands.len());
    for band in &gs.bands {
        results.push(band_biquadratic(band, lat_deg, lon_deg)?);
    }
    Some(results)
}

/// Biquadratic interpolation of a single [`GridBand`]; see
/// [`sample_grid_biquadratic`].
fn band_biquadratic(band: &GridBand, lat_deg: f64, lon_deg: f64) -> Option<f64> {
    let ext = &band.extent;
    let width = ext.cols() as i64;
    let height = ext.rows() as i64;
    // PROJ falls back to bilinear for grids too small for a 3×3 window.
    if width < 3 || height < 3 {
        return band_bilinear_fallback(band, lat_deg, lon_deg);
    }

    // Longitude normalization (single ±360° adjustment, as PROJ's normalizeX).
    let mut lon = lon_deg;
    let eps = (ext.lon_inc + ext.lat_inc) * REL_TOLERANCE_HGRIDSHIFT;
    if lon < ext.ll_lon - eps {
        lon += 360.0;
    } else if lon > ext.ur_lon + eps {
        lon -= 360.0;
    }

    // West-based column index and south-based row index (PROJ's grid frame).
    let x = (lon - ext.ll_lon) / ext.lon_inc;
    let y = (lat_deg - ext.ll_lat) / ext.lat_inc;
    let mut ix = x.floor() as i64;
    let mut iy = y.floor() as i64;
    let mut fx = x - ix as f64;
    let mut fy = y - iy as f64;

    // Boundary snap (PROJ `grid_interpolate`).
    if ix < 0 {
        if ix == -1 && fx > 1.0 - 10.0 * REL_TOLERANCE_HGRIDSHIFT {
            ix += 1;
            fx = 0.0;
        } else {
            return None;
        }
    } else if ix + 1 >= width {
        if ix + 1 == width && fx < 10.0 * REL_TOLERANCE_HGRIDSHIFT {
            ix -= 1;
            fx = 1.0;
        } else {
            return None;
        }
    }
    if iy < 0 {
        if iy == -1 && fy > 1.0 - 10.0 * REL_TOLERANCE_HGRIDSHIFT {
            iy += 1;
            fy = 0.0;
        } else {
            return None;
        }
    } else if iy + 1 >= height {
        if iy + 1 == height && fy < 10.0 * REL_TOLERANCE_HGRIDSHIFT {
            iy -= 1;
            fy = 1.0;
        } else {
            return None;
        }
    }

    // Shift the 3×3 window depending on which half-pixel we are in.
    if (fx <= 0.5 && ix > 0) || (ix + 2 == width) {
        ix -= 1;
        fx += 1.0;
    }
    if (fy <= 0.5 && iy > 0) || (iy + 2 == height) {
        iy -= 1;
        fy += 1.0;
    }
    if ix < 0 || iy < 0 || ix + 2 >= width || iy + 2 >= height {
        return None;
    }

    // Read the 3×3 window and qterp along x per south-row, then along y.
    // Grid rows are stored north-first, so south-based row j maps to
    // north_row = (height - 1) - (iy + j).
    let mut row_vals = [0.0f64; 3];
    for (j, slot) in row_vals.iter_mut().enumerate() {
        let north_row = (height - 1 - (iy + j as i64)) as usize;
        let mut cols = [0.0f64; 3];
        for (i, c) in cols.iter_mut().enumerate() {
            let v = band.get(north_row, (ix + i as i64) as usize)? as f64;
            if !v.is_finite() {
                return None;
            }
            *c = v;
        }
        *slot = quadratic_interpol(fx, cols[0], cols[1], cols[2]);
    }
    Some(quadratic_interpol(
        fy,
        row_vals[0],
        row_vals[1],
        row_vals[2],
    ))
}

/// Single-band bilinear fallback used by [`band_biquadratic`] for grids too
/// small for a 3×3 window. Reuses the shared [`sample_grid`] by wrapping the one
/// band in a temporary [`GridSet`] is avoided; instead we sample directly.
fn band_bilinear_fallback(band: &GridBand, lat_deg: f64, lon_deg: f64) -> Option<f64> {
    let tmp = GridSet {
        bands: vec![band.clone()],
        source: oxiproj_grids::GridSource::Memory,
    };
    sample_grid(&tmp, lat_deg, lon_deg).and_then(|v| v.first().copied())
}

/// Sample `gs` with the selected [`Interpolation`] kernel, returning one raw
/// value per band (arc-second / degree / metre, as stored).
fn sample_with(
    interp: Interpolation,
    gs: &GridSet,
    lat_deg: f64,
    lon_deg: f64,
) -> Option<Vec<f64>> {
    match interp {
        Interpolation::Bilinear => sample_grid(gs, lat_deg, lon_deg),
        Interpolation::Biquadratic => sample_grid_biquadratic(gs, lat_deg, lon_deg),
    }
}

/// Resolve the raw bytes of a grid file for a grid-based transformation.
///
/// This is the single resolution path shared by every grid-consuming transform
/// (`gridshift`, `hgridshift`, `vgridshift`, `xyzgridshift`, and — via the
/// engine's `+nadgrids=`/`+geoidgrids=` synthesis — the classic datum/geoid
/// shifts). It mirrors PROJ's `pj_open_lib` resource-resolution order:
///
/// 1. **In-memory registry** — bytes registered on the engine `Context` (via
///    `register_grid`) or any other [`GridRegistry`](crate::GridRegistry) the
///    caller passes. Registered bytes always win, keeping the in-memory
///    registration API authoritative (existing callers and tests rely on it).
/// 2. **Local resource directories** — the `PROJ_DATA` (modern) and `PROJ_LIB`
///    (legacy) directory lists, each of which may hold several
///    platform-path-separator-delimited directories, followed by the unified
///    downloadable-resource cache directory shared across the workspace
///    ([`oxiproj_grids::GridDiskCache::default_dir`]).
/// 3. **PROJ CDN** — only when this crate's off-by-default `network` feature is
///    enabled, the grid is fetched (and cached on disk for next time) via
///    [`oxiproj_grids::GridResolver`].
///
/// Returns [`ProjError::FileNotFound`] when the grid cannot be located in any
/// layer, preserving the previous "grid missing" error for callers.
pub(crate) fn resolve_grid_bytes(
    registry: Option<&dyn crate::GridRegistry>,
    name: &str,
) -> ProjResult<Vec<u8>> {
    // 1. In-memory registry (explicit registered bytes).
    if let Some(reg) = registry {
        if let Some(bytes) = reg.get_grid(name) {
            return Ok(bytes.to_vec());
        }
    }
    // 2. Local PROJ_DATA / PROJ_LIB directories and the unified disk cache.
    if let Some(bytes) = read_from_resource_dirs(name) {
        return Ok(bytes);
    }
    // 3. Network CDN, only with the `network` feature enabled.
    #[cfg(feature = "network")]
    {
        if let Some(bytes) = fetch_from_cdn(name) {
            return Ok(bytes);
        }
    }
    Err(ProjError::FileNotFound)
}

/// Read a grid named `name` from the first local resource directory that
/// contains it, or `None` if no directory holds it (or it cannot be read).
///
/// The searched directories — the `PROJ_DATA` and `PROJ_LIB` path lists, then
/// the unified disk-cache directory — come from the single shared source of
/// truth [`oxiproj_grids::resource_dirs`], the same list the factory's
/// no-network availability check consults, so "found locally" means the same
/// thing on both sides. No absolute path is ever hard-coded.
fn read_from_resource_dirs(name: &str) -> Option<Vec<u8>> {
    for dir in oxiproj_grids::resource_dirs() {
        let path = dir.join(name);
        if path.is_file() {
            if let Ok(bytes) = std::fs::read(&path) {
                return Some(bytes);
            }
        }
    }
    None
}

/// Fetch a grid from the PROJ CDN via [`oxiproj_grids::GridResolver`], caching
/// it in the unified disk cache for subsequent lookups. Compiled only with the
/// `network` feature; a clean miss or any error yields `None` so the caller
/// falls through to [`ProjError::FileNotFound`].
///
/// Building a grid op is eager: it resolves the grid bytes up front. To keep
/// operation selection offline-safe and deterministic by default — mirroring
/// PROJ's network-off default — the CDN is consulted ONLY when the process-wide
/// build-time network policy is explicitly enabled
/// ([`oxiproj_grids::network_build_allowed`], turned on by the factory solely
/// under `GridAvailability::AllowNetwork`). When it is off (the default) a
/// missing grid resolves to a clean [`ProjError::FileNotFound`] instead of
/// blocking on a network fetch.
#[cfg(feature = "network")]
fn fetch_from_cdn(name: &str) -> Option<Vec<u8>> {
    if !oxiproj_grids::network_build_allowed() {
        return None;
    }
    let dir = oxiproj_grids::GridDiskCache::default_dir().unwrap_or_else(std::env::temp_dir);
    let mut resolver = oxiproj_grids::GridResolver::new(dir).ok()?;
    if resolver.ensure(name).ok()? {
        return resolver.get_loaded(name).map(<[u8]>::to_vec);
    }
    None
}

/// Determines how the interpolated band values are interpreted as shifts.
#[derive(Debug)]
enum ShiftKind {
    /// 2-band NTv2: band0=lat-shift (arcsec), band1=lon-shift (arcsec).
    HorizontalArcSec,
    /// 1-band GTX / GeoTIFF vertical: band0=vertical shift (meters).
    Vertical,
    /// 2-band GeoTIFF horizontal offset applied via a semantics-aware
    /// [`HShiftPlan`] (role/unit/positive from GDAL_METADATA), matching PROJ.
    HorizontalGeoTiff,
    /// GeoTIFF `TYPE=GEOGRAPHIC_3D_OFFSET`: a semantics-aware horizontal offset
    /// ([`HShiftPlan`]) plus a vertical (`geoid_undulation`/`vertical_offset`/…)
    /// band in metres, applied together in a single pass (PROJ `gridshift.cpp`).
    Geographic3DOffset,
    /// 3-band: band0=lon-shift (deg), band1=lat-shift (deg), band2=dz (m).
    Xyz,
}

#[derive(Debug)]
struct GridShift {
    grids: Vec<GridSet>,
    kind: ShiftKind,
    /// Resolved horizontal-offset recipe; `Some` for
    /// [`ShiftKind::HorizontalGeoTiff`] and [`ShiftKind::Geographic3DOffset`].
    hplan: Option<HShiftPlan>,
    /// Index of the vertical (metre) band; `Some` only for
    /// [`ShiftKind::Geographic3DOffset`].
    idx_z: Option<usize>,
    /// Interpolation kernel (`+interpolation=bilinear|biquadratic`).
    interpolation: Interpolation,
    /// When set (`+no_z_transform`), the vertical component is left untouched.
    no_z_transform: bool,
}

impl Operation for GridShift {
    fn forward_4d(&self, c: Coord) -> ProjResult<Coord> {
        let v = c.v();
        // Empty grid list (all-optional grids absent): pass through unchanged,
        // matching PROJ's "only try the gridshift if at least one grid is
        // loaded" behavior.
        if self.grids.is_empty() {
            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_with(self.interpolation, gs, lat_deg, lon_deg) {
                return match self.kind {
                    ShiftKind::HorizontalArcSec => Ok(Coord::new(
                        v[0] + shifts[1] * ARCSEC_TO_RAD,
                        v[1] + shifts[0] * ARCSEC_TO_RAD,
                        v[2],
                        v[3],
                    )),
                    // The generic `gridshift` operator ADDS the vertical shift in
                    // the forward direction (PROJ `gridshift.cpp`
                    // `grid_apply_internal`: `out.z += shift.z`). This is the
                    // OPPOSITE of the dedicated `vgridshift` operator (which
                    // defaults `forward_multiplier = -1`); the two operators must
                    // not be conflated.
                    ShiftKind::Vertical => {
                        let z = if self.no_z_transform {
                            v[2]
                        } else {
                            v[2] + shifts[0]
                        };
                        Ok(Coord::new(v[0], v[1], z, v[3]))
                    }
                    ShiftKind::HorizontalGeoTiff => {
                        let plan = self.hplan.as_ref().ok_or(ProjError::UnsupportedOperation)?;
                        let (dlon, dlat) = geotiff_horizontal_offset(plan, &shifts);
                        Ok(Coord::new(v[0] + dlon, v[1] + dlat, v[2], v[3]))
                    }
                    ShiftKind::Geographic3DOffset => {
                        let plan = self.hplan.as_ref().ok_or(ProjError::UnsupportedOperation)?;
                        let (dlon, dlat) = geotiff_horizontal_offset(plan, &shifts);
                        let dz = self.geog3d_dz(&shifts);
                        Ok(Coord::new(v[0] + dlon, v[1] + dlat, v[2] + dz, v[3]))
                    }
                    ShiftKind::Xyz => {
                        let dz = if self.no_z_transform || shifts.len() < 3 {
                            0.0
                        } else {
                            shifts[2]
                        };
                        Ok(Coord::new(
                            v[0] + shifts[0].to_radians(),
                            v[1] + shifts[1].to_radians(),
                            v[2] + dz,
                            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);
        }
        match self.kind {
            ShiftKind::Vertical => {
                // No iteration needed for pure vertical shift. Inverse SUBTRACTS
                // the grid value (mirror of the forward addition; PROJ
                // `gridshift.cpp` `grid_apply_internal` with `isVerticalOnly`:
                // `out.z -= shift.z`).
                for gs in &self.grids {
                    if let Some(shifts) =
                        sample_with(self.interpolation, gs, v[1].to_degrees(), v[0].to_degrees())
                    {
                        let z = if self.no_z_transform {
                            v[2]
                        } else {
                            v[2] - shifts[0]
                        };
                        return Ok(Coord::new(v[0], v[1], z, v[3]));
                    }
                }
                Err(ProjError::OutsideGrid)
            }
            _ => {
                // Iterative inversion for horizontal components.
                let mut lon_r = v[0];
                let mut lat_r = v[1];
                'outer: for gs in &self.grids {
                    if sample_with(
                        self.interpolation,
                        gs,
                        lat_r.to_degrees(),
                        lon_r.to_degrees(),
                    )
                    .is_none()
                    {
                        continue;
                    }
                    for _ in 0..MAX_ITER {
                        let shifts = match sample_with(
                            self.interpolation,
                            gs,
                            lat_r.to_degrees(),
                            lon_r.to_degrees(),
                        ) {
                            Some(s) => s,
                            None => break 'outer,
                        };
                        let (new_lon, new_lat) = match self.kind {
                            ShiftKind::HorizontalArcSec => (
                                v[0] - shifts[1] * ARCSEC_TO_RAD,
                                v[1] - shifts[0] * ARCSEC_TO_RAD,
                            ),
                            ShiftKind::HorizontalGeoTiff | ShiftKind::Geographic3DOffset => {
                                let plan =
                                    self.hplan.as_ref().ok_or(ProjError::UnsupportedOperation)?;
                                let (dlon, dlat) = geotiff_horizontal_offset(plan, &shifts);
                                (v[0] - dlon, v[1] - dlat)
                            }
                            ShiftKind::Xyz => {
                                (v[0] - shifts[0].to_radians(), v[1] - shifts[1].to_radians())
                            }
                            // Vertical is handled by the outer match arm above;
                            // this arm is logically unreachable, but we return
                            // the input coordinates unchanged rather than panicking.
                            ShiftKind::Vertical => (v[0], v[1]),
                        };
                        let dlon = (new_lon - lon_r).abs();
                        let dlat = (new_lat - lat_r).abs();
                        lon_r = new_lon;
                        lat_r = new_lat;
                        if dlon < ITER_TOL && dlat < ITER_TOL {
                            break;
                        }
                    }
                    // Compute the dz correction for the kinds that carry one.
                    if let Some(shifts) = sample_with(
                        self.interpolation,
                        gs,
                        lat_r.to_degrees(),
                        lon_r.to_degrees(),
                    ) {
                        let dz = match self.kind {
                            ShiftKind::Xyz => {
                                if self.no_z_transform || shifts.len() < 3 {
                                    0.0
                                } else {
                                    shifts[2]
                                }
                            }
                            ShiftKind::Geographic3DOffset => self.geog3d_dz(&shifts),
                            _ => 0.0,
                        };
                        return Ok(Coord::new(lon_r, lat_r, v[2] - dz, v[3]));
                    }
                    return Ok(Coord::new(lon_r, lat_r, v[2], v[3]));
                }
                Err(ProjError::OutsideGrid)
            }
        }
    }

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

impl GridShift {
    /// Vertical (metre) shift for a [`ShiftKind::Geographic3DOffset`] grid: the
    /// value of the resolved `idx_z` band, honoring `+no_z_transform`.
    fn geog3d_dz(&self, shifts: &[f64]) -> f64 {
        if self.no_z_transform {
            return 0.0;
        }
        match self.idx_z {
            Some(i) => shifts.get(i).copied().unwrap_or(0.0),
            None => 0.0,
        }
    }
}

/// Parse `+interpolation=` into an [`Interpolation`], returning `None` when the
/// parameter is absent (so the grid's own `interpolation_method` metadata, then
/// the bilinear default, can apply) and erroring on unsupported values exactly
/// as PROJ `gridshift.cpp` does.
fn parse_interpolation(p: &TransParams) -> ProjResult<Option<Interpolation>> {
    match p.params.get_str("interpolation") {
        None => Ok(None),
        Some("bilinear") => Ok(Some(Interpolation::Bilinear)),
        Some("biquadratic") => Ok(Some(Interpolation::Biquadratic)),
        Some(_) => Err(ProjError::IllegalArgValue),
    }
}

/// Interpret a grid-level `interpolation_method` GDAL_METADATA string into an
/// [`Interpolation`], mirroring PROJ's fallback (`bilinear` when empty; error on
/// any other value).
fn interpolation_from_metadata(method: Option<&str>) -> ProjResult<Interpolation> {
    match method {
        None | Some("") | Some("bilinear") => Ok(Interpolation::Bilinear),
        Some("biquadratic") => Ok(Interpolation::Biquadratic),
        Some(_) => Err(ProjError::IllegalArgValue),
    }
}

/// Build a [`TransBuild`] wrapping a radians-in/radians-out [`GridShift`].
fn build(shift: GridShift) -> TransBuild {
    TransBuild::new(Box::new(shift), IoUnits::Radians, IoUnits::Radians)
}

/// Construct a `gridshift` transform from parsed parameters.
///
/// Splits `+grids=` on commas (honoring the per-entry `@` optional prefix),
/// resolves each grid, and determines the shift kind. For GeoTIFF grids the
/// grid-level `TYPE` GDAL_METADATA item drives dispatch
/// (`HORIZONTAL_OFFSET` / `GEOGRAPHIC_3D_OFFSET` /
/// `VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL` / `VERTICAL_OFFSET_VERTICAL_TO_VERTICAL`
/// / `ELLIPSOIDAL_HEIGHT_OFFSET`), matching PROJ `gridshift.cpp`'s
/// `checkGridTypes`; when `TYPE` is absent the band-count heuristic (NTv2 → 2-band
/// horizontal, 1-band → vertical, 3+ → xyz) is used. `+interpolation=` and
/// `+no_z_transform` are honored.
pub fn new(p: &TransParams) -> ProjResult<TransBuild> {
    let grid_name = p.params.get_str("grids").ok_or(ProjError::MissingArg)?;
    let interp_param = parse_interpolation(p)?;
    let no_z_transform = p.params.get_bool("no_z_transform");

    let resolved = resolve_grid_list(p.registry, grid_name)?;
    // All-optional list whose entries all failed: pass-through transform.
    if resolved.is_empty() {
        return Ok(build(GridShift {
            grids: Vec::new(),
            kind: ShiftKind::Vertical,
            hplan: None,
            idx_z: None,
            interpolation: interp_param.unwrap_or(Interpolation::Bilinear),
            no_z_transform,
        }));
    }

    // Determine the format from the first resolved grid; every entry in a
    // fallback list is parsed the same way (PROJ requires a homogeneous list).
    let (first_name, first_bytes) = &resolved[0];

    if is_tiff(first_bytes) {
        return build_geotiff_dispatch(&resolved, interp_param, no_z_transform);
    }

    // Non-GeoTIFF grids carry no interpolation_method metadata: honor the
    // explicit +interpolation, else default bilinear.
    let interpolation = interp_param.unwrap_or(Interpolation::Bilinear);

    // NTv2 (horizontal). Flatten every sub-grid of every entry into the
    // fallback list.
    if let Ok(grids) = read_ntv2(first_bytes, first_name) {
        if !grids.is_empty() && grids[0].bands.len() >= 2 {
            let mut all = grids;
            for (name, bytes) in &resolved[1..] {
                let more = read_ntv2(bytes, name)?;
                all.extend(more);
            }
            return Ok(build(GridShift {
                grids: all,
                kind: ShiftKind::HorizontalArcSec,
                hplan: None,
                idx_z: None,
                interpolation,
                no_z_transform,
            }));
        }
    }

    // GTX (vertical).
    if let Ok(gs) = read_gtx(first_bytes, first_name) {
        if gs.bands.len() == 1 {
            let mut all = vec![gs];
            for (name, bytes) in &resolved[1..] {
                all.push(read_gtx(bytes, name)?);
            }
            return Ok(build(GridShift {
                grids: all,
                kind: ShiftKind::Vertical,
                hplan: None,
                idx_z: None,
                interpolation,
                no_z_transform,
            }));
        }
    }

    Err(ProjError::FileNotFound)
}

/// Dispatch a GeoTIFF `+grids=` list: read the grid-level `TYPE` metadata of the
/// first grid and build the corresponding [`ShiftKind`], parsing every entry the
/// same way. When `TYPE` is absent, fall back to the band-count heuristic.
///
/// The interpolation kernel is the explicit `+interpolation=` when given,
/// otherwise the grid's own `interpolation_method` GDAL_METADATA item, otherwise
/// bilinear — matching PROJ `gridshift.cpp`'s `grid_interpolate`.
fn build_geotiff_dispatch(
    resolved: &[(String, Vec<u8>)],
    interp_param: Option<Interpolation>,
    no_z_transform: bool,
) -> ProjResult<TransBuild> {
    let (first_name, first_bytes) = &resolved[0];
    let first_gs = read_geotiff(first_bytes, first_name)?;
    let md = read_geotiff_gdal_metadata(first_bytes)?;
    let type_meta = md.as_ref().and_then(|m| m.get("TYPE").map(str::to_string));
    let interpolation = match interp_param {
        Some(i) => i,
        None => {
            interpolation_from_metadata(md.as_ref().and_then(|m| m.get("interpolation_method")))?
        }
    };

    // Helper: expand every `+grids=` entry into its full parent/child sub-grid
    // hierarchy (`read_geotiff_hierarchy` — every georeferenced IFD, ordered
    // child-first) and apply `map` to each resulting grid. A multi-IFD GeoTIFF
    // (e.g. de_geosn_NTv2_SN.tif: coarse RD83 parent + densified SN_* children)
    // therefore contributes every sub-grid, finest-first, so the
    // first-containing-grid-wins scan selects the finest sub-grid covering a
    // point — matching PROJ's `HorizontalShiftGrid::gridAt` recursion.
    let expand = |map: &dyn Fn(GridSet) -> ProjResult<GridSet>| -> ProjResult<Vec<GridSet>> {
        let mut all = Vec::new();
        for (name, bytes) in resolved {
            for gs in read_geotiff_hierarchy(bytes, name)? {
                all.push(map(gs)?);
            }
        }
        Ok(all)
    };

    match type_meta.as_deref() {
        Some("HORIZONTAL_OFFSET") => {
            let plan = resolve_hshift_plan(&first_gs).ok_or(ProjError::UnsupportedOperation)?;
            let grids = expand(&|gs| Ok(gs))?;
            Ok(build(GridShift {
                grids,
                kind: ShiftKind::HorizontalGeoTiff,
                hplan: Some(plan),
                idx_z: None,
                interpolation,
                no_z_transform,
            }))
        }
        Some("GEOGRAPHIC_3D_OFFSET") => {
            let plan = resolve_hshift_plan(&first_gs).ok_or(ProjError::UnsupportedOperation)?;
            let idx_z = resolve_vertical_band(&first_gs).ok_or(ProjError::FileNotFound)?;
            let grids = expand(&|gs| Ok(gs))?;
            Ok(build(GridShift {
                grids,
                kind: ShiftKind::Geographic3DOffset,
                hplan: Some(plan),
                idx_z: Some(idx_z),
                interpolation,
                no_z_transform,
            }))
        }
        Some("VERTICAL_OFFSET_GEOGRAPHIC_TO_VERTICAL")
        | Some("VERTICAL_OFFSET_VERTICAL_TO_VERTICAL")
        | Some("ELLIPSOIDAL_HEIGHT_OFFSET") => {
            let grids = expand(&|gs| vertical_gridset(gs))?;
            Ok(build(GridShift {
                grids,
                kind: ShiftKind::Vertical,
                hplan: None,
                idx_z: None,
                interpolation,
                no_z_transform,
            }))
        }
        // An explicit but unhandled TYPE (e.g. VELOCITY, which PROJ routes to
        // the deformation model, not generic gridshift) is an error, matching
        // PROJ's "Unhandled value for TYPE metadata item".
        Some(_) => Err(ProjError::UnsupportedOperation),
        // No TYPE metadata: fall back to the band-count heuristic of IFD0.
        None => match first_gs.bands.len() {
            2 => {
                let plan = resolve_hshift_plan(&first_gs).ok_or(ProjError::UnsupportedOperation)?;
                let grids = expand(&|gs| Ok(gs))?;
                Ok(build(GridShift {
                    grids,
                    kind: ShiftKind::HorizontalGeoTiff,
                    hplan: Some(plan),
                    idx_z: None,
                    interpolation,
                    no_z_transform,
                }))
            }
            1 => {
                let grids = expand(&|gs| Ok(gs))?;
                Ok(build(GridShift {
                    grids,
                    kind: ShiftKind::Vertical,
                    hplan: None,
                    idx_z: None,
                    interpolation,
                    no_z_transform,
                }))
            }
            _ => {
                let grids = expand(&|gs| Ok(gs))?;
                Ok(build(GridShift {
                    grids,
                    kind: ShiftKind::Xyz,
                    hplan: None,
                    idx_z: None,
                    interpolation,
                    no_z_transform,
                }))
            }
        },
    }
}

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

    fn build_ntv2(lat_shift: f32, lon_shift: f32) -> Vec<u8> {
        // Minimal 2×2 NTv2 covering lat 0-1°, lon -1-0°
        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]);
        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]);
        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());
            buf.extend_from_slice(&0.0f32.to_le_bytes());
        }
        buf
    }

    fn build_gtx(shift: f32) -> Vec<u8> {
        let mut buf = Vec::new();
        buf.extend_from_slice(&0.0f64.to_be_bytes()); // south_lat=0
        buf.extend_from_slice(&0.0f64.to_be_bytes()); // west_lon=0
        buf.extend_from_slice(&1.0f64.to_be_bytes()); // lat_inc=1
        buf.extend_from_slice(&1.0f64.to_be_bytes()); // lon_inc=1
        buf.extend_from_slice(&2i32.to_be_bytes()); // rows=2
        buf.extend_from_slice(&2i32.to_be_bytes()); // cols=2
        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"
        }
    }

    #[test]
    fn test_gridshift_detects_ntv2() {
        // NTv2 grid covers lat 0-1°, lon -1-0° (NTv2 positive-west convention)
        let data = build_ntv2(1800.0, 900.0);
        let mut map = std::collections::HashMap::new();
        map.insert("grid.gsb".to_string(), data);
        let reg = InMemReg(map);
        let tp = TestParams {
            grids: "grid.gsb".to_string(),
        };
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let p = crate::TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: Some(&reg),
        };
        let tb = new(&p).unwrap();
        // Point inside the grid: lat=0.5°, lon=-0.5°
        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();
        // NTv2 lat_shift=1800 arcsec → +0.5 deg (positive-north).
        // Raw lon_shift=900 arcsec is positive-WEST; read_ntv2 negates it to
        // east-positive (-0.25 deg), so lon moves WEST by 0.25 deg.
        let expected_lat = (0.5 + 0.5) * DEG_TO_RAD;
        let expected_lon = (-0.5 - 0.25) * DEG_TO_RAD;
        assert!((out.v()[1] - expected_lat).abs() < 1e-9, "lat");
        assert!((out.v()[0] - expected_lon).abs() < 1e-9, "lon");
    }

    #[test]
    fn test_gridshift_detects_gtx() {
        // GTX grid covers lat 0-1°, lon 0-1°
        let data = build_gtx(5.0);
        let mut map = std::collections::HashMap::new();
        map.insert("grid.gtx".to_string(), data);
        let reg = InMemReg(map);
        let tp = TestParams {
            grids: "grid.gtx".to_string(),
        };
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        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();
        // The generic `gridshift` operator ADDS the vertical grid value in the
        // forward direction (PROJ `gridshift.cpp` `out.z += shift.z`), which is
        // the OPPOSITE of the dedicated `vgridshift` operator.
        assert!(
            (out.v()[2] - 5.0).abs() < 1e-6,
            "z should be +5.0 (generic gridshift forward adds), got {}",
            out.v()[2]
        );
    }

    #[test]
    fn resolve_grid_bytes_registry_first_then_file_not_found() {
        // E4: registered in-memory bytes win (registry-first), and a grid
        // present nowhere maps to FileNotFound. The name is unique enough that
        // no local resource directory can hold it.
        let mut map = std::collections::HashMap::new();
        map.insert("reg.gsb".to_string(), vec![1u8, 2, 3]);
        let reg = InMemReg(map);
        assert_eq!(
            resolve_grid_bytes(Some(&reg), "reg.gsb").unwrap(),
            vec![1u8, 2, 3]
        );
        assert_eq!(
            resolve_grid_bytes(Some(&reg), "e4_absent_grid_9d3f7a.gsb").err(),
            Some(ProjError::FileNotFound)
        );
        assert_eq!(
            resolve_grid_bytes(None, "e4_absent_grid_9d3f7a.gsb").err(),
            Some(ProjError::FileNotFound)
        );
    }

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

    // ------------------------------------------------------------------
    // Semantics-aware GeoTIFF horizontal-offset application (offline; no
    // PROJ_DATA / real grid required).
    // ------------------------------------------------------------------

    use oxiproj_grids::{BandSemantics, GridBand, GridExtent, GridSource};

    /// Build a uniform 2-band GeoTIFF-style horizontal-offset grid covering
    /// lat [39,41], lon [-81,-79], with the given band semantics/values.
    fn geotiff_hshift_gridset(
        lat_role: BandRole,
        lat_unit: BandUnit,
        lat_val: f32,
        lon_role: BandRole,
        lon_unit: BandUnit,
        lon_positive: BandPositive,
        lon_val: f32,
    ) -> GridSet {
        let extent = GridExtent {
            ll_lat: 39.0,
            ll_lon: -81.0,
            ur_lat: 41.0,
            ur_lon: -79.0,
            lat_inc: 2.0,
            lon_inc: 2.0,
        };
        let band0 = GridBand {
            extent: extent.clone(),
            values: vec![lat_val; 4],
            semantics: BandSemantics {
                role: lat_role,
                unit: lat_unit,
                positive: BandPositive::Unknown,
            },
        };
        let band1 = GridBand {
            extent,
            values: vec![lon_val; 4],
            semantics: BandSemantics {
                role: lon_role,
                unit: lon_unit,
                positive: lon_positive,
            },
        };
        GridSet {
            bands: vec![band0, band1],
            source: GridSource::GeoTiff {
                path: "mem.tif".into(),
            },
        }
    }

    #[test]
    fn geotiff_hshift_no_band_swap_arcsec_positive_east() {
        // band0=latitude_offset (+3600 arcsec = +1°), band1=longitude_offset
        // (+3600 arcsec = +1°, positive east). At (-80,40): lat→41, lon→-79.
        let gs = geotiff_hshift_gridset(
            BandRole::LatitudeOffset,
            BandUnit::ArcSecond,
            3600.0,
            BandRole::LongitudeOffset,
            BandUnit::ArcSecond,
            BandPositive::East,
            3600.0,
        );
        let tb = build_horizontal_geotiff(gs).unwrap();
        let out = tb
            .operation
            .forward_4d(Coord::new(-80.0 * DEG_TO_RAD, 40.0 * DEG_TO_RAD, 0.0, 0.0))
            .unwrap();
        let v = out.v();
        assert!(
            (v[1].to_degrees() - 41.0).abs() < 1e-9,
            "lat += band0 (no swap): {}",
            v[1].to_degrees()
        );
        assert!(
            (v[0].to_degrees() - (-79.0)).abs() < 1e-9,
            "lon += band1 east: {}",
            v[0].to_degrees()
        );
    }

    #[test]
    fn geotiff_hshift_positive_west_negates_longitude() {
        // Same magnitudes, but longitude positive=west must subtract: lon→-81.
        let gs = geotiff_hshift_gridset(
            BandRole::LatitudeOffset,
            BandUnit::ArcSecond,
            0.0,
            BandRole::LongitudeOffset,
            BandUnit::ArcSecond,
            BandPositive::West,
            3600.0,
        );
        let tb = build_horizontal_geotiff(gs).unwrap();
        let out = tb
            .operation
            .forward_4d(Coord::new(-80.0 * DEG_TO_RAD, 40.0 * DEG_TO_RAD, 0.0, 0.0))
            .unwrap();
        let v = out.v();
        assert!(
            (v[0].to_degrees() - (-81.0)).abs() < 1e-9,
            "lon -= band1 (positive west): {}",
            v[0].to_degrees()
        );
    }

    #[test]
    fn geotiff_hshift_degree_unit_and_default_band_order() {
        // No DESCRIPTION roles (Unknown) => PROJ NTv2-inspired default:
        // band0=lat, band1=lon. Degree unit: +0.5° each. At (-80,40):
        // lat→40.5, lon→-79.5.
        let gs = geotiff_hshift_gridset(
            BandRole::Unknown,
            BandUnit::Degree,
            0.5,
            BandRole::Unknown,
            BandUnit::Degree,
            BandPositive::Unknown,
            0.5,
        );
        let tb = build_horizontal_geotiff(gs).unwrap();
        let out = tb
            .operation
            .forward_4d(Coord::new(-80.0 * DEG_TO_RAD, 40.0 * DEG_TO_RAD, 0.0, 0.0))
            .unwrap();
        let v = out.v();
        assert!((v[1].to_degrees() - 40.5).abs() < 1e-9, "lat default band0");
        assert!(
            (v[0].to_degrees() - (-79.5)).abs() < 1e-9,
            "lon default band1"
        );
    }

    #[test]
    fn geotiff_hshift_round_trips() {
        let gs = geotiff_hshift_gridset(
            BandRole::LatitudeOffset,
            BandUnit::ArcSecond,
            120.0,
            BandRole::LongitudeOffset,
            BandUnit::ArcSecond,
            BandPositive::East,
            -240.0,
        );
        let tb = build_horizontal_geotiff(gs).unwrap();
        let input = Coord::new(-80.0 * DEG_TO_RAD, 40.0 * DEG_TO_RAD, 7.0, 0.0);
        let fwd = tb.operation.forward_4d(input).unwrap();
        let inv = tb.operation.inverse_4d(fwd).unwrap();
        assert!((inv.v()[0] - input.v()[0]).abs() < 1e-11, "lon round-trip");
        assert!((inv.v()[1] - input.v()[1]).abs() < 1e-11, "lat round-trip");
    }

    #[test]
    fn resolve_hshift_plan_rejects_single_offset_channel() {
        // Only latitude_offset described (no longitude_offset) => None.
        let gs = geotiff_hshift_gridset(
            BandRole::LatitudeOffset,
            BandUnit::ArcSecond,
            1.0,
            BandRole::Unknown,
            BandUnit::ArcSecond,
            BandPositive::Unknown,
            1.0,
        );
        assert!(resolve_hshift_plan(&gs).is_none());
    }

    // ------------------------------------------------------------------
    // +interpolation / biquadratic sampler
    // ------------------------------------------------------------------

    struct InterpParams {
        interpolation: Option<String>,
        no_z: bool,
    }
    impl crate::TransParamLookup for InterpParams {
        fn get_str(&self, key: &str) -> Option<&str> {
            match key {
                "grids" => Some("grid.gtx"),
                "interpolation" => self.interpolation.as_deref(),
                _ => 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, key: &str) -> bool {
            key == "no_z_transform" && self.no_z
        }
        fn exists(&self, key: &str) -> bool {
            key == "grids"
                || (key == "interpolation" && self.interpolation.is_some())
                || (key == "no_z_transform" && self.no_z)
        }
    }

    fn interp_of(s: Option<&str>) -> ProjResult<Option<Interpolation>> {
        let ell = oxiproj_core::Ellipsoid::named("WGS84").unwrap();
        let tp = InterpParams {
            interpolation: s.map(str::to_string),
            no_z: false,
        };
        let p = TransParams {
            ellipsoid: &ell,
            params: &tp,
            registry: None,
        };
        parse_interpolation(&p)
    }

    #[test]
    fn parse_interpolation_maps_and_rejects() {
        // Absent => None (so grid metadata / bilinear default can apply).
        assert_eq!(interp_of(None).unwrap(), None);
        assert_eq!(
            interp_of(Some("bilinear")).unwrap(),
            Some(Interpolation::Bilinear)
        );
        assert_eq!(
            interp_of(Some("biquadratic")).unwrap(),
            Some(Interpolation::Biquadratic)
        );
        assert_eq!(
            interp_of(Some("bicubic")).err(),
            Some(ProjError::IllegalArgValue)
        );
    }

    #[test]
    fn interpolation_from_metadata_maps_and_rejects() {
        assert_eq!(
            interpolation_from_metadata(None).unwrap(),
            Interpolation::Bilinear
        );
        assert_eq!(
            interpolation_from_metadata(Some("bilinear")).unwrap(),
            Interpolation::Bilinear
        );
        assert_eq!(
            interpolation_from_metadata(Some("biquadratic")).unwrap(),
            Interpolation::Biquadratic
        );
        assert_eq!(
            interpolation_from_metadata(Some("nearest")).err(),
            Some(ProjError::IllegalArgValue)
        );
    }

    /// Build a single-band 5×5 grid over lon[0,4], lat[0,4] (1° spacing) whose
    /// value follows a separable quadratic `g(col, south_row)`; biquadratic
    /// interpolation must then reproduce `g` exactly at fractional positions.
    fn quadratic_gridset() -> GridSet {
        use oxiproj_grids::{BandSemantics, GridBand, GridExtent, GridSource};
        let n = 5usize;
        let g = |col: f64, srow: f64| {
            1.0 + 2.0 * col - 3.0 * srow + 0.5 * col * col + 0.25 * srow * srow
        };
        let mut values = Vec::with_capacity(n * n);
        for north_row in 0..n {
            let srow = (n - 1 - north_row) as f64;
            for col in 0..n {
                values.push(g(col as f64, srow) as f32);
            }
        }
        let extent = GridExtent {
            ll_lat: 0.0,
            ll_lon: 0.0,
            ur_lat: (n - 1) as f64,
            ur_lon: (n - 1) as f64,
            lat_inc: 1.0,
            lon_inc: 1.0,
        };
        GridSet {
            bands: vec![GridBand {
                extent,
                values,
                semantics: BandSemantics::default(),
            }],
            source: GridSource::Memory,
        }
    }

    #[test]
    fn biquadratic_reproduces_quadratic_exactly() {
        let gs = quadratic_gridset();
        let g = |col: f64, srow: f64| {
            1.0 + 2.0 * col - 3.0 * srow + 0.5 * col * col + 0.25 * srow * srow
        };
        // Query at (lon=2.3, lat=1.7) => col=2.3, south_row=1.7.
        let got = sample_grid_biquadratic(&gs, 1.7, 2.3).unwrap()[0];
        let expected = g(2.3, 1.7);
        assert!(
            (got - expected).abs() < 1e-6,
            "biquadratic exact for quadratic: got {got}, expected {expected}"
        );
        // Bilinear does NOT reproduce a curved quadratic — sanity that the two
        // kernels actually differ here.
        let bilinear = sample_grid(&gs, 1.7, 2.3).unwrap()[0];
        assert!(
            (bilinear - expected).abs() > 1e-3,
            "bilinear should differ from the quadratic truth"
        );
    }

    #[test]
    fn biquadratic_falls_back_to_bilinear_for_small_grids() {
        use oxiproj_grids::{BandSemantics, GridBand, GridExtent, GridSource};
        // 2×2 grid: too small for a 3×3 window, must match bilinear.
        let extent = GridExtent {
            ll_lat: 0.0,
            ll_lon: 0.0,
            ur_lat: 1.0,
            ur_lon: 1.0,
            lat_inc: 1.0,
            lon_inc: 1.0,
        };
        let gs = GridSet {
            bands: vec![GridBand {
                extent,
                values: vec![1.0, 2.0, 3.0, 4.0],
                semantics: BandSemantics::default(),
            }],
            source: GridSource::Memory,
        };
        let bq = sample_grid_biquadratic(&gs, 0.5, 0.5).unwrap()[0];
        let bl = sample_grid(&gs, 0.5, 0.5).unwrap()[0];
        assert!(
            (bq - bl).abs() < 1e-12,
            "2×2 biquadratic falls back to bilinear"
        );
    }

    // ------------------------------------------------------------------
    // resolve_grid_list (@-optional comma fallback)
    // ------------------------------------------------------------------

    #[test]
    fn resolve_grid_list_honors_optional_prefix_and_order() {
        let mut map = std::collections::HashMap::new();
        map.insert("a.gsb".to_string(), vec![1u8]);
        map.insert("b.gsb".to_string(), vec![2u8]);
        let reg = InMemReg(map);

        // Optional-missing then present: skips the first, keeps the second.
        let list = resolve_grid_list(Some(&reg), "@missing.gsb,b.gsb").unwrap();
        assert_eq!(list.len(), 1);
        assert_eq!(list[0].0, "b.gsb");

        // Order preserved for two present grids.
        let list = resolve_grid_list(Some(&reg), "a.gsb,b.gsb").unwrap();
        assert_eq!(
            list.iter().map(|(n, _)| n.as_str()).collect::<Vec<_>>(),
            ["a.gsb", "b.gsb"]
        );

        // Required missing aborts.
        assert_eq!(
            resolve_grid_list(Some(&reg), "a.gsb,missing.gsb").err(),
            Some(ProjError::FileNotFound)
        );

        // All-optional-absent => empty list (pass-through).
        assert!(resolve_grid_list(Some(&reg), "@x.gsb,@y.gsb")
            .unwrap()
            .is_empty());
    }

    // ------------------------------------------------------------------
    // Vertical band selection
    // ------------------------------------------------------------------

    #[test]
    fn resolve_vertical_band_prefers_geoid_role() {
        use oxiproj_grids::{BandSemantics, GridBand, GridExtent, GridSource};
        let extent = GridExtent {
            ll_lat: 0.0,
            ll_lon: 0.0,
            ur_lat: 1.0,
            ur_lon: 1.0,
            lat_inc: 1.0,
            lon_inc: 1.0,
        };
        let band = |role: BandRole| GridBand {
            extent: extent.clone(),
            values: vec![0.0; 4],
            semantics: BandSemantics {
                role,
                unit: BandUnit::Metre,
                positive: BandPositive::Unknown,
            },
        };
        // Band 1 is the geoid channel => selected and swapped to index 0.
        let gs = GridSet {
            bands: vec![band(BandRole::Unknown), band(BandRole::GeoidUndulation)],
            source: GridSource::Memory,
        };
        assert_eq!(resolve_vertical_band(&gs), Some(1));
        let reduced = vertical_gridset(gs).unwrap();
        assert_eq!(reduced.bands.len(), 1);
        assert_eq!(reduced.bands[0].semantics.role, BandRole::GeoidUndulation);

        // No described vertical channel => default to band 0.
        let gs2 = GridSet {
            bands: vec![band(BandRole::Unknown)],
            source: GridSource::Memory,
        };
        assert_eq!(resolve_vertical_band(&gs2), Some(0));
    }

    // ------------------------------------------------------------------
    // Time-gate helpers
    // ------------------------------------------------------------------

    #[test]
    fn time_gate_matches_proj_condition() {
        // Either bound zero => always apply.
        assert!(time_gate_applies(0.0, 0.0, 1234.0));
        assert!(time_gate_applies(2000.0, 0.0, 3000.0));
        assert!(time_gate_applies(0.0, 2010.0, 3000.0));
        // Both set: apply only when t < t_epoch and t_final > t_epoch.
        assert!(time_gate_applies(2000.0, 2010.0, 1995.0));
        assert!(!time_gate_applies(2000.0, 2010.0, 2005.0));
        assert!(!time_gate_applies(2000.0, 2010.0, 2000.0));
        // Degenerate bracket (t_final <= t_epoch) never applies for finite t.
        assert!(!time_gate_applies(2010.0, 2000.0, 1995.0));
    }

    #[test]
    fn civil_date_round_trips_known_days() {
        // 2020-02-29 is day 18321 since 1970-01-01 (a leap day).
        assert_eq!(civil_from_days(18_321), (2020, 2, 29));
        // tm_yday is 0-based: Jan 1 => 0, Feb 29 (leap) => 59.
        assert_eq!(day_of_year(2020, 1, 1), 0);
        assert_eq!(day_of_year(2020, 2, 29), 59);
        assert_eq!(day_of_year(2021, 3, 1), 59); // non-leap
    }
}