Kleine afwijking bij RDNAPTRANS2018 in PHP

Standaarden RDNAPTRANS
#1

Voor een website probeer ik een PHP-implementatie te maken van RDNAPTRANS2018. Mijn eerste versie had redelijk grote afwijkingen. Omdat het moeilijk is om te constateren in welke stap de afwijking is ontstaan, heb ik daarna de SAS-implementatie die ik heb gevonden op https://github.com/FVelli…ain/gm_rdnaptrans2018.sas herschreven naar PHP.

Die resultaten waren al een stuk beter, maar ik constateer nog steeds een klein verschil in de uitkomst. ik hoop dat iemand me kan duiden in welke stap het verschil optreedt.

Ik ben even uitgegaan van het eerste punt in het zelfvalidatie-bestand:
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

Tussenreultaten:
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

Voor de compleetheid, heb ik hieronder een compleet overzicht van alle variabelen die voorbij komen in mijn PHP-implementatie van de SAS-variant:
‘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

Ik hoop dat iemand me kan helpen en kan aangeven welke tussenresultaten kloppen en waar de afwijking begint.

Oh ja, de gebruikte waarden van het correctiegrid voor dit ene punt:
‘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)

Open source RDNAPTRANS2018 Java implementatie
#2

Het lijkt mis te gaan bij de “2.3 Iteratieve RD-correctie”, daar verwacht ik:

lat = 51.7295380258
lon =  4.7124590893

En ook bij “2.5 Hoogte-correctie” gaat het mis.

Beiden zouden verklaard kunnen worden door een eventuele fout in de bilineaire interpolatie van de correctiegrids of het op het uitlezen van gridwaarden op de verkeerde plek in het grid. Ik vermoed dat laatste.

#3

Bedankt Jochem. Je hebt me volgens mij weer een stuk verder geholpen. Het ging inderdaad fout met de lengtegraad. Die was 4.71… en het rdcorrectiegrid werd genomen voor 4.72 en 4.74 dus 1 vakje teveel naar het oosten.

De waarden zijn nu al veel beter.
‘gmv_RD_x_lon’ => float 108360.87903706 => na afronding in orde
‘gmv_RD_y_lat’ => float 415757.27455575 => na afronding 0.0001 te hoog
‘gmv_NAP_height’ => float 258.0043083297 => na afronding 0.0014 te laag

Van de y weet ik niet of dat binnen de marges valt…
Het lijkt dus nu alleen nog een probleem met de hoogte.

#4

De xy-coördinaten zijn voor dit punt nu goed. De marge om toestemming aan te kunnen vragen voor het gebruik van de merknaam RDNAPTRANS is 0.0010 m. Voor die aanvraag (laatste bladzijde van de pdf van de RDNAPTRANS-download) moeten echter alle 10000 punten van de validatieservice en ook de hoogte goed gaan. Is het probleem bij de hoogte niet ook dat je gridwaarden op de verkeerde plek in het grid uitleest? Welk gridbestand gebruik je?

#5

Naar aanleiding van je opmerking heb ik nu de hoogte-waarden maar nagekeken in de mysql-db en daar was dus iets fout gegaan met de import. Zojuist de import opnieuw uitgevoerd en nu klopt de hoogte ook in mijn PHP-script op basis van SAS:
‘gmv_RD_x_lon’ => float 108360.87903706
‘gmv_RD_y_lat’ => float 415757.27455575
‘gmv_NAP_height’ => float 258.00566518823

Nu dit script nog testen met de rest van de zelfvalidatie. En daarna mijn oorspronkelijke script ombouwen op basis van dit SAS-PHP-script. Hierdoor zijn de diverse onderdelen als de helmert-transformatie en de bilinaire interpolatie beter herbruikbaar.
Als dan de zelfvalidatie nog steeds slaagt, ga ik me aan de validatie-service wagen.
Maar ik heb nu in ieder geval een werkend script!.

5 likes
#6

Nu ik klaar ben met mijn PHP-script ombouwen, kom ik eindelijk grotendeels door de zelfvalidatie heen. Alleen wat wordt er verwacht bij de regels waar de NAP-hoogte een asterix () is?
Bijvoorbeeld het laatste punt:
ID:30019999
ETRS-Lat:40.961335618
ETRS-Lon:-1.461398258
ETRS-H:236.3071
RD-x:-426598.9295
RD-Y:-759790.1238
RD-h:

Het zijn zo te zien alle punten die buiten het correctiegrid liggen

#7

Mooi gemaakt! Wordt de code beschikbaar onder een open source licentie?

3 likes
#8

Buiten het grid voor de hoogtetransformatie is NAP ongedefinieerd Het beste is daarom om geen resultaat te geven in de output.

De pdf van de RDNAPTRANS-download zegt daarover:

A no result outside the bounds of the quasi-geoid grid model is best accompanied with a warning that the NAP coordinates are out of bounds.

Als je code toch een waarde terug moet geven dan is kun je het beste NaN gebruiken of zelf een no-data value zoals -9999.9999 definieren. Vervolgens kan een front-end op je code bij deze waarde dan de output leeg laten en bovenstaande waarschuwing geven.

Bij toepassingen waar ook buiten het grid toch echt een benadering van de NAP-hoogte nodig is zou de EGM-hoogte gebruikt kunnen worden met een waarschuwing er bij dat de hoogte geen NAP is.

#9

Is het gelukt met de validatieservice 100% te krijgen?

#10

Nee, helaas nog niet. De score bij ETRS89->RD komt op 34% en andersom op maar 2%.
Het gaat echt om hele kleine afwijkingen, want bij de zelfvalidatie krijg ik alleen maar afwijkingen op het achterste decimaal. Bijvoorbeeld de 1e van zelfvalidatie die fout gaat is:
punt: 30010020
ETRS89: 53.4731, 6.8867 (h=290.4741)
RD juist: x=254565.6835, y=610693.7409, h=250.1468
RD bij mij: x=254565.6835, y=610693.7408, h=250.1468

Voor nu heb ik mijn prioriteiten even verlegd naar een ander gedeelte van mijn site. Later ga ik nog verder zoeken waar dit hele kleine verschil vandaan komt.

V.w.b. open source. Is een idee, maar pas nadat ik het mysql-gedeelte heb vervangen door het inlezen van de meegeleverde tekstfiles.

#11

Dat de “laatste” decimaal (de 4e voor RD, NAP en ETRS89-hoogte in meters en de 9e voor ETRS89-lengte en -breedte in graden) een verschil geeft is normaal. Dit decimaal komt overeen met een tiende van een millimeter. Voor validatie van RDNAPTRANS is een afwijking van 0,0010 m of 0,000 000 010° toegestaan. Dit komt over een met een millimeter. Daarmee zou je resultaat dus gewoon door de Validatieservice moeten komen. Stuur anders je resultaten naar rd@kadaster.nl, dan kan ik kijken of ik een aanwijzing zie voor de oorzaak.

Let op, in de zelfvalidatieset zitten (in tegenstelling tot de Validatieservice) ook wat punten op honderden kilometers buiten de rechthoek om Nederland en EEZ, waar de transformatie sowieso niet meer op 1 millimeter nauwkeurig kan.

#12

De berekening ETRS89 naar RD lijkt nu voor een groot deel goed te gaan. Maar nu gaat het omgekeerde nog fout. Dus de omgekeerde stappen ben ik 1 voor 1 aan het controleren met behulp van de getallen die ik van de heenweg krijg.

nu zit bij die analyse met een vraag:
Als ik even naar de heenweg kijk richt ik me op de actie in paragraaf 2.2.1
:
ETRS to RD
Startpunt
‘lat’ => float 51.728601274
‘lon’ => float 4.712120126
‘h’ => float 301.7981

Stap 2.2.1
input Ellipsoidal geographic ETRS89 (na Helmert)
‘x’ => float 3944765.4696156
‘y’ => float 325182.34180738
‘z’ => float 4983710.9769915
Output:
‘lat’ => float 0.90285076437821
‘lon’ => float 0.082247920070159

Eindresultaat (RD en voldoet in de zelfvalidatie)
‘x’ => float 108360.87902938
‘y’ => float 415757.27450621
‘z’ => float 258.00566518823

Als ik dat dan vergelijk met de terugweg:
RD to ETRS (bevat fout)

Pak ik als RD-startpunt de output van de heenweg:
‘x’ => float 108360.87902938
‘y’ => float 415757.27450621
‘h’ => float 258.00566518823

Kom ik aan 1e deel van stap 3.3. De input voor die stap is nagenoeg gelijk aan het resultaat van 2.2.1 van de heenweg
‘lat’ => float 0.90285076437819
‘lon’ => float 0.082247920070158

Als ik deze da omreken naar Ellipsoidal geographic ETRS89 (voor Helmert)
Kom ik uit op:
‘x’ => float 3944766.1190569
‘y’ => float 325182.39534334
‘z’ => float 4983711.8029911

Dan is nu het verschil wel heel groot. Deze waarden moeten volgens mij nog door de Helmert-transformatie heen. En het eindresultaat lijkt dus nergens op.

Stap 3.3 is volgens mij de volgende berekening:
$e2=$f*(2-$f);

$Rn=$a/sqrt(1-$e2*pow(sin($phi),2));

$X=($Rn + $h)*cos($phi)*cos($labda);
$Y=($Rn + $h)*cos($phi)*sin($labda);
$Z=($Rn*(1-$e2) + $h)*sin($phi);

Waarbij a,f en h gebaseerd zijn op Bessel1841:
$a=6377397.155
$f=1/299.1528128
$h=0

Pas ik nu toch een verkeerde berekening toe in deze stap?

#13

Volgens mij is er niks mis met dit stuk van je berekening.

 
Met input RD:

x = 108360.8790 m
y = 415757.2745 m

 
krijg ik voor RD Bessel:

lat = 0.90285075908 rad
lon = 0.08224792697 rad
h   = 0.0000 m

 
voor pseudo-RD Bessel:

lat = 0.90285076438 rad (dat krijg jij ook)
lon = 0.08224792006 rad (dat krijg jij ook)
h   = 0.0000 m

X = 3945358.8717 m  
Y =  325207.7862 m
Z = 4984190.4405 m

 
en voor ETRS89:

X = 3944766.1191 m (dat krijg jij ook)
Y =  325182.3953 m (dat krijg jij ook)
Z = 4983711.8030 m (dat krijg jij ook)

lat = 51.7286012740 deg = 0.90283440968 rad
lon =  4.7121201255 deg = 0.08224201094 rad

 
Stukje achtergrond: Het is normaal dat het tussenresultaat XYZ voor ETRS89 bij de berekening de ene kant op anders is dan de andere kant op. Dit zijn 3D coördinaten (waarbij geen van de 3 assen verticaal is). Omdat we 2D coördinaten in deze berekening stoppen moet er een hoogte aangenomen worden. Hiervoor gebruiken we een vaste hoogte die een benadering is van een NAP-hoogte van 0 m (h = 0 m voor Bessel en h = 43 m voor ETRS89) in plaats van de werkelijke hoogte, zodat de berekening ook voor 2D coördinaten hetzelfde eindresultaat geeft. Die benaderde hoogte is de ene kant op een paar meter anders dan de ander kant op, maar dat heeft (tot een paar honderd kilometer buiten Nederland) geen significante invloed op de 2D coördinaten van het eindresultaat.

#14

Als het goed is moet ik dan alleen die Helmert transformatie doen. Daar kom ik uit op:
‘x’ => 3945358.8681625
‘y’ => 325207.79020563
‘z’ => 4984190.4394363

En vanaf daar weer terugrekenen naar ellipsiodaal op basis van GRS80:
‘lat’ => 0.90283440997204
‘lon’ => 0.082242012023223
hier wijk ik dus af van de antwoorden die je in de vorige post gaf.
Waardoor ik in graden dus ook anders uitkom:
‘lat’ => 51.728601290582
‘lon’ => 4.7121201875949

Zou je de waarden kunnen controleren/geven die ik zou moeten krijgen na de Helmert-transformatie?

Heel erg bedankt.

#15

Bovenstaande is overbodig. Ik zag dat ik de Helmert-transformatie had gedaan met de omgekeerde parameters van de ETRS2RD-transformatie. Maar de waardes wijken lichtjes af. Met de juiste waarden in de Helmert-Transformatie, krijg ik nu netjes de 100% in de validatie bij RD->ETRS89. YES!

Jammer dat ETRS89->RD nu nog geen 100% oplevert :frowning:

#16

Klopt. Als je de omgekeerde parameters zou willen gebruiken dan zou je ook de volgorde van het schalen, de 3 rotaties en de translatie om moeten keren, maar dat is ongebruikelijk.

Gefeliciteerd!

Stuur anders nogmaals je resultaten naar rd@kadaster.nl, dan kijk ik of ik zie wat er mis gaat.

Alternatieve BAG Python parser
#17

Zo, dankzij de hulp van Jochem een hele stomme fout verwijderd uit de inleesroutine voor de validatie. De berekening bleek al een tijdje goed te zijn.

Voor het gemak heb ik de 2 PHP-bestanden en de bestanden vanaf ngis geplaatst op github:

Moet nog wel even de bestanden misschien wat netter maken en verwijzingen aanbrengen naar de website van nsgi.

6 likes
#18

Ben ik toch nog een keer… Misschien een nieuw topic waard, maar ik ga voor nu toch de vraag eerst hier stellen.
Vanaf de ETRS89-coördinaten wil ik nog 1 stap verder naar ITRF2014. ook weer ene hoop ingelezen, maar ik kom er nog niet uit.
Volgens de discussie hier: RDNAPTRANS in Excel - #3 door Jochem - RDNAPTRANS - Geoforum

Bestaat de transformatie uit 3 stappen:

  • Eerst weer naar een geocentrische coördinatenop basis van GRS80.
    (Moet dat met uitgerekende hoogte of ook weer met de fictieve hoogte van 43m)
  • Daarna het tijdsverschil uitrekenen in jaren tussen de huidige datum en 1-1-2010 (Voor de voorbeeldberekening wordt er 2020-10-01 genomen als huidige datum)
  • Dan de juiste helmert parameters uitrekenen, met als basis de volgende waarden:
    array(‘tx’=>0.0547,‘ty’=>0.0522,‘tz’=>-0.0741,‘rx’=>deg_to_rad1(+0.001701/3600),‘ry’=>deg_to_rad1(+0.010290/3600),‘rz’=>deg_to_rad1((-0.016632)/3600),‘s’=>+0.00000000212);

En als snelheid/jaar: array(‘tx’=>0.00010,‘ty’=>0.00010,‘tz’=>-0.00190,‘rx’=>deg_to_rad1(+0.0000810/3600),‘ry’=>deg_to_rad1(+0.0004900/3600),‘rz’=>deg_to_rad1((-0.0007920)/3600),‘s’=>+0.000000000110);
Waar kan ik de bron vinden van deze waarden? Ik heb ze nu uit de excel-discussie en uit de sas-implementatie van wgs84toitrf2014

En na de helmert-transformatie, weer terug naar ITRF2014, ook weer m.b.v. GRS80.

En dan zou ik op het juiste antwoord moeten uitkomen (volgens de SAS-implementatie).
ETRS 51.999994778222, 4.9999922685, 42.98187019
WGS84 ≈ ITRS2014: 52, 5, 43

Ik heb allerlei manieren geprobeerd, maar ik kom op andere antwoorden
Uitgangspunt ETRS89
‘lat’ => float 51.999994778222
‘lon’ => float 4.9999922685
‘h’ => float 42.98187019

Geocentric:
‘x’ => float 3920013.6292309
‘y’ => float 342956.22040882
‘z’ => float 5002836.8721211

Helmert-parameters na toepassen epoch-correctie (10.75):
‘tx’ => float -0.145225
‘ty’ => float -0.147725
‘tz’ => float 3.724475
‘rx’ => float -7.768569583705E-7
‘ry’ => float -4.6995050568092E-6
‘rz’ => float 7.5959347040671E-6
‘s’ => float -2.177975E-7

Na Helmert:
‘x’ => float 3919988.5121045
‘y’ => float 342990.10530378
‘z’ => float 5002852.3928946

Ellipsiodal (radian):
‘lat’ => float 0.90757534651209
‘lon’ => float 0.087275462441856
‘h’ => float 41.6441

Zou ITRS2014 moeten zijn:
‘lat’ => float 52.000236945266
‘lon’ => float 5.0005156529709
‘h’ => float 41.6441

Is er weer iemand die me op weg kan helpen?

#19

Er staat nu een MIT (open source) license op!! Bedankt! :+1::heart:

#20

De originele bron voor de transformatie tussen ITRS en ETRS89 is EUREF Technical Note 1, maar dit is een nogal lastig leesbaar document. Ik zal je een conceptversie van mijn eigen beschrijving in de stijl van de RDNAPTRANS-pdf mailen.

Het gaat bij de eerste stap van ETRS89 geografisch naar geocentrisch al mis. Jij gebruikt de vaste hoogte van 43 m (zoals bij ETRS89 naar RD) terwijl je bij ETRS89 naar ITRS de werkelijke hoogte (in jouw voorbeeld 42.9819 m) moet gebruiken.

Met het berekenen van de 7 parameters gaat ook iets mis. Ik kom voor epoche 2020.75 uit op:

tx = -0.0558   m
ty = -0.0533   m
tz = +0.0945   m
rx = +0.002572 arcsec
ry = +0.015558 arcsec
rz = -0.025146 arcsec
s  = -0.00330  ppm

Hopelijk kom je hiermee weer een stap verder. Ik heb wel nog een meer fundamentele vraag: Waarom wil je de tijdsafhankelijke transformatie naar ITRS wil implementeren? In ITRS veranderen de coördinaten van alle punten in Nederland met 2,4 cm per jaar. Voor geo-databestanden is dat nogal onpraktisch. Daarom is er op Europees niveau afgesproken de ITRS-coördinaten van 1989 vast te houden, dat noemen we ETRS89. De transformatie tussen ITRF2014 naar ETRS89 is daardoor alleen bij specifieke toepassingen nodig, zie Handreiking CRS van Geonovum (paragraaf 3.3.2.2).

1 like