[???] RDNAPTRANS2018 berekening ETRS89->RD

[jvdmeer]

Zoals te zien is, heb ik de taal even leeggelaten in de topictitel, want die is niet heel erg relevant in dit geval.
Op de website van NSGI (https://www.nsgi.nl/rdnaptrans) zijn de berekeningen te downloaden om om te rekenen vanaf ETRS89 naar RD.

Ik heb hier eerst geprobeerd een PHP-implementatie van variant 1 te maken. En die lijkt te werken op het probleem na dat de uitkomsten niet kloppen.

Daarna een hoop rondgezocht op internet en een implementatie is SAS gevonden. Deze heb ik toen naar PHP vertaald. En nogmaals geprobeerd en ik kwam weer op verkeerde uitkomsten uit. Nu bestaat de berekening uit meerdere stappen, waarbij ik dus niet weet in welke stap ik de fout moet zoeken.

ik ben nu dus op zoek naar iemand die me de verschillende tussenresultaten kan leveren bij een ETRS89-coördinaat.

SAS-implementatie
De gevolgde SAS-implementatie is:
https://github.com/FVelli...ain/gm_rdnaptrans2018.sas

Die schijnt correct te zijn, maar ik heb geen SAS tot mijn beschikking om daar tussenresultaten uit te halen.

Mijn controle-coördinaat
Voor mijn controle gebruik ik de eerste coördinaat in het bestand: Z001_ETRS89andRDNAP.txt
https://github.com/FVelli...n/Z001_ETRS89andRDNAP.txt

ETRS89:
lat: 51.728601274
lon: 4.712120126
hoogte: 301.7981

Verwacht resultaat:
x:108360.8790
y:415757.2745
h:258.0057

Mijn resultaat:
'gmv_RD_x_lon' => float 108360.88437748
'gmv_RD_y_lat' => float 415757.27429712
'gmv_NAP_height' => float 258.01342490227

Het gaat maar om kleine verschillen, maar de afwijking is te groot om geaccepteerd te worden als correct. Het verwachte resultaat is wat geaccepteerd wordt.

Benodigd
De resultaten van de verschillende stappen die in de pdf worden genoemd.
ik heb de namen van de waarden genomen zoals die genoemd worden in de SAS-implementatie:

2.1 Omrekenen naar radialen
Mijn resultaat:
'gmv_ETRS89_lat_rad' => float 0.90283440968263
'gmv_ETRS89_lon_rad' => float 0.08224201094819

2.2 Datum transformatie (stap 1)
'gmv_geoc_cart_ETRS89_X' => float 3945358.2222048
'gmv_geoc_cart_ETRS89_Y' => float 325207.73269047
'gmv_geoc_cart_ETRS89_Z' => float 4984189.6144722

2.2 Helbert-transformatie naar pseudo Bessel (stap 2)
'gmv_geoc_cart_ETRS89_X' => float 3945358.2222048
'gmv_geoc_cart_ETRS89_Y' => float 325207.73269047
'gmv_geoc_cart_ETRS89_Z' => float 4984189.6144722

'gmv_geog_ellips_psdo_Bessel_lat' => float 0.90285076437824
'gmv_geog_ellips_psdo_Bessel_lon' => float 0.082247920070159

2.3 Iteratieve RDcorrectie
'gmv_geog_ellips_real_Bessel_lat' => float 51.72953802384
'gmv_geog_ellips_real_Bessel_lon' => float 4.7124591666284

2.4 Omzetting naar RD
'gmv_RD_x_lon' => float 108360.88437748
'gmv_RD_y_lat' => float 415757.27429712

2.5 Hoogte-correctie
'gmv_NAP_height' => float 258.01342490227

Waarbij dus de laatste drie fout zijn.

Wie kan mij op weg helpen met de juiste tussenresultaten bij de genoemde ETRS89-coördinaten?
Zodat ik specifiek kan zoeken waar de fout zit.

Voor de compleetheid, de complete lijst met variabelen zoals die uit mijn PHP-versie van het SAS-script komen:
De waarden in de correctioncache, zijn diegene die zijn opgezocht in een mysql-db voor stap 2.3 en 2.5.

code:

'correctionCache' => 
    array (size=4)
      51.7375,4.7200 => 
        array (size=3)
          'LatCor' => string '0.000000320' (length=11)
          'LonCor' => string '-0.000000440' (length=12)
          'Height' => string '43.7790' (length=7)
      51.7375,4.7400 => 
        array (size=3)
          'LatCor' => string '0.000000320' (length=11)
          'LonCor' => string '-0.000000520' (length=12)
          'Height' => string '43.7699' (length=7)
      51.7250,4.7200 => 
        array (size=3)
          'LatCor' => string '0.000000298' (length=11)
          'LonCor' => string '-0.000000414' (length=12)
          'Height' => string '43.7949' (length=7)
      51.7250,4.7400 => 
        array (size=3)
          'LatCor' => string '0.000000296' (length=11)
          'LonCor' => string '-0.000000493' (length=12)
          'Height' => string '43.7855' (length=7)
  'gmv_GRS80_a' => int 6378137
  'gmv_GRS80_f' => float 0.0033528106811823
  'gmv_ETRS89_h0' => int 43
  'gmv_epsilon_GRS80_threshold' => float 2.0E-11
  'gmv_tX' => float -565.7346
  'gmv_tY' => float -50.4058
  'gmv_tZ' => float -465.2895
  'gmv_D_dec' => float 5.3876388888889
  'gmv_D_rad' => float 0.0940320375196
  'gmv_alpha' => float -1.9151255475293E-6
  'gmv_beta' => float 1.6036473018269E-6
  'gmv_gamma' => float -9.0954585716021E-6
  'gmv_delta' => float -4.07242E-6
  'gmv_tX_inv' => float 565.7381
  'gmv_tY_inv' => float 50.4018
  'gmv_tZ_inv' => float 465.2904
  'gmv_alpha_inv' => float 1.9151400919398E-6
  'gmv_beta_inv' => float -1.6036279092796E-6
  'gmv_gamma_inv' => float 9.0954634197389E-6
  'gmv_delta_inv' => float 4.07244E-6
  'gmv_B1841_a' => float 6377397.155
  'gmv_B1841_f' => float 0.0033427731821748
  'gmv_epsilon_B1841_threshold' => float 2.0E-11
  'gmv_RD_Bessel_h0' => int 0
  'gmv_phi_min' => int 50
  'gmv_phi_max' => int 56
  'gmv_labda_min' => int 2
  'gmv_labda_max' => int 8
  'gmv_phi_delta' => float 0.0125
  'gmv_labda_delta' => float 0.02
  'gmv_c0' => int 0
  'gmv_epsilon_RD_threshold' => float 1.0E-9
  'gmv_k_amersfoort' => float 0.9999079
  'gmv_x0_amersfoort' => int 155000
  'gmv_y0_amersfoort' => int 463000
  'gmv_phi0_amersfoort' => float 0.91029672689324
  'gmv_labda0_amersfoort' => float 0.0940320375196
  'lmv_height_bln' => int 1
  'gmv_ETRS89_lat_dec' => float 51.728601274
  'gmv_ETRS89_lon_dec' => float 4.712120126
  'gmv_ETRS89_height' => float 301.7981
  'gmv_ETRS89_lat_rad' => float 0.90283440968263
  'gmv_ETRS89_lon_rad' => float 0.08224201094819
  'gmv_geoc_cart_ETRS89_X' => float 3945358.2222048
  'gmv_geoc_cart_ETRS89_Y' => float 325207.73269047
  'gmv_geoc_cart_ETRS89_Z' => float 4984189.6144722
  'gmv_geoc_cart_RD_Bessel_X' => float 3944765.4696156
  'gmv_geoc_cart_RD_Bessel_Y' => float 325182.34180738
  'gmv_geoc_cart_RD_Bessel_Z' => float 4983710.9769915
  'gmv_geog_ellips_psdo_Bessel_lat' => float 0.90285076437824
  'gmv_geog_ellips_psdo_Bessel_lon' => float 0.082247920070159
  'lmv_phi0' => float 51.729538329034
  'lmv_phi' => float 51.729538329034
  'lmv_phi1' => float 51.72953802384
  'lmv_labda0' => float 4.7124586937494
  'lmv_labda' => float 4.7124586937494
  'lmv_labda1' => float 4.7124591666284
  'gmv_geog_ellips_real_Bessel_lat' => float 51.72953802384
  'gmv_geog_ellips_real_Bessel_lon' => float 4.7124591666284
  'gmv_RD_x_lon' => float 108360.88437748
  'gmv_RD_y_lat' => float 415757.27429712
  'gmv_NAP_height' => float 258.01342490227

[ Voor 7% gewijzigd door jvdmeer op 06-08-2022 00:20 . Reden: Zag een copy-paste fout ]

[downtime]

Misschien zeg ik nu iets doms maar is de afwijking niet verwacht als je berekeningen met floats vertaalt van de ene math library naar de ander? Oftewel: Wie zegt dat de berekening in PHP exact dezelfde uitkomst als in SAS moet hebben?

[eheijnen]

Misschien heb je iets aan de excel sheet op deze pagina:
https://titans.ladesk.com...TRS89-conversie-met-Excel

[RayNbow]

eheijnen schreef op zaterdag 6 augustus 2022 @ 10:07:
Misschien heb je iets aan de excel sheet op deze pagina:
https://titans.ladesk.com...TRS89-conversie-met-Excel

Zo'n Excel sheet zal in ieder geval niet over RDNAPTRANS2018 gaan:

Updated: Jan 14, 2017

---

@jvdmeer is er een reden waarom je niet gebruikmaakt van PROJ?

Edit: En heb je de validatieservice al geprobeerd?

[ Voor 16% gewijzigd door RayNbow op 06-08-2022 10:26 ]

[jvdmeer]

Die Excelsheet is inderdaad net te oud, maar slaat ook stap 2.3 over. De correctie op basis van het correctiegrid.

En PROJ had ik ook gevonden, maar ik heb bij mijn hoster en probleem met het aanroepen van uitvoerbare bestanden. En zuiver ik kon zien, kon ik proj niet eenvoudig integreren in PHP.
Daarom toen de keus gemaakt om de berekening zelf te implementeren.
In principe is de precisie voor normaal gebruik al ruim voldoende. Maar het lijkt me leuk als ik dat stempeltje rdnaptrans2018 er op kan plakken.

Maar ik zal ook nog eens kijken naar proj4php: https://github.com/proj4php/proj4php

De validatieservice heb ik nog niet geprobeerd, omdat ik nu al faal op het eerste punt van de zelfvalidatie.

[ Voor 8% gewijzigd door jvdmeer op 06-08-2022 11:12 ]

[jvdmeer]

downtime schreef op zaterdag 6 augustus 2022 @ 09:49:
Misschien zeg ik nu iets doms maar is de afwijking niet verwacht als je berekeningen met floats vertaalt van de ene math library naar de ander? Oftewel: Wie zegt dat de berekening in PHP exact dezelfde uitkomst als in SAS moet hebben?

Vanuit de eisen van rdnaptrans2018 wordt verwacht dat je in ieder geval tot het 4e decimaal de correcte uitkomsten geeft. Het hoeft dus wat mij betreft niet exact hetzelfde te zijn, maar als ik kan zien wat de tussenresultaten zijn, kan ik specifieker zoeken waar het fout gaat.
Ik kan dan ook eventueel een tussenresultaat injecteren halverwege het proces, om te controleren of het vervolg goed gaat.

[MadGazelle]

Ik denk dat het de moeite waard kan zijn om je topic te crossposten naar het RDNAPTRANS subforum op geoforum.nl. Daar zijn namelijk mensen van het NSGI actief. Bijvoorbeeld zie dit topic van de implementatie van RDNAPTRANS2018 in Java.

[Umbrah]

In plaats van het zelf te re-coden, zou je evt. een externe library aan kunnen roepen. GDAL is ongeveer goud qua FOSS geodetische zaken:

https://gdal.org/programs/gdaltransform.html

RDNAP is geen formele EPSG richtlijn, maar desondanks wel herkenbaar, dus dat zal geen probleem moeten opleveren.

Ik ben absoluut geen specialist in PHP (en ik zal ook zeggen, eigenlijk ook niet in GDAL in alle detail -- ik houd me hoofdzakelijk bezig met ESRI/FME-achtigen), maar er zijn wrappers die GDAL + PHP koppelen: https://github.com/geo6/php-gdal-wrapper

Alternatief: is het iets wat je client-side kan doen? Wederom: ben meer een esri-aan dan een OSGeo-er, maar in javascript zijn er zat libraries die je in staat stellen om zoiets toe te passen: https://developers.arcgis...-geometry-projection.html

[Janoz]

jvdmeer schreef op zaterdag 6 augustus 2022 @ 00:08:

Verwacht resultaat:
x:108360.8790
y:415757.2745
h:258.0057

Mijn resultaat:
'gmv_RD_x_lon' => float 108360.88437748
'gmv_RD_y_lat' => float 415757.27429712
'gmv_NAP_height' => float 258.01342490227

Het gaat maar om kleine verschillen, maar de afwijking is te groot om geaccepteerd te worden als correct. Het verwachte resultaat is wat geaccepteerd wordt.

Als ik even snel kijk in naar de significantie, dan lijkt het er heel erg op dat je berekening gebruik maakt van 32bit floats. Bij 32 bit floats wordt 1 bit gebruikt voor sign, 8 voor exponent en 23 voor de mantisse (eigenlijk 24 omdat de eerste bit altijd 1 is en daarom weggelaten kan worden). 10836088 is precies wat nog in die 24 bits past. De verschillen die jij ziet zijn dus gewoon afrondingsfouten. Je berekening is waarschijnlijk correct.

Om het werkend te krijgen zul je php dos op 1 of andere manier zover moeten krijgen dat er 64bit floats gebruikt worden. De php manual is mij niet helemaal duidelijk. Daar wordt aangegeven dat het platform afhankelijk is en geven ze aan dat het meestal 64bit is (wat blijkbaar hier niet het geval is).

Je zou kunnen kijken of BCMath helpt.

[jvdmeer]

Ik heb me voor dit probleem ook maar aangemeld op geoforum.nl. En via Jochem heb ik de fout in mijn berekening gevonden. Deze zat in stap 2.3 en 2.5.

Zie ook de discussie aldaar:
https://geoforum.nl/t/kle...dnaptrans2018-in-php/7443