Betrouwbare bron voor proj4 definitie van RD New / EPSG:28992

Ik ben op zoek naar een actuele en betrouwbare bron voor de proj4 definitie string van RD New. Op de github van @antonbakker met voorbeelden zie ik een projectiestring. Waar komt deze vandaan? De proj definitie is in elk geval anders dan bv http://epsg.io/28992. Wellicht dat @lhuisman of @antonbakker mij kunnen helpen?

https://www.nsgi.nl/geodetische-infrastructuur/producten/coordinatentransformatie

De link onder het kopje Parameters en formules RDNAPTRANS leidt naar het aanvragen van een zipje met technische documentatie, waaronder projectie strings

1 like

Michiel,

De PROJ.-string waar je naar verwijst geeft een benadering zonder het correctiegrid en met +towgs84-parameters van een verouderde RDNAPTRANS™. De website epsg,io geeft ook een benadering zonder correctiegrid met andere ook verouderde +towgs84-parameters. De verschillen tussen de +wgs84-parameters zijn in de ordegrootte van een centimeter, de fout door geen correctiegrid te gebruiken is veel groter. De officiële website epsg.org van IOGP vermeldt wel de juiste parameters en het correctiegrid, maar geen PROJ.-strings. De officiële RDNAPTRANS™2018-publicatie op nsgi.nl waar Robin naar verwijst bevat de correctiegridbestanden en (in bijlage 3 van de pdf) de officiële PROJ.6-string.

Voor software die op versie PROJ.4 is blijven steken werkt dat niet. De enige echte oplossing hiervoor is een nieuwere PROJ. gaan gebruiken. Vanuit de NSGI raden we het gebruik van PROJ.4 ten zeerste af! (vanwege enkele bugs met fouten van enkele cm op zee). Desalniettemin kan RDNAPTRANS™2018 wel met PROJ.4 (zie cs2cs-commando’s hieronder), want voorlopig negeert de validatieservice de bugs in PROJ.4 nog.

Groeten, Jochem Lesparre (Rijksdriehoeksmeting)
 

RDNAPTRANS™2018 implementation variant 2 with cs2cs of PROJ.4
From ETRS89 to RD and NAP
cs2cs +proj=lonlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs +to +proj=sterea +lat_0=52.156160556 +lon_0=5.387638889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +wktext +nadgrids=rdtrans2018.gsb +geoidgrids=naptrans2018.gtx +no_defs -f %.4f input.txt > output.txt

From RD and NAP to ETRS89
cs2cs +proj=sterea +lat_0=52.156160556 +lon_0=5.387638889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +wktext +nadgrids=rdtrans2018.gsb +geoidgrids=naptrans2018.gtx +no_defs +to +proj=lonlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs -f %.9f input.txt > output.txt

1 like

@Michiel_F ik gebruik zelf vaak epsg.io, deze website is gemaakt door het maptiler.com (niet gelieerd aan EPSG zelf). In de footer van de site kan je vinden welke versie van de EPSG database uitgeserveerd wordt, op dit moment v9.8. De epsg.io website gebruikt dus de EPSG database als bron van wijsheid. Je kan op de EPSG website ook terugvinden welke datum een release is gedaan, maar dat vereist het aanmaken van een account. EPSG heeft overigens ook een API, en dan heb je geen account nodig, zie https://apps.epsg.org/api/v1/VersionHistory?pageSize=500.

Overigens matcht de gebruikte projectiestring in mijn voorbeeld niet met wat epsg.io nu oplepelt. Ik vermoed dat die waardes uit een eerdere voorbeeld heb over gecopy-paste, hier zou ik aanraden de waarde over te nemen van epsg.io. Ik zal mijn voorbeeld ook even aanpassen.

Wat betreft het advies van @Jochem om proj4 niet te gebruiken, dat lijkt mij in dit geval een moeilijk haalbare kaart. Er zijn voor zover ik weet namelijk nog geen PROJ.6 implementaties in Javascript. Alle Javascript GIS libs leunen voorlalsnog op PROJ.4. Dit geldt ook voor mijn voorbeeld, de Proj4Leaflet library is gebaseerd op proj4js.

Mijn advies om PROJ.4 niet meer te gebruiken mag je eventueel voorlopig negeren. Maar het gebruik van de PROJ.4-string zonder correctiegrid, zoals Anton voorstelt, is echt onverstandig. Misschien heb je de nauwkeurigheid van transformatie met het correctiegrid nu niet nodig, maar vroeg of laat wordt het dan uiteindelijk per ongeluk toch ook gebruikt voor toepassingen waar die hoge nauwkeurigheid wel nodig is.

Deze PROJ.-string zou gewoon moeten werken in PROJ.4:
+proj=sterea +lat_0=52.156160556 +lon_0=5.387638889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +wktext +nadgrids=rdtrans2018.gsb +geoidgrids=naptrans2018.gtx +no_defs
Alle parameterwaarden en de correctiegridbestandsnamen hierin zijn afkomstig uit RDNAPTRANS™2018.pdf van de NSGI. Michel vroeg om een actuele betrouwbare bron, de NSGI is als wettelijke autoriteit voor coördinaten in Nederland de meest actuele en betrouwbare bron (al zeg ik het zelf ; )

@Jochem gaan deze proj strings met correctiegrid ook in de EPSG database landen ? Ik snap dat NSGI de wettelijke autoriteit is die hier over gaat. Maar het overgrote deel van de software in het GIS domein gebruikt de EPSG definities als de standaard. Het opnemen van de proj strings met correctie grid in de EPSG database zou de adaptatie zeker ten goede komen.

@Anton, RDNAPTRANS™2018 is gepubliceerd in augustus 2019 en staat inclusief correctiegrid sinds januari 2020 in EPSG (epsg.org). Het duurt vervolgens alleen ook nog een tijd voor het ook in allerhande software en dergelijke staat (zoals epsg.io, als die überhaupt updates doet).

1 like

Dank voor de antwoorden. Zoals gezegd, Proj6 is nog niet beschikbaar voor javascript applicaties. Ik ben blij dat er nu een “officiele” proj4 definitiestring is voor RD. Ik zal deze gaan gebruiken.
Overigens is een nauwkeurigheid van ~1 meter voldoende.

1 like

Liep hier toevallig ook net tegenaan: BlenderGIS gebruikt epsg.io via een API call. Daaruit komt de “klassieke” EPSG:28992 string met de 6-parameter transformatie zoals die bijna overal gebruikt wordt, bijv voorbeeld Anton. …alleen is m.i. de +towgs84 waarde incorrect. De RDNAPTRANS grid-transformatie is natuurlijk nauwkeuriger. Maar vraag mij af hoe dat in bijv proj.js (JavaScript) in browser moet werken. Je hebt toch dan correctie-grid nodig? Maar goed, gaat hier even om de proj-string uit Amersfoort / RD New - Netherlands - Holland - Dutch - EPSG:28992 voor EPSG:28992. Deze 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 

Als je deze gebruikt zul je merken dat je objecten een paar honderd meter naar het noord-oosten zijn verschoven (vanwege incorrecte +towgs84) op bijv BRT achtergrond kaart in QGIS:
image

De enige goede string, die ik bij @antonbakker zie (en denk ik weer moet terugveranderen naar zijn oorspronkelijke waarde) en ik ook al jaren in mijn viewers als expliciete proj-string vereist is, gebruik, is:

+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel  +towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812 +units=m +no_defs

Dus met

+towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812

Dit issue komt al bijna 20 jaar steeds weer opduiken. Zie bijv [1]. Op epsg.org kan ik overigens geen proj-string voor EPSG:28992 vinden. Gaat de RDNAPTRANS-correctie-grid-gedreven transformatie dan ook EPSG:28992 gebruiken? Lijkt mij niet handig…

[1] http://blog.openstreetmap.nl/index.php/2012/01/21/rd/

Correctiegrid in PROJ.4
Om RDNAPTRANS™2018 te gebruiken moet inderdaad zowel bij PROJ.4 als bij latere PROJ.-versies het bestand rdtrans2018.gsb in de default map staan, of je moet het volledige path in de PROJ.-string vermelden. Of dat werkt (en zo ja hoe) met JavaScript in een browser weet ik niet.

Fout van > 100 m
De fout van > 100 m die jij hier ziet ligt niet aan de parameters (het zijn er overigens 7 i.p.v. 6) die jij gebruikt. De verschillen tussen de +towgs84=parameters die je vermeldt zijn in de ordegrootte van een centimeter.

Een fout van met een ordegrootte van meer dan 100 m in de richting die we hier zien, wordt meestal veroorzaakt door het gebruiken van +towgs84=0,0,0,0,0,0,0. Het lijkt er dus op dat jij (bijna) goede parameters gebruikt, maar dat de achtergrondkaart waarmee je vergelijkt fout ligt.

EPSG:28992
EPSG (.org) publiceert inderdaad geen PROJ.-strings, maar alleen GML en WKT-strings. Bovendien is in EPSG een CRS-code zoals 28992 alleen de definitie van het CRS, en dit omvat NIET de transformatie. Transformaties hebben in EPSG eigen codes. De transformatie conform RDNAPTRANS™2018 heeft transformatiecode 9282, de benaderde transformatie zonder correctiegrid heeft code 9281 (was 4830).

Het vervelende is dat software die de EPSG-registry gebruikt meestal alleen de CRS-codes weergeeft en dan zelf een transformatie kiest zonder te laten zien welke er gebruikt wordt. Daardoor is het inderdaad heel verwarrend dat 28992 zowel voor goed als voor fout getransformeerde coördinaten gebruikt wordt. Extra verwarrend is dat als je uitsluitend met fout getransformeerde RD-coördinaten werkt in GIS-software, je dan niet merkt dat de verkeerde transformatie toegepast wordt, omdat al je data op dezelfde manier verkeerd ligt.

Vanuit EPSG-logica is het slordig van die software om alleen de CRS-code weer te geven zonder de transformatie-code te vermelden. EPSG geeft namelijk nooit een nieuwe CRS-code bij een nieuwe transformatie. De RD-coördinaten van punten in het veld zijn namelijk door de nieuwe transformatie niet veranderd, dus blijft RD de CRS-code 28992 houden. Alleen voor RD-coördinaten van 50 jaar gelden, van voor de verschuiving van de oorsprong van Amersfoort naar Frankrijk, (RD old) is een andere CRS-code.

De manier waarop je hier als gebruiker mee om moet gaan is dus om een CRS-code als 28992 alleen gebruiken als label dat de coördinaten in RD zijn. Of en hoe die coördinaten getransformeerd zijn (met welke precisie), dat moet je als gebruiker zelf bijhouden in de meta-data.

1 like

@Jochem bedankt voor je uitgebreide antwoord. Dit soort issues had ik al lang niet meer gehad en was een deja-vu.

Achtergrond: ben bezig met genereren reliëfkaart (hillshading) in Blender met BlenderGIS plugin met brondata AHN3 DTM TIFF bestanden. Het werken met deze nieuwe tools was al uitdagend genoeg… Op mijn eigen systeem (Mac OS, Homebrew) heb ik Proj 7.2.0 met GDAL 3.2.0 dus heel recent (bijv cs2cs klopt precies). De BlenderGIS (Python) plugin gebruikt niet automatisch Proj voor transformaties maar een API van epsg.io. Die kan ook via de website aangeroepen, voorbeeld [1]. Resultaat is daar m.i. correct. Het resultaat van de hillshading is weer een TIFF die ik met gdal-translate naar GeoTIFF omzet. Die bekijk ik vervolgens in QGIS (laatste versie 3.14) op BRT achtergrondkaart via PDOK Plugin. Daar zag ik dus die verschuiving. Dacht aanvankelijk dat dit lag aan de +towgs84 waarden van epsg.io. Maar je hebt gelijk, die minimale verschillen maken niet uit. Mogelijk veranderen die in de tijd (? moving continents). Maar de deja-vu was uit de tijd dat +towgs84 soms geheel verdwenen was uit de EPSG database van proj (of GDAL, of PostGIS, vaak wel/niet aanwezig meerdere plaatsen).

Lang verhaal kort: heb de BlenderGIS plugin kunnen configureren om mijn lokale pyproj (Python bindings), dus Proj 7.2.0 te gebruiken en vanaf dat moment waren de resultaten correct, pasten precies op BRT. Mijn +towgs84 verhaal eerder boven klopt dus niet (ook t.a.v. @antonbakker).

Dit allemaal even los van RDNAPTRANS2018 en rdtrans2018.gsb. Ik zie mijzelf toch als redelijk ervaren, maar vond het lastig informatie te vinden. Bij Kadaster/NSGI [3] moest ik een formulier invullen. Proj 7.2.0 heeft de grid-correctie file niet automatisch (?), vind deze uiteindelijk in Ubuntu package [2]. De website van NSGI [3] werd ik niet veel wijzer van.

Lijkt mij dat dit voor de beginnende en dus zelfs gevorderde GIS-er toch wel imponerende materie is. Wat we m.i. nodig hebben is iets van een handreiking doc en m.n. voor Proj, wie weet iets voor Geonovum? Alweer meer dan 4 jaar geleden organiseerden wij (Frank Steggink @fsteggink via OSGeo.nl) een mini-seminar “RD en Open Source Software” [4], waarin m.n. Lennard Huisman e.e.a. helder uitlegde ([5] slides) . Misschien weer tijd voor (virtuele) opfris-seminar met laatste stand van zaken! Vanuit OSGeo.nl faciliteren wij dat graag!

[1] Transform coordinates - GPS online converter
[2] https://ubuntu.pkgs.org/20.04/ubuntu-multiverse-amd64/proj-rdnap_2008+2018-4_all.deb.html
[3] https://www.nsgi.nl/geodetische-infrastructuur/producten/coordinatentransformatie
[4] Verslag: mini-seminar RD en Open Source Software – okt 2016 | OSGeo.nl - Wegwijs in Open Geo
[5] http://io.osgeo.nl/sitecontent/events/RDMiniSeminar2016/LennardHuisman.pdf

Formulier
Dat de RDNAPTRANS™2018-download achter een formulier zit, is een erfenis uit het verleden. Het maakt het wel mogelijk om iedereen die de informatie downloadt te vragen of ze per email op de hoogte gehouden willen worden bij updates. Ik had niet het idee dat dit een groot probleem is voor gebruikers, maar ik zou dit inderdaad liever makkelijker maken. Wat denk je van alleen een ja/nee-vinkje voor op de hoogte gehouden worden en dat de download-knop beschikbaar komt zodra je een emailadres ingevuld hebt of nee aangevinkt hebt? Of vind je dat het helemaal zonder vinkjes zetten te downloaden zou moeten zijn? Overigens hebben we geen strikte copyrights meer op de gridbestanden, dus anderen kunnen de gridbestanden zonder formulier herdistrubueren.

RDNAPTRANS™2018-documentatie
De informatie van de download is sowieso vooral bedoeld voor (landmeetkundige) software-ontwikkelaars en niet zo zeer voor individuele (GIS-)gebruikers. De Rijksdriehoeksmeting en NSGI-partners hebben ook niet de personele capaciteit om deze grote groep gebruikers direct te ondersteunen. We zetten daarom in op het ondersteunen van de software-leveranciers om RDNAPTRANS™2018 beschikbaar te maken voor hun gebruikers en dat zij de ondersteuning die daarvoor nodig is kunnen leveren. Maar ook daarvoor is onze capaciteit om proactief ondersteuning bij de implementatie RDNAPTRANS™2018 te bieden beperkt.

Grip op RD en ETRS89
Om wat extra voortgang te boeken bij het breed geïmplementeerd krijgen van RDNAPTRANS™ is in 2016(!) al in het Geonovum-projectrapport Grip op RD en ETRS89 een lijst aanbevelingen gedaan, maar het is tot op heden nog niet gelukt daar financiering voor te krijgen.

Ondertussen hebben we enkele van de aanbevelingen wel al opgepakt, zoals:

  • Aanbieden van een validatieservice
  • Introductie van RDNAPTRANS™2018
  • Opname in de EPSG-registry

Maar de meeste andere plannen laten nog op zich wachten, zoals:

  • Technische ondersteuning:
    • Certificeringsfunctie in de validatieservice (om toestemming voor het gebruik van de merknaam RDNAPTRANS™ aan te vragen).
    • Repository voor RDNAPTRANS™-implementaties in verschillende programmeertalen waar eenieder implementaties kan toevoegen of verbeteren
  • Advies en richtlijnen voor o.a.:
    • Welk CRS waar te gebruiken
    • Hoe om te gaan met met meerdere CRS voor rasterdata en lange lijnen in vectordata
  • Afdwingen gebruik met o.a.:
    • Regelgeving
    • Product-, dienst- of inkoopspecificaties
  • Bewustwording en voorlichting met:
    • Publicaties
    • Bijeenkomsten

Een presentatie van mij bij een opfris-seminar over RD en NAP in open source software past daar dus zeker bij. Onze personele capaciteit zoiets te organiseren is beperkt, dus als OS Geo NL daarin wat kan betekenen, graag! Vaak hebben wij wel budget voor de onkosten van een fysieke bijeenkomst (als dat weer verantwoord is).

PROJ.7
Als het goed is worden de gridbestanden van RDNAPTRANS™2018 nu meegeleverd met PROJ. (als .gsb/.gtx en/of GeoTiff), maar ik heb dat nog niet zelf getest. Ik zal dat binnenkort eens gaan doen, zodat als PROJ.7 standaard met QGIS meegeleverd gaat worden ik weet hoe dat werkt.

Groeten, Jochem (Rijksdriehoeksmeting, NSGI)

2 likes

Gebruikers kunnen in QGIS wel een verwijzing opgeven naar het gridbestand. Op de website van de NSGI is dat niet zomaar te vinden als losse download. Met wat zoeken vind je het gridbestand op wat obscuur ogende websites, dus heel gebruiksvriendelijk is het niet.

Als het straks met QGIS default wordt meegeleverd is het nog mooier, maar anders zou een rechtstreeks linkje op NSGI naar het gridbestand best fijn zijn.

Overigens geweldig wat je hier aan informatie geeft, zeer waardevol! :+1:

1 like

Het zelf moeten configureren van RDNAPTRANS™2018 als custom CRS in QGIS is sowieso niet erg gebruiksvriendelijk. Ik heb wel een conceptversie van een handleiding daarvoor gemaakt, maar die is nog niet gepubliceerd. Voor de meeste gebruikers wordt het waarschijnlijk toch wachten tot het standaard in QGIS zit.

In QGIS 3.20 is het correctiegrid nu aanwezig. Uit de betreffende readme onder PROJ:

# nl_nsgi_README.txt

The files in this section result from the conversion of datasets originating
from [NSGI](https://www.nsgi.nl/)

## Included grids

### Netherlands: RD and NAP -> ETRS89

*Source*: [NSGI](https://www.nsgi.nl/geodetische-infrastructuur/producten/coordinatentransformatie)  
*Format*: GeoTIFF converted from NTv2 & GTX  
*License*: [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/)  

Het GSB bestand is dus niet meer nodig, er zijn nu GeoTiffs beschikbaar met deze data. In QGIS:

In eerdere versies was hier een button aanwezig waar je een pad op kon geven naar het benodigde GSB bestand. Dat hoeft dus niet meer.

2 likes

De GeoTIFFs in QGIS zijn door PROJ gemaakt op basis van het NTv2-grid (.gsb) en VDatum-grid (.gtx) van de NSGI. Wij van de NSGI hebben de transformatie met de GeoTIFFs gecontroleerd. De meta-data is nog niet helemaal goed (hier en daar staan er nog verkeerde omschrijvingen in), maar de coördinaten die je bij transformatie tussen RDNAP en ETRS89 krijgt zijn wel 100% correct volgens onze validatieservice (1 mm).

Let op: Door de manier waarop de EPSG Registry georganiseerd is, gaat de transformatie tussen RDNAP en ITRF2014, WGS84(G1762) e.d. helaas nog niet automatisch op 1 mm goed in QGIS.

Let op! Om de hoogte helemaal goed volgens RDNAPTRANS™2018 te berekenen moet in PROJ4 de string die ik hierboven noem gebruikt worden. Vanaf PROJ6 moet voor een string in PROJ4-stijl een ander gridbestand voor de quasi-geoïde gebruikt worden:

+proj=sterea +lat_0=52.156160556 +lon_0=5.387638889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +nadgrids=rdtrans2018.gsb +units=m +geoidgrids=nlgeo2018.gtx +vunits=m +no_defs +type=crs