[PROJ] Correcting metadata in GeoTIFF files

Lesparre, Jochem Jochem.Lesparre at kadaster.nl
Wed Jan 26 11:54:17 PST 2022 

Hi Even,

Thanks again for your quick reply, but I think you misunderstood the problem I tried to describe in my last email. I'll try to explain our problem a bit more comprehensive:

We distribute two grid files for the quasi-geoid of the Netherlands: NLGEO2018 and NAPTRANS2018. Both expect ETRS89 (EPSG:4937) ellipsoidal height as input and give the same NAP national levelling system height (EPSG:5709) as output. However, the latlon coordinates of NLGEO2018 are in ETRS89 and the latlon coordinates of NAPTRANS2018 are in the national triangulation system "Amersfoort" (EPSG:4289). We prefer the use of the NLGEO2018. We supply NAPTRANS2018 only because this is the type of VDatum grid file that was needed for 3D transformation with PROJ4 (see [1] in my previous email).

Creating a Geodetic TIFF Grid file for NLGEO2018 is straight-forward [A]. However, I can't create a Geodetic TIFF Grid file for NAPTRANS2018. If I use GEOGRAPHIC_TO_VERTICAL the --interpolation-crs is ignored [B]. If I use VERTICAL_TO_VERTICAL, I get AssertionError due to the use of a 3D instead of a vertical --source-crs [C]. This confuses me, since our use case would probably not be the only VDatum grid with this problem.


[A] Creating a Geodetic TIFF Grid file for NLGEO2018 is straight-forward:
python vertoffset_grid_to_gtiff.py --type GEOGRAPHIC_TO_VERTICAL --area-of-use "European Netherlands including EEZ of the North Sea" --source-crs EPSG:4937 --target-crs EPSG:5709 --copyright "Nederlandse Samenwerking Geodetische Infrastructuur (NSGI). Creative Commons Attribution 4.0 International licence https://creativecommons/licenses/by/4.0/" nlgeo2018.gtx nlgeo2018.tif

Result: OK


[B] Creating a Geodetic TIFF Grid file for NAPTRANS2018 as I expected it to be done:
python vertoffset_grid_to_gtiff.py --type GEOGRAPHIC_TO_VERTICAL --area-of-use "European Netherlands including EEZ of the North Sea" --source-crs EPSG:4937 --target-crs EPSG:5709 --interpolation-crs EPSG:4289 --copyright "Nederlandse Samenwerking Geodetische Infrastructuur (NSGI). Creative Commons Attribution 4.0 International licence https://creativecommons/licenses/by/4.0/" naptrans2018.gtx naptrans2018.tif

Result: interpolation-crs ignored


[C] NAPTRANS2018 using VERTICAL_TO_VERTICAL instead
python vertoffset_grid_to_gtiff.py --type VERTICAL_TO_VERTICAL --area-of-use "European Netherlands including EEZ of the North Sea" --source-crs EPSG:4937 --target-crs EPSG:5709 --interpolation-crs EPSG:4289 --copyright "Nederlandse Samenwerking Geodetische Infrastructuur (NSGI). Creative Commons Attribution 4.0 International licence https://creativecommons/licenses/by/4.0/" naptrans2018.gtx naptrans2018.tif

Result: AssertionError due to 3D instead of vertical source-crs

Regards, Jochem



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



Le 25/01/2022 à 15:43, Lesparre, Jochem a écrit :
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.
You should use --type GEOGRAPHIC_TO_VERTICAL if it is a geoid file (the name is perhaps a bit misleading as the values are to go from vertical CRS to geographic CRS but this is the EPSG naming convention for the transformation)

The latter doesn't have a vertical EPSG code. What is the WKT string of ETRS89 ellipsoidal height?
You should use EPSG:4937 for ETRS89 geographic 3D


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:

1.       Are the simpler EPSG codes preferred over WKT strings?
Yes


2.       Is it OK to use a WKT string for the source and a EPSG code for the target,
Yes



3.       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])?
You likely need to remove the ID["EPSG",4289] node. And it is a bit weird to still use Amesfoort as the datum, since, as far as I understand, if 2 geographic 2D CRS that have same axis, units, datum and prime meridian, they should be the same... That said that will have little consequences since PROJ ignores that anyway.


4.

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><mailto:even.rouault at spatialys.com>
Sent: maandag 24 januari 2022 21:17
To: Lesparre, Jochem <Jochem.Lesparre at kadaster.nl><mailto:Jochem.Lesparre at kadaster.nl>; PROJ at lists.osgeo.org<mailto: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<https://secure-web.cisco.com/1Rhnb2SuM-ZisQ5-ehSorBEaEAymp3uWVHsbgnlOs1Pg5MgVrsPjC7O8qTjD4zT80_6Xn1RawtVUWqlKkCdVIMiapAnSV6eAnPZUsK2JJhMZsWvP4SEVK9tdjcphLVMKBb8TXeu7NFSS6QXoxNn9xjRjduRUV1SOtN0x7PgayMSReCgEDVClEN8TQJxMvxRfFhCldHWr6FP2eavQotePbR4Ywg-kQn1stVzboPN0PeuHNKOVfJWhmqp4AIe58aZyLwgTXIENtf7myj6SP_fzPT0RKDzdxx5Bn7Qnh0FcCeUhlH8n2FTZNj6O-eHXprtRCq8RPUiUODDiRzvSnHATXvrDLhbXeryc5xMaVfHnlhsA/https%3A%2F%2Fwww.kadaster.com%2Fgeneral-terms-and-conditions>].

--

http://www.spatialys.com<http://secure-web.cisco.com/1LU8njMfYmMFaFPapILcBM73tvtUOgBd-EEc7zqR2Vk9u3neho5DJev426ic0f46uR_eq5VF0Bqp-RFLefFqoLhaF2inDhGAZgnd5N6ho4piq66Ib4z9Ofr5VTjgAv2QzcdUtWHHHgmkCtRTNBrW46MafOphoV-gY7zxxMHnmefFFNYAHWH2dTmXr73fCLBLmrXOKOzWA8RNSBPR33Lyktcr0Fh3qpSVbr_HPcjjd4QDQWHP__uoaASTHBzE3OXCuYDzm5BeqrVuM3T9qe_P1bBy78Apc6PItZR4oQfxICykNd2aay1UZnhTBVWzaWp6AD3K8lx-f7LYv50X_XY8Oz8hPfe-0Rt2IObmC5axZy3Y/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/20220126/e58699a4/attachment-0001.html>

More information about the PROJ mailing list