Puntenwolk CRS converteren naar RDnew+NAP

Toepassingen
#1

Beste GIS experts,

Ik heb een gedetailleerde puntenwolk - enkele honderden gigabytes aan LAS data die onze interne pipeline oplevert in EPSG:4978, deze dataset is visueel geinspecteerd in Cesium en de hoogte boven het aardoppervlak lijkt aannemelijk.

Nu wil onze klant deze gegevens graag in RDnew+NAP (EPSG:7415) maar tot dusver lijken mijn pogingen om het om te zetten niet succesvol. Mijn huidige aanpak maakt gebruik van PDAL om een projectie te maken naar EPSG:7415.

Als verificatie stap combineer ik de resulterende puntenwolk met gegevens van AHN3, waar de punten van mijn dataset tientallen meters boven de AHN data uitkomen.

Heeft iemand suggesties voor een werkwijze?

#2

EPSG:4978 is geocentrisch WGS84 (X,Y,Z alle drie in meters, dus niet in graden). Ik vraag me af of je data echt in EPSG:4978 is, aangezien een horizontaal vlak daarin in een hoek van ca. 52 graden met de Z-as staat en je de hoogtes boven het aardoppervlak dus niet zo makkelijk geïnspecteerd kunt hebben.

WGS84 heeft per definitie een onnauwkeurigheid van enkele meters, dus een benaderde transformatie is goed genoeg voor je toepassing? Hoe kom je aan WGS84-coördinaten, stand-alone GPS? Of heb je een of andere referentie of correctiesignaal gebruikt? In dat laatste geval zijn je coördinaten niet in WGS84 maar in ETRS89 of ITRS.

De beste manier om te transformeren is met de open source bibliotheek PROJ.org (wat jij waarschijnlijk indirect al gebruikt), waarvan command line applicaties (cs2cs en cct) met het programma QGIS meegeleverd wordt (OSGeo4W Shell). Daarmee kan je eenvoudig een benaderde transformatie doen, maar de exacte transformatie is ook mogelijk.

Bij de eenvoudige benaderde transformatie naar EPSG:7415 wordt de ellipsoïdische hoogte van WGS84 gelijkgesteld aan NAP. Dat gaat niet goed omdat het verschil tussen ellipsoïdische WGS84-hoogte en NAP vrij groot is (ca. 40 m in Noord-Groningen tot 46 m in Zuid-Limburg). Een simplistische oplossing is om 43 meter van al je hoogtes af te trekken. Kom je dan wel ongeveer goed uit? Als je iets nauwkeurigers wil, dan kan ik je uitleggen hoe dat met PROJ kan, maar dan moet ik eerst weten of je coördinaten echt in EPSG:4978 zijn en welke nauwkeurigheid je nodig hebt voor de transformatie naar RDNAP.

3 likes
#3

De brongegevens komen van een Riegl laserscanner met een GNSS-RTK ontvanger. Daarna wordt expliciet aangegeven om de punten in WGS84 geocentric te rapporteren, wat helpt om de punten in cesium te brengen (specificatie)

Wij gebruiken voor de transformatie op dit moment de open source bibliotheek PDAL, die voor zover ik kan zien inderdaad intern gebruik maakt van PROJ.

Een afwijking van 43 meter klinkt wel ongeveer goed, ik ga dat in ieder geval proberen maar ik ben zeker ook geinteresseerd in de nauwkeurige aanpak.

#4

Als je RTK gebruikt dan zijn je coördinaten in het coördinatensysteem van je referentiestation. De coördinaten van een referentiestation zijn normaal gesproken niet in WGS84. Welk referentiestation of welke correctiedienst heb je gebruikt?

#5

We maken gebruik van Leica Smartnet, de ruwe data die we daar vandaan krijgen is in een lokaal stelsel ten opzichte van de scanner en wordt in een nabewerkingsstap mbv RiScan geogerefereerd.

Ik zal navragen wat de precieze instellingen zijn, je hebt een goed punt dat hier wellicht een verkeerde waarde is gebruikt.

#6

Waarom zou er een verkeerde waarde zijn gebruikt? Hij staat toch gewoon netjes in ECEF en WGS84?

#7

Ik denk ook dat er waarschijnlijk geen verkeerde instelling gebruikt is maar dat de RTK-coördinaten van Leica Smartnet die je krijgt in werkelijkheid in ETRS89 zijn, ondanks dat er in er in de data misschien staat dat het WGS84 is.

In dat geval is de juiste manier van transformeren:

input
   |
EPSG-code corrigeren van EPSG:4978 naar EPSG:4936
   |
ETRS89 X,Y,Z (EPSG:4936)
   |
   | conversie
   |
ETRS89 lon,lat,h (EPSG:4937)
   |
   | RDNAPTRANS™
   |
RD x,y (EPSG:28992); NAP H (EPSG:5709)
   |
   | combineren
   |
RDNAP x,y,H (EPSG:7415)

Met de nieuwste versies van PROJ gaat de transformatie van ETRS89 naar RDNAP automatisch op de millimeter goed, maar gezien de hoogtefout die je beschrijft vermoed ik dat je een oude PROJ-versie gebruikt. Welke PROJ-versie gebruik je?

In dit topic heb ik de PROJ4-string voor de nauwkeurige transformatie RDNAPTRANS™2018 gegeven:

Als dat niet werkt in jouw implementatie (bijvoorbeeld door het ontbreken van gridbestanden) dan is het ook nog wel mogelijk de hoogte beter te krijgen door de horizontale coördinaten en de hoogte afzonderlijk te transformeren met PROJ en voor beide een benadering te gebruiken.

1 like
#8

Bedankt voor de suggestie, het lijkt er op dat ik inderdaad een verouderde PROJ versie geinstalleerd had. Bij de nieuwste versie ging de transformatie goed en kwam de resulterende puntenwolk op de correcte hoogte terecht.

#9

Ik denk wel dat je inputcoördinaten in ETRS89 in plaats van WGS84 zijn en dat je dus de code EPSG:4978 moet corrigeren naar EPSG:4936. Voor de hoogte heeft dat maar een paar centimeter effect, maar voor de xy-coördinaten is het verschil tussen ETRS89 en WGS84 inmiddels 0,8 meter.

PS: Ook als je de EPSG-code voor WGS84 gebruikt krijg je waarschijnlijk de transformatie alsof je input ETRS89 is. Dat komt omdat PROJ WGS84 en ETRS89 aan elkaar gelijkstelt zolang je bovenstaande “ensemble”-codes gebruikt in plaats van EPSG-codes die specificeren welke realisatie van WGS84 en ETRS89 je wil gebruiken. Ik denk dat dat in jouw geval dus juist is, maar mocht je input toch echt WGS84 zijn dan moet je andere codes gebruiken en een epoche specificeren…

1 like