Bug in conversion from "OGC:CRS84" to "EPSG:3395"

[natassaf]

Below there is a test in rust that fails. The proj version used is the latest as of today, proj = "0.27.0"

use proj::Proj;
use approx::relative_eq;

pub fn merc_to_coords(input_point: (f64, f64)) -> (f64, f64) {
    let input_point = (input_point.0 as f64, input_point.1 as f64);
    let from = "EPSG:3395";
    let to = "OGC:CRS84";
    let converter = Proj::new_known_crs(from, to, None).unwrap();
    converter.convert(input_point).unwrap()
}

pub fn coords_to_merc(input_point: (f64, f64)) -> (f64, f64) {
    let from = "OGC:CRS84";
    let to = "EPSG:3395";
    let converter = Proj::new_known_crs(&from, &to, None).unwrap();
    let output_point = converter
        .convert(input_point)
        .expect("convertion from coords to merc failed");
    (output_point.0, output_point.1)
}

#[test]
fn test_geoprojector_bug() {
    let p1_merc = (-26471262.0, 3462994.0);
    let p1_coords = merc_to_coords(p1_merc);
    let p1_merc_reconstructed = coords_to_merc(p1_coords);

    let expected_p1_coords = (122.2046103800736, 29.849799842396617);
    assert!(relative_eq!(p1_coords.0, expected_p1_coords.0, max_relative=1e-6)); // pass
    assert!(relative_eq!(p1_coords.1, expected_p1_coords.1, max_relative=1e-6)); // pass

    // p1_merc_reconstructed is (13603754.685578486, 3462993.9999999995) instead of the initial p1_merc
    assert!(relative_eq!(p1_merc.0, p1_merc_reconstructed.0, max_relative=1e-6)); // fails
    assert!(relative_eq!(p1_merc.1, p1_merc_reconstructed.1, max_relative=1e-6)); // fails
}

[lnicola]

I suspect they're equivalent. I couldn't make cs2cs work with OGC:CRS84 never mind, I made a typo, but I guess they're almost the same thing anyway, but consider:

$ echo "-26471262.0 3462994.0" | cs2cs EPSG:3395 EPSG:4326 -f "%.8f"
29.84979984 122.20460755 0.00000000
$ echo "29.84979984 122.20460755" | cs2cs EPSG:4326 EPSG:3395 -f "%.8f"
13603754.68505784   3462993.99969395 0.00000000
$ echo "13603754.68505784 3462993.99969395" | cs2cs EPSG:3395 EPSG:4326 -f "%.8f"
29.84979984 122.20460755 0.00000000
# same point

Keep in mind that your point appears to be outside of the EPSG:3395 bounds:

$ projinfo EPSG:3395
PROJ.4 string:
+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs

WKT2:2019 string:
PROJCRS["WGS 84 / World Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["World Mercator",
        METHOD["Mercator (variant A)",
            ID["EPSG",9804]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Very small scale conformal mapping."],
        AREA["World between 80°S and 84°N."],
        BBOX[-80,-180,84,180]],
    ID["EPSG",3395]]

$ echo "-80 -180" | cs2cs EPSG:4326 EPSG:3395 -f "%.8f"
-20037508.34278924  -15496570.73972372 0.00000000

# -26471262.0 < -20037508.34278924

[urschrei]

Have you confirmed that the result is correct when run with cs2cs? someone else got there first.

[lnicola]

Yeah, see my edit:

$ echo "-26471262.0 3462994.0" | cs2cs EPSG:3395 EPSG:4326 -f "%.8f"
29.84979984 122.20460755 0.00000000
$ echo "13603754.68505784 3462993.99969395" | cs2cs EPSG:3395 EPSG:4326 -f "%.8f"
29.84979984 122.20460755 0.00000000

[lnicola]

Note for random readers: spatialreference.org (e.g.) has the projected bounds for the area of use, don't do it the hard way like I did.

[natassaf]

Thank you guys! Indeed I believe the issue is that the initial EPSG:3395 coordinates is out of bounds!