Skip to content

Netherlands RD (EPSG:28992) : inaccuracies of 100m when using PROJJSON or WKT2 #549

@algil

Description

@algil

The datum shift rotation parameters for this transformation in the proj4js source code are in microradians. They should be in arcseconds.

e.g.
1/ proj4.defs("EPSG:28992", {.....}); // {....} is the PROJJSON or WKT2 for 28892
proj4("EPSG:28992", {x: lon y: lat});
2/ Using desktop PROJ do:
echo lat lon | cs2cs WGS84 EPSG:28992

The resulting x/y coordinates differ by about 100 meters for points in Amsterdam.

3/ The usual +proj definition generally used in the past for Netherlands RD is :
+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs +type=crs
If this is used in PROJ or proj4js the result matches that in 2/ i.e they are correct.

Analysis:-
The newer official values for the datum shift can be found in https://berthub.eu/articles/RDNAPTRANS2018_v220627.pdf. It specifies that the shift should be achieved by grid based method (eg using .tif files) but it does give the shift parameters for both ETRS89-->RD (at page 8) and RD-->ETRS89 (page 24/25) which actually differ slightly, however it says (bottom of page 7) either could be used for both directions without any noticeable error for most software apps. The difference between WTS84 and ETRS89 is less than 1 meter so can also be ignored for most purposes. The RD-->ETRS89 parameters are given as :
𝑡𝑋 = +565.7381 m translation in direction of X axis
𝑡𝑌 = +50.4018 m translation in direction of Y axis
𝑡𝑍 = +465.2904 m translation in direction of Z axis
𝛼 =+1.91514∙10−6 rad = +0.395026″ rotation angle around X axis
𝛽 =−1.60363∙10−6 rad = −0.330772″ rotation angle around Y axis
𝛾 =+9.09546∙10−6 rad = +1.876074″ rotation angle around Z axis
𝛿 =+4.07244∙10−6 scale difference (dimensionless)

So in +proj form (allowing for the rotation direction they used) :- +towgs84 = 565.7381,50.4018,465.2904,-0.395026,0.330772,-1.876074,4.07244
NB !!! Using arcseconds for the rotation units !!

The full string would be :-
+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.7381,50.4018,465.2904,-0.395026,0.330772,-1.876074,4.07244 +units=m +no_defs +type=crs

If this is used in the the examples above the results all agree with the correct PROJ value. PROJ by default uses the grid based method but if told to use the +towgs84 it uses that and gives the correct result , as does proj4js.

CAUSE:
In the proj4js source code Datum.js at line 161 is the towgs84 used for PROJJSON/WKT2: (the RD base_crs has ID 4289)
EPSG_4289: {
towgs84: '565.7381,50.4018,465.2904,-1.91514,1.60363,-9.09546,4.07244'
},

The parameters are the newer RDNAPTRANS2018 values but the rotations are entered in microradians. They should be in arcseconds.

I amended my source code to use arcseconds and retried and now got the correct result. i.e. the line becomes:-
towgs84: '565.7381,50.4018,465.2904,-0.395026,0.330772,-1.876074,4.07244'

I also eventually worked out how to use the BoundCRS method of specifying datum shifts and and that also gives correct results if the correct towgs84 parms are used.
EPSG.io for EPSG:4289 also gives the Proj4js string with the towgs84 rotation in microradians - so it looks like they also are not respecting the units the rotations are recorded in.

By the way, it looks like the entries in Datum.js for the following might also suffer from the same issue:
EPSG_4231
EPSG_4660

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions