[PROJ] Correcting metadata in GeoTIFF files

Tue Jan 25 06:43:32 PST 2022

Thanks.

I encountered a new problem:
We use a horizontal grid as well as a vertical grid for the transformation from ETRS89 (EPSG:4937) to our national CRS: RD (EPSG:28992) with NAP height (EPSG:5709). With PROJ4-style strings [1], an interpolation CRS (EPSG:4289) is needed for the Geodetic TIFF grid of the geoid, next to the source (EPSG:4937) and target (EPSG:5709) CRS. To specify this, I have to use --type VERTICAL_TO_VERTICAL in the Python script for the conversion from the VDatum grid. This results in the description "vertical offset" instead of "geoid undulation". This is a bit strange, as the grid is in fact a geoid. Problematic is that (according to the documentation [2]) the source and target CRS both must be a vertical CRS, while one is our national levelling height system NAP, but the other is ellipsoidal ETRS89 height. The latter doesn't have a vertical EPSG code. What is the WKT string of ETRS89 ellipsoidal height? This seems way too complicated for a Geodetic TIFF grid file of a VDatum geoid grid file of a type that must have been very common under PROJ4.

Next to this, I have some new questions on the use of WKT strings instead of EPSG codes for the specification of a CRS in a Geodetic TIFF grid file:

  *   Are the simpler EPSG codes preferred over WKT strings? Is it OK to use a WKT string for the source and a EPSG code for the target, or would it be better to use a WKT string for both source and target in such case?
  *   Can I use the WKT string given by EPSG.org for a CRS with an EPSG code and edit the CRS name in this WKT string to use an alias instead of the official EPSG name (see example [3])?

Regards, Jochem

[1] PROJ4-style transformation from ETRS89 (EPSG:4937) to national horizontal and vertical CRS of the Netherlands:
cs2cs +init=epsg:4937 +to +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=naptrans2018.gtx +vunits=m +no_defs

[2] https://github.com/OSGeo/PROJ-data/blob/master/grid_tools/README.md

[3] WKT string of EPSG:4289 where I replaced the official EPSG name "Amersfoort" to the official Dutch name "RD Bessel":
GEOGCRS["RD Bessel",
   DATUM["Amersfoort",
      ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
         LENGTHUNIT["metre",1,
            ID["EPSG",9001]
         ],
         ID["EPSG",7004]
      ],
      ID["EPSG",6289]],
   CS[ellipsoidal,2,
      ID["EPSG",6422]
   ],
   AXIS["latitude (Lat)",north],
   AXIS["longitude (Lon)",east],
   ANGLEUNIT["degree",0.0174532925199433,
      ID["EPSG",9102]
   ],
   ID["EPSG",4289]
]

From: Even Rouault <even.rouault at spatialys.com>
Sent: maandag 24 januari 2022 21:17
To: Lesparre, Jochem <Jochem.Lesparre at kadaster.nl>; PROJ at lists.osgeo.org
Subject: Re: [PROJ] Correcting metadata in GeoTIFF files

Jochem,
Le 24/01/2022 à 20:09, Lesparre, Jochem via PROJ a écrit :

Dear list,

I'm from Netherlands Partnership Geodetic Infrastructure (NSGI). We published NTv2 (.gsb) and VDatum (.gtx) grid files for the national CRS of the Netherlands (RD coordinates with NAP height) in 2019. PROJ converted these to GeoTIFF (.tif) in 2020 [1]. Unfortunately, some wrong information was stored in the metadata of the GeoTIFF files in the conversion. We are now creating new GeoTIFF files with corrected metadata, using the Python scrips ntv2_to_gtiff.py [2] and vertoffset_grid_to_gtiff.py [3].

For most changes in the metadata this is straightforward, but there is one more difficult issue:

We have two variants of the horizontal transformation:

1.       Conventional 1-step transformation (variant 2): geographic coordinates of national CRS ---[nl_nsgi_rdtrans2018.tif]---> ETRS89

2.       Better 2-step transformation (variant 1): geographic coordinates of national CRS ---[nl_nsgi_rdcorr2018.tif]---> corrected geographic coordinates of national CRS ---[7_parameter_transformation]---> ETRS89

Since there is no separate EPSG code for the corrected geographic coordinates of the national CRS, we want to use the same EPSG code for both corrected and uncorrected coordinates. Or will it give problems when the source and target CRS of a GeoTIFF file are the same?

Instead of the target_crs_code metadata item, you could include a target_crw_wkt with a WKT2 string. That can be done through the ntv2_to_gtiff.py script as it can accept a WKT2 CRS string as the value for --target-crs.

Next to this, I have some other questions:

  1.  Shouldn't a GeoTIFF grid file for a vertical transformation have an accuracy band like a grid for a horizontal transformation? Is it possible with the Python script to create an accuracy band in the GeoTIFF from a VDatum grid file?

You could possibly add with GDAL standard tools (gdal_translate, etc) a band with the accuracy and a description of "geoid_undulation_accuracy" (you may need to go through a VRT to manually add the Units to the band). The validate_vertical_offset_geographic_to_vertical() method of the check_gtiff_grid.py script would likely have to be updated so that it doesn't emit an information message about the new band not being recognized. And the spec at https://proj.org/specifications/geodetictiffgrids.html<https://secure-web.cisco.com/1IPsQrcS7XiA25_SIbwFzb-QZAKsmLdiZCviwwLfntX_LP5RL4MhT86HqkIDLHoWnxjGPsOj9qatKSQRQ6s-3OgTu2NWmUxO08eXFQDCyBNoiM6PgcO9f94P0uPY5Tod4uKHZ2DYmeBgeH9t6R_uf9AxSLI6uZ7h4ReNQ8fupVoj7z9uP0to6Nlr1Kti7lJycLJvZgYL1hGJywSBwFWky-TKX6tQWP90b26Q0Iub853RFhi4Meej8qTYaHcMfPtczIIxD62IAkCtQyIv3XMdTneUYyFMRszve23cNIudpKQzeuZzt3yS_yl6bQIPnPMO6/https%3A%2F%2Fproj.org%2Fspecifications%2Fgeodetictiffgrids.html> as well

  1.  Can I add a recommended_interpolation_method with the Python scripts?
Not currently.  Either enhance them or add the metadata item manually through GDAL tools (gdal_translate -mo recommended_interpolation_method=foo in.tif out.tif -co COMPRESS=DEFLATE -co INTERLEAVE=BAND -co PREDICTOR=3). Note that PROJ itself will ignore that item and always apply bilinear interpolation.

  1.  How should I supply the corrected GeoTIFF files to PROJ?
Pull request against https://github.com/OSGeo/PROJ-data<https://github.com/OSGeo/PROJ-data/tree/master/nl_nsgi>

  1.  Is it necessary for PROJ to use new files names to distinguish them from the old version with incorrect metadata?
No

  1.

Regards, Jochem

[1] https://github.com/OSGeo/PROJ-data/tree/master/nl_nsgi

[2] https://github.com/OSGeo/PROJ-data/blob/master/grid_tools/ntv2_to_gtiff.py

[3] https://github.com/OSGeo/PROJ-data/blob/master/grid_tools/vertoffset_grid_to_gtiff.py

J. Lesparre

Netherlands Partnership Geodetic Infrastructure (NSGI.nl)

Disclaimer:
De inhoud van deze e-mail is vertrouwelijk en uitsluitend bestemd voor de geadresseerde(n).
Gebruik, openbaarmaking, vermenigvuldiging, verspreiding en/of verstrekking van deze informatie aan derden is niet toegestaan.
Op al onze producten en diensten zijn onze algemene leveringsvoorwaarden van toepassing
[https://www.kadaster.nl/algemene-leveringsvoorwaarden].

Disclaimer:
This email and any files transmitted with it are confidential and intended solely for the use of the individual or entity to whom they are addressed.
If you are not the intended recipient, you are notified that disclosing, copying, distributing or taking any action in reliance on the contents of this information is strictly prohibited.
Our general terms and conditions of delivery apply to all our products and services
[https://www.kadaster.com/general-terms-and-conditions<https://secure-web.cisco.com/1fU8sg4FrMFO4TyIW9ZofFSwdqxB0Kag8iZvKOseM_86a_Jx5Z2YV1Hrtt_e8mtlHLMQ0H8zu1i3lanHmsFjggZd1QXIC7i_zohRumwvppYst-xwoiyIzVpRIaUe7MImO2MSCS5aoznYYHb6E-LdzL2tarUM3ByI1fBzPf6EX-yQBpfd_vLLLUaYH2he7-jy3HzpLC2iu2ee6zOAfHGtKhEMeppr62m_f186VKNfXWlLwAxhGXGppzPfI_g4RpQ9OMJiT8koLCKdzzpwTMdFYDY2DaE08KWCvO9WNaB0ciQrOfxgdxfjiQ802OXi4mhKe/https%3A%2F%2Fwww.kadaster.com%2Fgeneral-terms-and-conditions>].

_______________________________________________

PROJ mailing list

PROJ at lists.osgeo.org<mailto:PROJ at lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/proj

--

http://www.spatialys.com<http://secure-web.cisco.com/1fXZS9zkl1X2JBd5DEK4YAq1gPWoDfSCpDSaFlOf5HpAUYrV7SIEw12ok2xJ1dLlvazjXyElbKEjEJLnTBo5NsMUQgMzUh1-_47mgVQpSA2Wih89heEKMsI5iXG5QbxXLAERdl1W-2up80jCvknOD8vY83_Wh-5g43WZIwWQNwaJ38eTscgBCOajh8Bgq8Xmw06STQekuuToXJWmYwQ6woQiNYlE3tHnjcvuhJOGF6aedm4xTRlC4hjYeQ9b8T3HPmwawXXFsXFnlPiuOSR06jIeZ0FCvBXtbexj0ZloQCmzM8W3sBGc11dNGXri3ia7w/http%3A%2F%2Fwww.spatialys.com>

My software is free, but my time generally not.

Disclaimer:
De inhoud van deze e-mail is vertrouwelijk en uitsluitend bestemd voor de geadresseerde(n).
Gebruik, openbaarmaking, vermenigvuldiging, verspreiding en/of verstrekking van deze informatie aan derden is niet toegestaan.
Op al onze producten en diensten zijn onze algemene leveringsvoorwaarden van toepassing
[https://www.kadaster.nl/algemene-leveringsvoorwaarden].

Disclaimer:
This email and any files transmitted with it are confidential and intended solely for the use of the individual or entity to whom they are addressed.
If you are not the intended recipient, you are notified that disclosing, copying, distributing or taking any action in reliance on the contents of this information is strictly prohibited.
Our general terms and conditions of delivery apply to all our products and services
[https://www.kadaster.com/general-terms-and-conditions].
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20220125/596f691d/attachment-0001.html>