oxiproj-engine 0.1.2

//! End-to-end differential tests for the newly-wired draft projections
//! (cea, eqc, somerc, moll, sinu, robin), run through the oxiproj-engine
//! public API (create + forward/inverse). Reference values are hardcoded
//! from PROJ 9.7.0 and exercise a-scaling, lon_0, and false easting/northing
//! end to end.

use oxiproj_core::{Coord, DEG_TO_RAD};
use oxiproj_engine::create;

/// Forward a (lon_deg, lat_deg) point through `proj_str`, returning (x_m, y_m).
fn fwd(proj_str: &str, lon_deg: f64, lat_deg: f64) -> (f64, f64) {
    let pj = create(proj_str).expect("create proj");
    let out = pj
        .forward(Coord::new(
            lon_deg * DEG_TO_RAD,
            lat_deg * DEG_TO_RAD,
            0.0,
            0.0,
        ))
        .expect("forward");
    let v = out.v();
    (v[0], v[1])
}

/// Inverse an (x_m, y_m) point through `proj_str`, returning (lon_rad, lat_rad).
fn inv(proj_str: &str, x: f64, y: f64) -> (f64, f64) {
    let pj = create(proj_str).expect("create proj");
    let out = pj.inverse(Coord::new(x, y, 0.0, 0.0)).expect("inverse");
    let v = out.v();
    (v[0], v[1])
}

#[test]
fn somerc_bessel_ch1903() {
    let proj = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel";
    let (lon, lat) = (8.0, 47.0);
    let (x, y) = fwd(proj, lon, lat);
    assert!(
        (x - 642617.5280747189).abs() < 1e-6,
        "somerc_bessel_ch1903 x = {x}"
    );
    assert!(
        (y - 205442.8138998982).abs() < 1e-6,
        "somerc_bessel_ch1903 y = {y}"
    );
    let (lo, la) = inv(proj, x, y);
    assert!(
        (lo - lon * DEG_TO_RAD).abs() < 1e-9,
        "somerc_bessel_ch1903 inv lon = {lo}"
    );
    assert!(
        (la - lat * DEG_TO_RAD).abs() < 1e-9,
        "somerc_bessel_ch1903 inv lat = {la}"
    );
}

#[test]
fn cea_wgs84() {
    let proj = "+proj=cea +lat_ts=30 +lon_0=0 +ellps=WGS84";
    let (lon, lat) = (12.0, 55.0);
    let (x, y) = fwd(proj, lon, lat);
    assert!((x - 1157835.3630107583).abs() < 1e-6, "cea_wgs84 x = {x}");
    assert!(
        (y - 6_005_522.387_616_127).abs() < 1e-6,
        "cea_wgs84 y = {y}"
    );
    let (lo, la) = inv(proj, x, y);
    assert!(
        (lo - lon * DEG_TO_RAD).abs() < 1e-9,
        "cea_wgs84 inv lon = {lo}"
    );
    assert!(
        (la - lat * DEG_TO_RAD).abs() < 1e-9,
        "cea_wgs84 inv lat = {la}"
    );
}

#[test]
fn cea_sphere() {
    let proj = "+proj=cea +lat_ts=30 +lon_0=0 +R=1";
    let (lon, lat) = (12.0, 55.0);
    let (x, y) = fwd(proj, lon, lat);
    assert!((x - 0.181379936423).abs() < 1e-6, "cea_sphere x = {x}");
    assert!((y - 0.945875306555).abs() < 1e-6, "cea_sphere y = {y}");
    let (lo, la) = inv(proj, x, y);
    assert!(
        (lo - lon * DEG_TO_RAD).abs() < 1e-9,
        "cea_sphere inv lon = {lo}"
    );
    assert!(
        (la - lat * DEG_TO_RAD).abs() < 1e-9,
        "cea_sphere inv lat = {la}"
    );
}

#[test]
fn eqc_wgs84() {
    let proj = "+proj=eqc +lat_ts=30 +lon_0=0 +ellps=WGS84";
    let (lon, lat) = (12.0, 55.0);
    let (x, y) = fwd(proj, lon, lat);
    // PROJ 9.8.0 ELLIPSOIDAL eqc: x = nu1*cos(lat_ts)*a*lam, y = a*(M(phi)-M(phi0)). This differs from the PROJ 9.7.0 binary (spherical); documented version change.
    assert!((x - 1157835.3630107583).abs() < 1e-6, "eqc_wgs84 x = {x}");
    assert!(
        (y - 6_097_230.313_124_184).abs() < 1e-6,
        "eqc_wgs84 y = {y}"
    );
    let (lo, la) = inv(proj, x, y);
    assert!(
        (lo - lon * DEG_TO_RAD).abs() < 1e-9,
        "eqc_wgs84 inv lon = {lo}"
    );
    assert!(
        (la - lat * DEG_TO_RAD).abs() < 1e-9,
        "eqc_wgs84 inv lat = {la}"
    );
}

#[test]
fn eqc_sphere_plate_carree() {
    let proj = "+proj=eqc +lon_0=0 +R=1";
    let (lon, lat) = (12.0, 55.0);
    let (x, y) = fwd(proj, lon, lat);
    assert!(
        (x - 0.209439510239).abs() < 1e-6,
        "eqc_sphere_plate_carree x = {x}"
    );
    assert!(
        (y - 0.959931088597).abs() < 1e-6,
        "eqc_sphere_plate_carree y = {y}"
    );
    let (lo, la) = inv(proj, x, y);
    assert!(
        (lo - lon * DEG_TO_RAD).abs() < 1e-9,
        "eqc_sphere_plate_carree inv lon = {lo}"
    );
    assert!(
        (la - lat * DEG_TO_RAD).abs() < 1e-9,
        "eqc_sphere_plate_carree inv lat = {la}"
    );
}

#[test]
fn moll_sphere() {
    let proj = "+proj=moll +lon_0=0 +R=1";
    let (lon, lat) = (12.0, 55.0);
    let (x, y) = fwd(proj, lon, lat);
    assert!((x - 0.133156601284).abs() < 1e-6, "moll_sphere x = {x}");
    assert!((y - 1.001323735770).abs() < 1e-6, "moll_sphere y = {y}");
    let (lo, la) = inv(proj, x, y);
    assert!(
        (lo - lon * DEG_TO_RAD).abs() < 1e-9,
        "moll_sphere inv lon = {lo}"
    );
    assert!(
        (la - lat * DEG_TO_RAD).abs() < 1e-9,
        "moll_sphere inv lat = {la}"
    );
}

#[test]
fn moll_wgs84_a() {
    let proj = "+proj=moll +lon_0=0 +ellps=WGS84";
    let (lon, lat) = (12.0, 55.0);
    let (x, y) = fwd(proj, lon, lat);
    assert!((x - 849291.0454429336).abs() < 1e-6, "moll_wgs84_a x = {x}");
    assert!(
        (y - 6_386_579.968_095_269).abs() < 1e-6,
        "moll_wgs84_a y = {y}"
    );
    let (lo, la) = inv(proj, x, y);
    assert!(
        (lo - lon * DEG_TO_RAD).abs() < 1e-9,
        "moll_wgs84_a inv lon = {lo}"
    );
    assert!(
        (la - lat * DEG_TO_RAD).abs() < 1e-9,
        "moll_wgs84_a inv lat = {la}"
    );
}

#[test]
fn sinu_wgs84() {
    let proj = "+proj=sinu +lon_0=0 +ellps=WGS84";
    let (lon, lat) = (12.0, 55.0);
    let (x, y) = fwd(proj, lon, lat);
    assert!((x - 767929.5515727085).abs() < 1e-6, "sinu_wgs84 x = {x}");
    assert!(
        (y - 6_097_230.313_124_184).abs() < 1e-6,
        "sinu_wgs84 y = {y}"
    );
    let (lo, la) = inv(proj, x, y);
    assert!(
        (lo - lon * DEG_TO_RAD).abs() < 1e-9,
        "sinu_wgs84 inv lon = {lo}"
    );
    assert!(
        (la - lat * DEG_TO_RAD).abs() < 1e-9,
        "sinu_wgs84 inv lat = {la}"
    );
}

#[test]
fn sinu_sphere() {
    let proj = "+proj=sinu +lon_0=0 +R=1";
    let (lon, lat) = (12.0, 55.0);
    let (x, y) = fwd(proj, lon, lat);
    assert!((x - 0.120129567914).abs() < 1e-6, "sinu_sphere x = {x}");
    assert!((y - 0.959931088597).abs() < 1e-6, "sinu_sphere y = {y}");
    let (lo, la) = inv(proj, x, y);
    assert!(
        (lo - lon * DEG_TO_RAD).abs() < 1e-9,
        "sinu_sphere inv lon = {lo}"
    );
    assert!(
        (la - lat * DEG_TO_RAD).abs() < 1e-9,
        "sinu_sphere inv lat = {la}"
    );
}

#[test]
fn robin_sphere() {
    let proj = "+proj=robin +lon_0=0 +R=1";

    // Forward checks against PROJ 9.7.0 reference (table tol 1e-4).
    let (x1, y1) = fwd(proj, 12.0, 55.0);
    assert!(
        (x1 - 0.148422341990).abs() < 1e-4,
        "robin_sphere x(12,55) = {x1}"
    );
    assert!(
        (y1 - 0.915371909463).abs() < 1e-4,
        "robin_sphere y(12,55) = {y1}"
    );

    let (x2, y2) = fwd(proj, 170.0, 85.0);
    assert!(
        (x2 - 1.440881763768).abs() < 1e-4,
        "robin_sphere x(170,85) = {x2}"
    );
    assert!(
        (y2 - 1.319980067271).abs() < 1e-4,
        "robin_sphere y(170,85) = {y2}"
    );

    // Round-trip (forward then inverse of the freshly-computed output) for two points.
    for &(lon, lat) in &[(12.0_f64, 55.0_f64), (-60.0_f64, 30.0_f64)] {
        let (x, y) = fwd(proj, lon, lat);
        let (lo, la) = inv(proj, x, y);
        assert!(
            (lo - lon * DEG_TO_RAD).abs() < 1e-9,
            "robin_sphere roundtrip lon({lon},{lat}) = {lo}"
        );
        assert!(
            (la - lat * DEG_TO_RAD).abs() < 1e-9,
            "robin_sphere roundtrip lat({lon},{lat}) = {la}"
        );
    }
}