Gdalwarp pipeline voor RDNAPTRANS2018

Standaarden RDNAPTRANS
#1

Beste mensen,

Ik probeer om met gdalwarp en RDNAPTRANS2018 geotifs te projecteren. Geotif in ETRS89: xy en met rasterwaarde z (hoogte).

Ik heb de documentatie van www.nsgi.nl/geodetische-infrastructuur/producten/coordinatentransformatie gelezen en de gridbestanden gedownload (rdcorr2018.gsb) en kan het cct voorbeeld van de 3D transformatie in Appendix 3 reproduceren. Maar ik ben nu op zoek naar de manier om dit met gdalwarp te doen voor mijn DEM.

Ik heb GDAL 3.0.4 en Proj 7.0.0 geinstalleerd op linux.

Wie kan me helpen?

Mvg.
Eelke

#2

Ik weet niet zeker of ik je vraag goed begrijp. Zoals jij het formuleert, interpreteer ik het als dat je een Geotif in ETRS89 hebt en dat je die om wil zetten naar een Geotif in RD&NAP.

  • Is het raster in ongeprojecteerd ETRS89 of in een of andere projectie van ETRS89?
  • Is de rasterwaarde in je Geotif ellipsoĂŻdische ETRS89-hoogte, NAP, of EGM?
  • Is je inputdata ĂĽberhaupt alleen in ETRS89 beschikbaar? Voor wereldwijde data is dat onlogisch en Nederlandse data is meestal (ook) in RD&NAP beschikbaar.
  • Wil je als output weer een Geotif of is ascii of ander formaat ook goed?
#3

  • ongeprojecteerd

  • inderdaad ellipsoĂŻdische ETRS89-hoogte (~40m boven NAP)

  • Deze DEM is output van photogrammetry; hier is epsg:4258 gekozen als output CRS

  • graag weer als geotif (maar met een extra gdal stap kan volgens mij ascii → geotif vrij makkelijk)

Hieronder voor de volledigheid de gdalsrsinfo output

PROJ.4 : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

OGC WKT2:2018 :
BOUNDCRS[
SOURCECRS[
GEOGCRS[“ETRS89”,
DATUM[“European Terrestrial Reference System 1989”,
ELLIPSOID[“GRS 1980”,6378137,298.257222101004,
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”,4258]]],
TARGETCRS[
GEOGCRS[“WGS 84”,
DATUM[“World Geodetic System 1984”,
ELLIPSOID[“WGS 84”,6378137,298.257223563,
LENGTHUNIT[“metre”,1]]],
PRIMEM[“Greenwich”,0,
ANGLEUNIT[“degree”,0.0174532925199433]],
CS[ellipsoidal,2],
AXIS[“latitude”,north,
ORDER[1],
ANGLEUNIT[“degree”,0.0174532925199433]],
AXIS[“longitude”,east,
ORDER[2],
ANGLEUNIT[“degree”,0.0174532925199433]],
ID[“EPSG”,4326]]],
ABRIDGEDTRANSFORMATION[“Transformation from ETRS89 to WGS84”,
METHOD[“Position Vector transformation (geog2D domain)”,
ID[“EPSG”,9606]],
PARAMETER[“X-axis translation”,0,
ID[“EPSG”,8605]],
PARAMETER[“Y-axis translation”,0,
ID[“EPSG”,8606]],
PARAMETER[“Z-axis translation”,0,
ID[“EPSG”,8607]],
PARAMETER[“X-axis rotation”,0,
ID[“EPSG”,8608]],
PARAMETER[“Y-axis rotation”,0,
ID[“EPSG”,8609]],
PARAMETER[“Z-axis rotation”,0,
ID[“EPSG”,8610]],
PARAMETER[“Scale difference”,1,
ID[“EPSG”,8611]]]]

#4

In dat geval zou je dus met GDAL de Geotif om kunnen zetten naar ascii lat,lon,h; dat met PROJ volgens RDNAPTRANS transformeren naar x,y,NAP; en met GDAL die ascii weer als Geotif schrijven? Maar misschien weet iemand met meer kennis van GDAL of kan dat handiger kan.

#5

Let op: Na het omzetten van 3D ETRS89 naar RD&NAP is het geen regelmatig grid meer. Je zult het dus ook moeten resamplen en daarbij voor een resamplingmethode moeten kiezen.

#6

Dank voor je suggestie Jochem. gdalwarp doet ook resampling (waarbij uit verschillende methodes gekozen kan worden). Ik wacht even andere suggesties af, want ik vermoed dat het in een enkele gdal stap moet kunnen.

1 like
#7

Omdat er nog niemand met ervaring met gdalwarp gereageerd heeft, heb ik zelf even naar de documentatie gekeken. Het makkelijke commando is (met default nearest-neigbour resampling!):

gdalwarp -s_srs EPSG:4258 -t_srs EPSG:7415 input.tif output.tif

Maar omdat RDNAPTRANS2018 pas na PROJ.7 in de EPSG registry is gekomen, geeft dat waarschijnlijk fouten tot 25 cm horizontaal en misschien enkele meters verticaal (als er verticaal al getransformeerd wordt).

Je zult voor RDNAPTRANS2018 dus zelf een PROJ- of WKT-string op moeten geven. Ik denk dat het dan zoiets moet worden:

gdalwarp -ct '+proj=pipeline +step=...' input.tif output.tif

Waar je tussen de aanhalingstekens de PROJ-string moet zetten uit de RDNAPTRANS-download (bijlage 3): Aanvragen download RDNAPTRANS2018
Als je data niet al te ver buiten Nederland en EEZ valt, dan kun je de kortere PROJ-string voor variant 2 gebruiken (paragraaf A3.2).

1 like
#8

Bedankt.

Ik heb met Variant1 dit geprobeerd:

gdalwarp -ct "+proj=pipeline \
+step +proj=axisswap +order=2,1,3,4 \
+step +proj=vgridshift +grids=nlgeo2018.gtx \
+step +proj=axisswap +order=1,2,4,-3 \
+step +proj=vgridshift +inv +grids=egm96_15.gtx \
+step +proj=cart +ellps=GRS80 \
+step +proj=helmert +x=-565.7346 +y=-50.4058 +z=-465.2895 +rx=-0.395023 +ry=0.330776 +rz=-1.876073 +s=-4.07242 +convention=coordinate_frame +exact \
+step +proj=cart +inv +ellps=bessel \
+step +proj=hgridshift +inv +grids=rdcorr2018.gsb,null \
+step +proj=sterea +lat_0=52.156160556 +lon_0=5.387638889 \
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel \
+step +proj=unitconvert +z_in=m +z_out=1e32 \
+step +proj=axisswap +order=1,2,-4,-3" test_4258.tif test_rdnap.tif

en met variant2:

gdalwarp -ct "+proj=pipeline \
+step +proj=axisswap +order=2,1,3,4 \
+step +proj=hgridshift +inv +grids=rdtrans2018.gsb \
+step +proj=vgridshift +grids=naptrans2018.gtx \
+step +proj=sterea +lat_0=52.156160556 +lon_0=5.387638889 \
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel" test_4258.tif test_rdnap.tif

In beide gevallen krijg ik deze fout:
ERROR 1: PROJ: proj_create: Error -38: failed to load datum shift file
ERROR 6: Cannot instantiate pipeline +proj=pipeline +step +proj=axisswap +order=2,1,3,4 … etc

Ik voer het commando uit in de folder waar de benodigde gridshift bestanden zijn opgeslagen. Bij Variant1 heb ik ook het egm96_15.gtx bestand in de folder staan.

#9

Misschien helpt het om het hele pad naar de gridbestanden op te geven? Bijvoorbeeld:

C:\Programs\QGIS3.14\share\proj\rdcorr2018.gsb

of om te specificeren dat het in de huidige folder staat:

.\rdcorr2018.gsb
#10

gdalinfo rdtrans2018.gsb en gdalinfo naptrans2018.gtx geeft onderstaande metadata en ik denk daarom dat deze bestanden wel gelezen worden. Ik zoek verder.

Driver: NTv2/NTv2 Datum Grid Shift
Files: rdtrans2018.gsb
Size is 61, 61
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            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",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (1.950000000000000,56.049999999999997)
Pixel Size = (0.100000000000000,-0.100000000000000)
Metadata:
  CREATED=22-11-18
  GS_TYPE=SECONDS
  MAJOR_F=6377397.155
  MAJOR_T=6378137
  MINOR_F=6356078.963
  MINOR_T=6356752.314
  PARENT=NONE
  SUB_NAME=NL_EEZ
  SYSTEM_F=BSSL1841
  SYSTEM_T=ETRS89
  UPDATED=
  VERSION=NTv2.0
Subdatasets:
  SUBDATASET_0_NAME=NTv2:0:rdtrans2018.gsb
  SUBDATASET_0_DESC=NL_EEZ
  SUBDATASET_1_NAME=NTv2:1:rdtrans2018.gsb
  SUBDATASET_1_DESC=NL
Corner Coordinates:
Upper Left  (   1.9500000,  56.0500000) (  1d57' 0.00"E, 56d 3' 0.00"N)
Lower Left  (   1.9500000,  49.9500000) (  1d57' 0.00"E, 49d57' 0.00"N)
Upper Right (   8.0500000,  56.0500000) (  8d 3' 0.00"E, 56d 3' 0.00"N)
Lower Right (   8.0500000,  49.9500000) (  8d 3' 0.00"E, 49d57' 0.00"N)
Center      (   5.0000000,  53.0000000) (  5d 0' 0.00"E, 53d 0' 0.00"N)
Band 1 Block=61x1 Type=Float32, ColorInterp=Undefined
  Description = Latitude Offset (arc seconds)
Band 2 Block=61x1 Type=Float32, ColorInterp=Undefined
  Description = Longitude Offset (arc seconds)
Band 3 Block=61x1 Type=Float32, ColorInterp=Undefined
  Description = Latitude Error
Band 4 Block=61x1 Type=Float32, ColorInterp=Undefined
  Description = Longitude Error
#11

Doorgaan met de discussie van Gdalwarp pipeline voor RDNAPTRANS2018:

Ik vroeg me af of het je gelukt is @Eelke. Het werkte bij mij om het hele pad te gebruiken voor de rasters. Waar ik mee zit is de omgekeerde transformatie, namelijk van RD/NAP naar ETRS84 (3D). Dan kan ik dezelfde methode gebruiken met een -I vlag wanneer ik proj gebruik (voor enkele coördinaten). In GDAL bestaat de -I vlag niet, hebben jullie daar een idee over? Over hoe ik de transformatie kan omkeren voor een raster?

#13

Deze PROJ-string specificeert alleen de projectie en de geoĂŻde, maar niet het RD-correctiegrid en de datumtransformatie. Ik verwacht dat het daardoor een horizontale fout van 100 meter geeft!

#14

Oei, excuus daarvoor. Ik moet zoeken naar de juiste string en testen. Ik verwijder mijn foutieve suggestie.

1 like
#15

Bedankt voor de suggestie, het is helaas precies het omkeren van de transformatie wat niet voor GDAL gedocumenteerd lijkt te zijn.

#16

Werkt het omkeren van source en target in het makkelijke commando (met default nearest-neigbour resampling) niet gewoon? Dus:

gdalwarp -s_srs EPSG:7415 -t_srs EPSG:4258 input.tif output.tif

Voorwaarde daarvoor is dat je een GDAL-versie met een recente PROJ-versie (en dus EPSG-database) gebruikt, anders krijg je waarschijnlijk fouten tot 25 cm horizontaal en misschien enkele meters verticaal (als er verticaal al getransformeerd wordt). Anders moet het met een inverse PROJ-string.

gdalwarp -ct '+proj=pipeline +step=...' input.tif output.tif

Waar ik de inverse PROJ-string tussen de aanhalingstekens dan speciaal voor deze toepassing zou moeten opstellen. Dus een -I optie zou handiger zijn…

#17

Hoi Jochem,

Mijn doel is van RD/NAP (EPSG:7415) naar EPSG:4326 om te zetten. De achterliggende gedachte is een recente DEM te genererenvan AHN data met pixelgrootte van 10-20 meter.

Als ik het dus goed begrijp is

gdalwarp -s_srs EPSG:7415 -t_srs EPSG:4258 input.tif output.tif

hetzelfde als de volgende (RDNAPTRANS) transformatie met een -I vlag

+proj=pipeline
+step +proj=unitconvert +xy_in=deg +xy_out=rad
+step +proj=axisswap +order=2,1
+step +proj=vgridshift +grids=nlgeo2018.tif
+step +proj=hgridshift +inv +grids=rdtrans2018.tif
+step +proj=sterea +lat_0=52.156160556 +lon_0=5.387638889
+k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel

Indien dat zo is, was mijn veronderstelling dat de gdalwarp met -s_srs en -t_srs niet het goede antwoord zou geven onjuist (voor de latere GDAL versies). Kan ik dan stellen dat (andere -t_srs)

gdalwarp -s_srs EPSG:7415 -t_srs EPSG:4326 input.tif output.tif

net zo goed werkt? Met in het achterhoofd dat dat CRS een (veel) beperktere nauwkeurigheid heeft (ca. 2 meter).

#18

Net niet helemaal want EPSG:4258 is 2D en EPSG:7415 is 3D. Voor 3D RDNAPTRANS van RDNAP naar ETRS89 moet je EPSG:4937 gebruiken:

gdalwarp -s_srs EPSG:7415 -t_srs EPSG:4937 input.tif output.tif

En voor 2D RDNAPTRANS moet je EPSG:28992 gebruiken:

gdalwarp -s_srs EPSG:28992 -t_srs EPSG:4258 input.tif output.tif

Dan zou het goed moeten gaan. Mits je een recente GDAL-versie gebruikt met een PROJ-versie met een recente EPSG-database. De beste manier om te checken of de transformatie zeker weten goed gaat is testen met de NSGI-coördinatentransformatievalidatieservice. Edit: Dat zou je dan moeten testen met een ander GDAL-commando of met PROJ, want de validatieservice werkt met coördinaten van punten en niet met GeoTIFFs.

Transformeren naar EPSG:4326 (2D WGS 84) zou ik niet doen. Het advies van NSGI en Geonovum in de Handreiking CRS is namelijk om WGS 84 (voor de meeste toepassingen) gelijk te stellen aan ETRS89 met een nultransformatie. Als je GDAL vraagt naar WGS 84 te transformeren, loop je het risico dat GDAL toch een of andere (foutieve) transformatie uit gaat voeren.

#19

Bedankt Jochem, dat is compleet helder. Overigens is de Handreiking CRS ontzettend goed om erbij te hebben, bedankt!

Als laatst, heb je een idee om welke (GDAL, EPSG-database) versies het zou gaan?

#20

De laatste wijziging n.a.v. RDNAPTRANS™2018 in de EPSG-database is EPSG v10.013, 2021-02-05. Vanaf PROJ 9.0.1 zijn deze wijzigingen ook opgenomen in PROJ. Maar mogelijk gaat het met eerdere versies ook al goed.

1 like