Projection system parameter_Geodesy

mariamsajadian · 31 mei 2023 om 06:57

Hi There,

We have a partner account with Oracle and we are seeking for accurate projection system parameters for the Netherlands. We have attempted using there parameters in QGIS and the website https://epsg.org/crs_28992/Amersfoort-RD-New.html. However, when using the Python library (pyproject), I encountered a significant error of approximately 150 meters and 270 meters, respectively.

Additionally, we conducted a check using the Oracle spatial database, which resulted in an error of around 40 meters.

I would greatly appreciate it if you could provide me with information on how to obtain precise geodetic solutions for the Netherlands. Could you also inform me about the organization responsible for geodetic solutions in the country?

Kind regards,
Maryam

Jochem · 31 mei 2023 om 08:20

I am from the responsible agency for the national projected coordinate system (RD). We are working together with other agencies in the Netherlands Partnership Geodetic Infrastructure (NSGI). Information on the RD parameters is available at the NSGI website: nsgi.nl/rdnaptrans The documentation includes PROJ strings that can be used by pyproj.

Based on the size of your errors, I suspect you forgot to do the datum transformation: EPSG transformation code: 9282

mariamsajadian · 31 mei 2023 om 09:44

Hi Jochem,
Thanks for your message. My colleague is suspect to a different datum parameter set, or some different additional transformation, such as possibly an offset matrix.
I have a control point projected into WGS using python code
p1 = 142636.01, 411623.18
P1 _WGS 84 = 5.207958042165996 51.693304297668696
My co-worker tried to update sdo_cs.transform function using below configurations.
Here is his email: "

I have been trying several configurations. None of them replicate what you are expecting, at this point:

– WKT according to current epsg.io:
PROJCS[“Amersfoort / RD New”, GEOGCS [ “Amersfoort”, DATUM [“Amersfoort (EPSG ID
5006289)”, SPHEROID [“Bessel 1841 (EPSG ID 7004)”, 6377397.155, 299.1528128], 5
65.4171, 50.3319, 465.5524, 1.9342, -1.6677, 9.1019, 4.0725], PRIMEM [ “Greenwic
h”, 0.000000000 ], UNIT [“Decimal Degree”, 0.0174532925199433]], PROJECTION [“Ob
lique Stereographic”], PARAMETER [“Latitude_Of_Origin”, 52.1561605555555556], PA
RAMETER [“Central_Meridian”, 5.3876388888888889], PARAMETER [“Scale_Factor”, 0.9
999079], PARAMETER [“False_Easting”, 155000.0], PARAMETER [“False_Northing”, 463
000.0], UNIT [“Meter”, 1.0]]

SQL> select
2 sdo_cs.transform(
3 sdo_geometry(
4 2001,
5 5028992,
6 sdo_point_type(142636.01, 411623.18, null),
7 null,
8 null),
9 4326)
10 from
11 dual;

SDO_CS.TRANSFORM(SDO_GEOMETRY(2001,5028992,SDO_POINT_TYPE(142636.01,411623.18,NU

SDO_GEOMETRY(2001, 4326, SDO_POINT_TYPE(5.21068153, 51.6938635, NULL), NULL, NULL)
– Compared to 1)“POINT(5.207958042165996 51.693304297668696)”, this is ~270m off

– WKT according to your email:
PROJCS[“Amersfoort / RD New”, GEOGCS [ “Amersfoort”, DATUM [“Amersfoort (EPSG ID
5006289)”, SPHEROID [“Bessel 1841 (EPSG ID 7004)”, 6377397.155, 299.1528128]],
PRIMEM [ “Greenwich”, 0.000000000 ], UNIT [“Decimal Degree”, 0.0174532925199433]
], PROJECTION [“Oblique Stereographic”], PARAMETER [“Latitude_Of_Origin”, 52.156
1605555555556], PARAMETER [“Central_Meridian”, 5.3876388888888889], PARAMETER ["
Scale_Factor", 0.9999079], PARAMETER [“False_Easting”, 155000.0], PARAMETER [“Fa
lse_Northing”, 463000.0], UNIT [“Meter”, 1.0]]

SQL> select
2 sdo_cs.transform(
3 sdo_geometry(
4 2001,
5 5028992,
6 sdo_point_type(142636.01, 411623.18, null),
7 null,
8 null),
9 4326)
10 from
11 dual;

SDO_CS.TRANSFORM(SDO_GEOMETRY(2001,5028992,SDO_POINT_TYPE(142636.01,411623.18,NU

SDO_GEOMETRY(2001, 4326, SDO_POINT_TYPE(5.20879089, 51.694765, NULL), NULL, NULL)
– Compared to 1)“POINT(5.207958042165996 51.693304297668696)”, this is ~150m off

– Tfm based on current Oracle Spatial WKT:
SQL> select
2 sdo_cs.transform(
3 sdo_geometry(
4 2001,
5 28992,
6 sdo_point_type(142636.01, 411623.18, null),
7 null,
8 null),
9 4326)
10 from
11 dual;

SDO_CS.TRANSFORM(SDO_GEOMETRY(2001,28992,SDO_POINT_TYPE(142636.01,411623.18,NULL

SDO_GEOMETRY(2001, 4326, SDO_POINT_TYPE(5.20837986, 51.6932492, NULL), NULL, NULL)
– Compared to 1)“POINT(5.207958042165996 51.693304297668696)”, this is ~40m off

Expected Tfm:
Wkt in python:

1)“POINT(5.207958042165996 51.693304297668696)”

So, I suspect that the source you compare with may use either a different datum parameter set, or some different additional transformation, such as possibly an offset matrix. Do you know whom to ask to find out whether they truly use the following WKT:

PROJCRS[“Amersfoort / RD New”,
BASEGEOGCRS[“Amersfoort”,
DATUM[“Amersfoort”,
ELLIPSOID[“Bessel 1841”,6377397.155,299.1528128,
LENGTHUNIT[“metre”,1]]],
PRIMEM[“Greenwich”,0,
ANGLEUNIT[“degree”,0.0174532925199433]],
ID[“EPSG”,4289]],
CONVERSION[“RD New”,
METHOD[“Oblique Stereographic”,
ID[“EPSG”,9809]],
PARAMETER[“Latitude of natural origin”,52.1561605555556,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8801]],
PARAMETER[“Longitude of natural origin”,5.38763888888889,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8802]],
PARAMETER[“Scale factor at natural origin”,0.9999079,
SCALEUNIT[“unity”,1],
ID[“EPSG”,8805]],
PARAMETER[“False easting”,155000,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8806]],
PARAMETER[“False northing”,463000,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8807]]],
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[“Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone.”],
BBOX[50.75,3.2,53.7,7.22]],
ID[“EPSG”,28992]]
Proj4
+proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m
"

Are we missing something?

Kind regards,
Maryam

Jochem · 31 mei 2023 om 10:28

Before I respond to your question, I have to repeat to NEVER use epsg.io! It is an unofficial site that has serious bugs and uses an outdated copy of the database from the official epsg.org.

Thanks for the additional information. This confirms my suspicion that problem is that you do only the projection and not the datum transformation (although the first WKT string does specify outdated approximate datum transformation parameters).

If you use (pyproj with) a uptodate PROJ version and an uptodate EPSG database version you can get the correct WKT string with the projinfo command:

\>projinfo -s epsg:4937 -t epsg:28992``

...

CONCATENATEDOPERATION["Inverse of Amersfoort to ETRS89 (9) + RD New",
    SOURCECRS[
        GEOGCRS["ETRS89",
            ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
                MEMBER["European Terrestrial Reference Frame 1989"],
                MEMBER["European Terrestrial Reference Frame 1990"],
                MEMBER["European Terrestrial Reference Frame 1991"],
                MEMBER["European Terrestrial Reference Frame 1992"],
                MEMBER["European Terrestrial Reference Frame 1993"],
                MEMBER["European Terrestrial Reference Frame 1994"],
                MEMBER["European Terrestrial Reference Frame 1996"],
                MEMBER["European Terrestrial Reference Frame 1997"],
                MEMBER["European Terrestrial Reference Frame 2000"],
                MEMBER["European Terrestrial Reference Frame 2005"],
                MEMBER["European Terrestrial Reference Frame 2014"],
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]],
                ENSEMBLEACCURACY[0.1]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,3],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["ellipsoidal height (h)",up,
                    ORDER[3],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",4937]]],
    TARGETCRS[
        PROJCRS["Amersfoort / RD New",
            BASEGEOGCRS["Amersfoort",
                DATUM["Amersfoort",
                    ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                ID["EPSG",4289]],
            CONVERSION["RD New",
                METHOD["Oblique Stereographic",
                    ID["EPSG",9809]],
                PARAMETER["Latitude of natural origin",52.1561605555556,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8801]],
                PARAMETER["Longitude of natural origin",5.38763888888889,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8802]],
                PARAMETER["Scale factor at natural origin",0.9999079,
                    SCALEUNIT["unity",1],
                    ID["EPSG",8805]],
                PARAMETER["False easting",155000,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8806]],
                PARAMETER["False northing",463000,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8807]]],
            CS[Cartesian,2],
                AXIS["easting (X)",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1]],
                AXIS["northing (Y)",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",28992]]],
    STEP[
        COORDINATEOPERATION["Inverse of Amersfoort to ETRS89 (9)",
            SOURCECRS[
                GEOGCRS["ETRS89",
                    ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
                        MEMBER["European Terrestrial Reference Frame 1989"],
                        MEMBER["European Terrestrial Reference Frame 1990"],
                        MEMBER["European Terrestrial Reference Frame 1991"],
                        MEMBER["European Terrestrial Reference Frame 1992"],
                        MEMBER["European Terrestrial Reference Frame 1993"],
                        MEMBER["European Terrestrial Reference Frame 1994"],
                        MEMBER["European Terrestrial Reference Frame 1996"],
                        MEMBER["European Terrestrial Reference Frame 1997"],
                        MEMBER["European Terrestrial Reference Frame 2000"],
                        MEMBER["European Terrestrial Reference Frame 2005"],
                        MEMBER["European Terrestrial Reference Frame 2014"],
                        ELLIPSOID["GRS 1980",6378137,298.257222101,
                            LENGTHUNIT["metre",1]],
                        ENSEMBLEACCURACY[0.1]],
                    PRIMEM["Greenwich",0,
                        ANGLEUNIT["degree",0.0174532925199433]],
                    CS[ellipsoidal,3],
                        AXIS["geodetic latitude (Lat)",north,
                            ORDER[1],
                            ANGLEUNIT["degree",0.0174532925199433]],
                        AXIS["geodetic longitude (Lon)",east,
                            ORDER[2],
                            ANGLEUNIT["degree",0.0174532925199433]],
                        AXIS["ellipsoidal height (h)",up,
                            ORDER[3],
                            LENGTHUNIT["metre",1]],
                    ID["EPSG",4937]]],
            TARGETCRS[
                GEOGCRS["Amersfoort",
                    DATUM["Amersfoort",
                        ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
                            LENGTHUNIT["metre",1]]],
                    PRIMEM["Greenwich",0,
                        ANGLEUNIT["degree",0.0174532925199433]],
                    CS[ellipsoidal,2],
                        AXIS["geodetic latitude (Lat)",north,
                            ORDER[1],
                            ANGLEUNIT["degree",0.0174532925199433]],
                        AXIS["geodetic longitude (Lon)",east,
                            ORDER[2],
                            ANGLEUNIT["degree",0.0174532925199433]],
                    ID["EPSG",4289]]],
            METHOD["Inverse of HORIZONTAL_SHIFT_GTIFF"],
            PARAMETERFILE["Latitude and longitude difference file","nl_nsgi_rdtrans2018.tif"],
            OPERATIONACCURACY[0.001],
            ID["INVERSE(DERIVED_FROM(EPSG))",9282],
            REMARK["Replaces Amersfoort to ETRS89 (7) (tfm code 7000). Horizontal component of official transformation RDNAPTRANS(TM)2018."]]],
    STEP[
        CONVERSION["RD New",
            METHOD["Oblique Stereographic",
                ID["EPSG",9809]],
            PARAMETER["Latitude of natural origin",52.1561605555556,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",5.38763888888889,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8802]],
            PARAMETER["Scale factor at natural origin",0.9999079,
                SCALEUNIT["unity",1],
                ID["EPSG",8805]],
            PARAMETER["False easting",155000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8806]],
            PARAMETER["False northing",463000,
                LENGTHUNIT["metre",1],
                ID["EPSG",8807]],
            ID["EPSG",19914]]],
    USAGE[
        SCOPE["unknown"],
        AREA["Netherlands - onshore, including Waddenzee, Dutch Wadden Islands and 12-mile offshore coastal zone."],
        BBOX[50.75,3.2,53.7,7.22]]]

...

The grid file nl_nsgi_rdtrans2018.tif that you need, is included (without the prefix nl_nsgi_ in the name) in the RDNAPTRANS download you already requested.

If you need physical heights in the national height system you need to use EPSG:7415 instead of EPSG:28992.

mariamsajadian · 31 mei 2023 om 12:41

Thanks! it was very useful.
what is projinfo ? is it a library?

Jochem · 31 mei 2023 om 12:57

projinfo is an executable of the PROJ library (and pyproj is an python interface for PROJ). You can use projinfo in the OSGeo4W Shell of QGIS.

mariamsajadian · 31 mei 2023 om 14:58

Hi Jochem,

My co-worker find the zip file very helpful.
He asked if you have ASCII version of rdtrans2018.gsb.

Here is his email: "

Your previous email was very helpful, with the forwarded info from the forum and the zip file.

The discrepancy happens, since your expected result assumes an NTv2 transformation (not merely a 7-parameter datum transformation), and we do support NTv2. I can help you with setting up this NTv2 transformation.

While the email includes the rdtrans2018.gsb, we actually need the ASCII version of that file, which may typically be called .gsa or .txt. But the included rdtrans2018.txt is not quite the correct format. The format we are looking for (ASCII version of .gsb) would have a header with approximately 22 lines, before the data.

Do you have access to this ASCII version of rdtrans2018.gsb? "

Kind regards,
Maryam

Jochem · 31 mei 2023 om 15:15

I think that I might have the type of ascii file your coworker mentions. I have made it temporarily available here: https://gnss-data.kadaster.nl/misc/rdnaptrans2018/rdtrans2018.asc

mariamsajadian · 31 mei 2023 om 15:18

Thanks a lot!
Kind regards,
Maryam

system · 27 november 2023 om 15:18

Dit topic is 180 dagen na het laatste antwoord automatisch gesloten. Nieuwe antwoorden zijn niet meer toegestaan.