Transformation error

Lounge
1

Hi everyone,

I'm Dave, and I'm a developer from the U.K.

I have implemented Proj4 in C++ into our software but I'm getting transformation errors compared to

https://epsg.io/transform#s_srs=23037&t_srs=32637&x=564928.5884000&y=4262281.7699000
whether we use the database or well known text.

We're using WKT2 2019 but if I convert to WKT1 and compare my WKT to ESPG (see below) I notice that
I'm not including the TOWGS84 clause.
If I manually add the clause to my WKT, I get the required transformation.

Does anyone have any suggestions as to why there is a discrepancy with my WKT compared to ESPG? I'm very new to all this, so would really appreciate some far more expert advise.

Thank you in advance 

Dave

----------------------------------------------------------
My WKT2
----------------------------------------------------------
PROJCRS["ED50 / UTM zone 37N",
BASEGEOGCRS["ED50",
    DATUM["European Datum 1950",
        ELLIPSOID["International 1924",6378388,297,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4230]],
CONVERSION["UTM zone 37N",
    METHOD["Transverse Mercator",
        ID["EPSG",9807]],
    PARAMETER["Latitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8801]],
    PARAMETER["Longitude of natural origin",39,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["Scale factor at natural origin",0.9996,
        SCALEUNIT["unity",1],
        ID["EPSG",8805]],
    PARAMETER["False easting",500000,
        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["unknown"],
    AREA["Europe - 36°E to 42°E and ED50 by country"],
    BBOX[29.18,36,79.09,42]],
ID["EPSG",23037]]
----------------------------------------------------------
PROJCRS["WGS 84 / UTM zone 37N",
BASEGEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]],
CONVERSION["UTM zone 37N",
    METHOD["Transverse Mercator",
        ID["EPSG",9807]],
    PARAMETER["Latitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8801]],
    PARAMETER["Longitude of natural origin",39,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["Scale factor at natural origin",0.9996,
        SCALEUNIT["unity",1],
        ID["EPSG",8805]],
    PARAMETER["False easting",500000,
        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["unknown"],
    AREA["World - N hemisphere - 36°E to 42°E - by country"],
    BBOX[0,36,84,42]],
ID["EPSG",32637]]
----------------------------------------------------------
My WKT1
----------------------------------------------------------
PROJCS["ED50 / UTM zone 37N",
GEOGCS["ED50",
    DATUM["European_Datum_1950",
        SPHEROID["International 1924",6378388,297,
            AUTHORITY["EPSG","7022"]],
        AUTHORITY["EPSG","6230"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4230"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",39],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","23037"]]
----------------------------------------------------------     
WKT1 from https://epsg.io/23037
----------------------------------------------------------         
PROJCS["ED50 / UTM zone 37N",
GEOGCS["ED50",
    DATUM["European_Datum_1950",
        SPHEROID["International 1924",6378388,297,
            AUTHORITY["EPSG","7022"]],
        TOWGS84[-87,-98,-121,0,0,0,0],
        AUTHORITY["EPSG","6230"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4230"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",39],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","23037"]]
2

Dave, welcome here but I think you are making a mistake. This question is clearly appropriate for the Proj4 mailing list. Here we only indirectly use proj4.

3

I do not understand why the TOWGS84 clause is missing, and I’m not sure if I can help you.

I thought that EPSG and WKT are just two ways of representing the same thing, so they should be completely identical. Using EPSG-codes is just a super handy and short way of saying the same.

Remember also that tranforming ED50 <-> WGS84 may be a pain in the bum. I’ve dealt with this several times. There are several transformation parameters that you may use, depending on e.g. how accurate it should be and in which area you are working. I haven’t grasped the full concept of this, but I’ve had my share of pain.

You’ll see that the ED50 transformation involves three parameters in your case;

TOWGS84[-87,-98,-121,0,0,0,0]

There are more variants of this. E.g. in Norwegian areas, the recommended transformation involves a 7-parameter Bursa-Wolfe transformation. Furthermore, there is a different set of parameteres depending if you are north or south of 62d N (epsg: 1612 and epsg:1613, respectively). That’s how I’ve interpreted it at least. Confusing stuff.

Look up Geomatics Guidance Note 10 Geodetic Transformations Offshore Norway to read more about my example.

You can test different transformations here: https://mygeodata.cloud/cs2cs/ and define the transformation string to test several approaches and see how they match your “wanted final coordinates”. The transformation string is given as Proj.4 text which is just another way of telling how you want to process your coordinates, completely the same as an epsg-code or wkt.

4

Ah, ah. It’s worst than you think.

ED50 and WGS84 are datums. Datum’s involve an ellipsoid describing the earth shape and an origin. Just to make you more dizzy, there is really no such a thing as a WGS84. There are six of them, which may differ as much as 1 meter (see for example).

Now, when we want to transform between datum we must use a transform, which may be a 3 or 7 parameters transform, or use a grid shift for better accuracy (have to look into the Proj4 docs for this). Those towgs transforms have been improving along time so it’s normal, and confusing, to find different parameters to achieve the same transform. Furthermore, local datum exist which try to achieve a local best fit of the ellipsoid to local coordinate system, this also leads to different towgs parameters. Not pretty, I agree.

Recent Proj procedures advise against using towgs clause (and have dropped them from WKT definitions) and use the use the new piping system (haven’t dive into it myself).

5

Thanks Joaquim. I’m going to have nightmares tonight.

Sometimes it feels like you have to have a Ph.D in geodesy just to look at this stuff.

Just imagine all the confused e.g. geologists doing this without having a clue…

6

Australians changed their referencing system because due to plate kinematics the previous reference was already 1.8 meters NorthEast of Antarctica.

7

Yikes!