[PROJ] Correcting metadata in GeoTIFF files

Lesparre, Jochem Jochem.Lesparre at kadaster.nl
Mon Jan 31 03:45:20 PST 2022


Hi Even,

> Looking at the 'epsg' file that was delivered in PROJ < 6, there's no single entry having both +nadgrids and +geoidgrids in it.

That just means that countries did not take the effort to register their national projected 2D CRS and national levelling height system as a compound CRS in EPSG. Or that the precise transformation grids are not registered in EPSG. However, I expect that most countries in the world will have at least one such compound CRS in use which would need a horizontal and vertical grid.

> I'm not completely sure if the fact that the order of those operations was the one you expected was a design decision or just the result of a 50% 50% dice throw.

We noticed that PROJ4 used to do the horizontal transformation first, so we created our naptrans2008.gtx based on “Amersfoort” instead of ETRS89 horizontally. We didn’t like this, but this was the only way of getting the hight transformation exactly right.

> Hum but now I'm lost, if doing 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"

- with recent PROJ, it  does ... +step +inv +proj=vgridshift +grids=naptrans2018.gtx +multiplier=1  +step +inv +proj=hgridshift +grids=nl_nsgi_rdtrans2018.tif , so vertical adjustment then horizontal

- but my understanding of the PROJ <6 code is that it does in the reverse order.

playing with cs2cs & cct, I can actually confirm this. There's a difference in vertical values of the order of a few millimeters, whereas if I modify the pipeline to do horizontal then vertical, I get the PROJ < 6 results with a sub-millimeter precision.

I didn’t realise this, but it seems you are right. It looks like the behaviour of PROJ has changed. Although I like the new order (first vertical then horizontal) better, it is very confusing that this has changed. I wonder if this was a bug, or an intended improvement?

> Who is right, who is wrong:  I give up. At that point we probably just need to accept the discrepancy as nobody raised it. Bottom line: use cct instead of cs2cs to unambiguously control the order of operations.

The procedure in both orders are a valid approach, you just need a slightly different geoid grid file for each.

> Actually I found an older thread about that : https://lists.osgeo.org/pipermail/proj/2019-May/008526.html and it already gave me tons of headaches...

I wonder how to prevent other less-experienced users from such headaches.

> using the EPSG:7415 code gives the expected result (but I now see that it is because it uses the nlgeo2018 grid)

$ echo 53 6 0 | PROJ_NETWORK=ON PROJ_LIB=data bin/cs2cs -f "%.4f" EPSG:4937 EPSG:7415
196139.4316    557179.0261 -41.6949

whereas using the PROJ.4 string for EPSG:7415 results (in PROJ >= 6) in the reverse order for horizontal & vertical adjustments

$ echo 53 6 0 | PROJ_NETWORK=ON PROJ_LIB=data bin/cs2cs -f "%.4f" 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 +type=crs"
196139.4316    557179.0261 -41.6982

So, if we accept the change of behaviour in PROJ >=6, we only need the naptrans2008 and 2018 grids for PROJ <6 users. PROJ >=6 users can use the nlgeo2008 and 2018 grids instead.
> On reflection, it might be better to have the GeoTIFF SRS be the interpolation CRS (like in the VERTICAL_TO_VERTICAL case, but in that case there's no other sane choice), an explicit source_crs_code=7931 and interpolation_crs_wkt=....    The https://proj.org/specifications/geodetictiffgrids.html<https://secure-web.cisco.com/16SSmij_qIqbZT7vkQpUgv3KX-6ecNxzEzEXIoOqSbC8wsdbTwfrbE_FHii5vpGRXpy7p7mDJPOV9wwZH5frOT8wENoH3KeWcaJNKsvOcWR102emBb5lTF8g_1Uy7RktJ2UUmj3z6nwbZkOZd3WQ6QeCVKSIRuOgLOe2qshPpMN0mjpyNmZCCOtzD-WFnTN9GJpjcLhgy0uPONd4jA9BqtfzV3_rO7ZZlnOt2NVivF9vXtsez9PnwTiLsFYEcz-qUn4KR6xmRX7nIWT4NeqSkOr2Nx-PEGX3Xx45MPDXKpKzAO5l7xVmeKxcee4UBtr2x8EfeJmv6rPQ8qjX_kOFdafhFPEwhG30LbhIKo8taaFs/https%3A%2F%2Fproj.org%2Fspecifications%2Fgeodetictiffgrids.html> page should be amended to document that, and the validation script too.  And perhaps the conversion script too if that's thought to be a not-so-rare use case.

That’s what I was thinking. However, you won’t need a  interpolation_crs_wkt=.... in that case, as it would be the same as the GeoTIFF SRS, right? Since the naptrans2008 and 2018 grids of the Netherlands are the only use case, and these are only needed for PROJ <6, then we might not need GeoTIFFs for these grids, I think. Maybe, this option could be removed entirely?

> I guess that's OK, adding just a bit more confusion to an already super confusing topic

Good point. I’ll reconsider if I find it worth the additional confusion.

> Note: I'm not sure if it is a formatting issue when writing the email, but the value of the interpolation_crs_wkt item is wrong as it lacks double quoting around strings.
Thanks!

Regards, Jochem



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/20220131/65894da8/attachment-0001.html>


More information about the PROJ mailing list