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_VERTICALin 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 thespecification of a CRS in a Geodetic TIFF
grid file:
* Are the simpler EPSG codes preferred over WKT strings?
Yes
* Is it OK to use a WKT string for the source and a EPSG code for
the target,
Yes
* 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.
*
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
<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 <[email protected]>
*Sent:* maandag 24 januari 2022 21:17
*To:* Lesparre, Jochem <[email protected]>; [email protected]
*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
2. 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.
3. 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>
4. Is it necessary for PROJ to use new files names to distinguish
them from the old version with incorrect metadata?
No
5.
Regards, Jochem
[1] https://github.com/OSGeo/PROJ-data/tree/master/nl_nsgi
<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
<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
<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
<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
[email protected] <mailto:[email protected]>
https://lists.osgeo.org/mailman/listinfo/proj
<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].
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
PROJ mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/proj
