1

I'm using PROJ with a WKT projection string to try and match a coordinate system being employed in a 3rd party solution.

I seem to be able to match the 3rd party coordinates when projecting EPSG:3857 coordinates to my WKT definition, but when I try to convert back again the Y value is correct but the X values is completely wrong.

I'm pretty new trying to use WKT projection strings so I'm wondering if there's something wrong with my definition?

This is my custom WKT definition:

PROJCS["PacificIslands",
    GEOGCS["NZGD2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4167"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",-21.2],
    PARAMETER["central_meridian",-159.75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","100001"]]

If I take the following EPSG:3857 coordinates:

X: 22295668.360
Y: -2414760.379

Projecting them using the WKT given above:

X: 3674.005
Y: 909.963

Then if I project them back to EPSG:3857 again I get:

X: -17779348.325
Y: -2414760.379

Any idea what's going on?

I've tried this in both the C++ PROJ 8.0.0 library and using the C# ProjNet library and get the same result.

Here's my C++ PROJ example code in case it's of any use:

#include <stdio.h>
#include <proj.h>

int main(void) {
    PJ_CONTEXT* C;
    PJ* P;

    /* or you may set C=PJ_DEFAULT_CTX if you are sure you will     */
    /* use PJ objects from only one thread                          */
    C = proj_context_create();

    P = proj_create_crs_to_crs(C,
        "EPSG:3857",
        "PROJCS[\"PacificIslands\",GEOGCS[\"NZGD2000\",DATUM[\"New_Zealand_Geodetic_Datum_2000\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6167\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4167\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",-21.2],PARAMETER[\"central_meridian\",-159.75],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AUTHORITY[\"EPSG\",\"100001\"]]",
        NULL);

    PJ_COORD coordFrom, coordTo, coordReversed;
    coordFrom = proj_coord(22295668.3602, -2414760.3791, 0, 0);

    coordTo = proj_trans(P, PJ_FWD, coordFrom);
    printf("ORIGINAL x: %.3f, y: %.3f\n", coordFrom.enu.e, coordFrom.enu.n);
    printf("CONVERTED x: %.3f, y: %.3f\n", coordTo.enu.e, coordTo.enu.n);

    // Reverse again
    coordReversed = proj_trans(P, PJ_INV, coordTo);
    printf("REVERSED BACK x: %.3f, y: %.3f\n", coordReversed.enu.e, coordReversed.enu.n);

    /* Clean up */
    proj_destroy(P);
    proj_context_destroy(C); /* may be omitted in the single threaded case */
    return 0;
}
Improve this question
2
  • 6
    It looks like your original EPSG:3857 coordinates are outside the usual bounds of [-20037508.3428, 20037508.3428] meters (-180 to +180 degrees), as though they wrapped around the antemeridian. In fact, when we calculate 22295668.360 - (2 * 20037508.3428), we get -17779348.3256, the same as the output value. The difference between the two equals the circumference of the Earth (40075016.6856 meters).
    – FSimardGIS
    Commented Mar 10, 2021 at 3:43
  • @FSimardGIS thank you for that, that was the issue! I had sourced my coordinates from bboxfinder.com but it provides incorrect coordinates if you scroll past the 180th meridian before selecting your area. Thank you :)
    – Gavin
    Commented Mar 10, 2021 at 19:57

1 Answer 1

Reset to default
0

Works fine for me with gdaltransform when using valid EPSG:3857 coordinates as @FSimardGIS noticed. I had to remove the last line with a fake EPSG code from your WKT but after that:

gdaltransform -t_srs proj.prj -s_srs epsg:3857
Enter X Y [Z [T]] values separated by space, and press Return.
19333064.690113 -5216915.262195
-2193080.84767025 -2703101.3572082 0

>gdaltransform -s_srs proj.prj -t_srs epsg:3857
Enter X Y [Z [T]] values separated by space, and press Return.
-2193080.84767025 -2703101.3572082 0
19333064.690113 -5216915.262195 0
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.