Different (~100m) results from transforming WGS84 -> Belgian lambert 72 depending on definition used

[piemonkey]

I've used the bug report template as this at the very least is unexpected behaviour for an unskilled user, if not an actual bug.

When digging into differences of ~100m between transformations from WGS84 to Belgian Lambert 72 given by different tools, I've spotted that the result given by cs2cs for these transformations depends on which definition of the Belgian Lambert 72 CRS is used.

Example of problem

Different Belgian Lambert 72 coordinates obtained when using the EPSG code or the OGC WKT2 definition compared to the OGC WKT definition.

## Using epsg codes
echo 4.3839 50.8428 | cs2cs +init=epsg:4326 +to +init=epsg:31370
# 151066.89	170265.84 0.00

## Using OGC WKT2 definition from epsg.io
echo 4.3839 50.8428 | cs2cs +init=epsg:4326 +to 'PROJCRS["BD72 / Belgian Lambert 72",BASEGEOGCRS["BD72",DATUM["Reseau National Belge 1972",ELLIPSOID["International 1924",6378388,297,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4313]],CONVERSION["Belgian Lambert 72",METHOD["Lambert Conic Conformal (2SP)",ID["EPSG",9802]],PARAMETER["Latitude of false origin",90,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",4.36748666666667,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",51.1666672333333,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",49.8333339,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",150000.013,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",5400088.438,LENGTHUNIT["metre",1],ID["EPSG",8827]]],CS[Cartesian,2],AXIS["easting (X)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing (Y)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["Belgium - onshore."],BBOX[49.5,2.5,51.51,6.4]],ID["EPSG",31370]]'
# 151066.89	170265.84 0.00

## Using OGC WKT definition from epsg.io
echo 4.3839 50.8428 | cs2cs +init=epsg:4326 +to 'PROJCS["BD72 / Belgian Lambert 72",GEOGCS["BD72",DATUM["Reseau_National_Belge_1972",SPHEROID["International 1924",6378388,297],TOWGS84[-106.8686,52.2978,-103.7239,-0.3366,0.457,-1.8422,-1.2747]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4313"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",90],PARAMETER["central_meridian",4.36748666666667],PARAMETER["standard_parallel_1",51.1666672333333],PARAMETER["standard_parallel_2",49.8333339],PARAMETER["false_easting",150000.013],PARAMETER["false_northing",5400088.438],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","31370"]]'         
# 151124.61	170295.56 0.00

Problem description

This seems to be related to the inclusion of the TOWGS84 parameters in the WKT definition (which technically makes it a description of a transform rather than a CRS?) as this is reflected in the proj strings generated by projinfo for the different definitions (which are identical except for the inclusion of a +towgs84=-106.8686,52.2978,-103.7239,-0.3366,0.457,-1.8422,-1.2747 term when passing the WKT definition to projinfo). I have no idea why this isn't being accounted for by the transform that proj seems to match for the operation, which is as far as I can tell, the correct one.

Expected Output

I would expect that the output of the transformation to match between different definition formats, to within the stated accuracy of the CRSs used.

Environment Information

• PROJ version (proj 9.4.1)
• Arch linux 6.9.3-arch1-1

Installation method

• pacman

[rouault]

This seems to be related to the inclusion of the TOWGS84 parameters in the WKT definition (which technically makes it a description of a transform rather than a CRS?) as this is reflected in the proj strings generated by projinfo for the different definitions (which are identical except for the inclusion of a +towgs84=-106.8686,52.2978,-103.7239,-0.3366,0.457,-1.8422,-1.2747 term when passing the WKT definition to projinfo). I have no idea why this isn't being accounted for by the transform that proj seems to match for the operation, which is as far as I can tell, the correct one.

Contrary to what its name suggests, epsg.io is not an official reference, and in particular, for this use case, it got the values of the TOWGS84 clause wrong.
It uses the parameters of https://epsg.io/15929, which is a Coordinate Frame rotation transformation, unmodified in a TOWGS84[] node, which uses the Position Vector transformation convention! The 2 conventions are very similar, but differ by the sign of the rotation parameters. So the TOWGS84[] clause should be corrected to -106.8686,52.2978,-103.7239,0.3366,-0.457,1.8422,-1.2747 to get the expected values

You'd be better using the WKT strings from https://spatialreference.org/ which are derived from PROJ.

[piemonkey]

Thanks @rouault for your incredibly swift reply.

They really seem to think they're an official reference

The EPSG.io website is built around the official EPSG database maintained by OGP Geomatics Committee (http://www.epsg.org/). The database comprises of very detailed geodetic parameters from a range of sources and authorities.

I'll try to get them to update these values so others don't fall into the trap of believing them and assuming their details are correct.

[jorisvandenbossche]

@piemonkey there are already various issues about it (@rouault commented on maptiler/epsg.io#194, and maptiler/epsg.io#190 seems another duplicate)

Also another issue of Even about it being misleading people in believing the site is a reference: maptiler/epsg.io#171

[piemonkey]

Glad to see people have already tried to correct this. I hope something eventually gets through...

As for spatialreference.org, it seems they don't include the TOWGS84 parameters at all. Similar to projinfo -o WKT1:GDAL epsg:31370. For me, I have the transformation I need now (thanks again!), but I wonder if it would be possible to include those parameters somehow. I'm looking specifically at WKT1 strings as the tool I'm using (proj4js) only accepts WKT1 or proj4.

[klokan]

We are going to fix this bug in EPSG.io. Thank you for spotting it!

Keen to discuss what other actions can be taken to make the experience better - to reduce the pains on OSGEO/PROJ side @rouault - a meetup during FOSS4GE?

[rouault]

@rouault - a meetup during FOSS4GE?

yes, we can try to find a moment to discuss