OSGeo / PROJ-data

Repository for proj datum grids (for use by PROJ 7 or later)
Other
72 stars 33 forks source link

Add grids for DKLAT, DKMSL and GLLMSL #121

Closed kbevers closed 5 months ago

kbevers commented 5 months ago

Transformation grids added to EPSG recently.

rouault commented 5 months ago

The grids should also be mentionned in copyright_and_licenses.csv

kbevers commented 5 months ago

The grids should also be mentionned in copyright_and_licenses.csv

Corrected

rouault commented 5 months ago

For DKLAT and DKMSL grids, I see the metadata of the grid register them as transformations from ETRS89 to a vertical CRS of depth type (so with vertical axis = down). This is I believe the first time we use that, and was slightly confusing at the first sight.

A similar case was for the Norway no_kv_CD_above_Ell_ETRS89_v2023b.tif grid that is ultimately used also with a vertical CRS of depth type but the metadata of the grid indicates a transformation to the CD Norway height CRS (vertical axis = up), so the grid itself is registered in the usual way

For DKLAT and DKMSL, this is probably not that much as an issue in practice since AFAICS, the values in the grid are still geoid undulations, and "projinfo -k operation EPSG:10564" seems to return a pipeline that does the vertical axis inversion, and lead to correct values (I assume you've tested cs2cs with EPSG codes?)

Perhaps https://proj.org/en/9.4/specifications/geodetictiffgrids.html should be clarified for type=geoid_undulation to indicate that: for clarity a target vertical CRS is of height type is preferred, but that a target vertical CRS of depth type may be accepted, with the understanding that the values in the grid are actually to transform to a related vertical CRS of height type. Although that feels hacky...

kbevers commented 5 months ago

the values in the grid are still geoid undulation

They aren't, but they are used in similar manner. And for the MSL also close to the geoid but not quite the geoid. Both DKLAT and DKMSL are meant for offshore use where the chosen axis makes sense. Basically the use case is to be able to measure the depth to the sea floor using GNSS and an echo sounder. This has previosly been done using the DVR90 geoid model and manually(ish) converting to depths instead of heights.

I did struggle a bit to find the best way to convert these to GTG. The options seemed to be either geoid undulation or hydroid height. This would best be described as a "hydroid depth". As far as I remember I did try to use the hydroid option but it caused another problem (which I can't remember on top of my head - I did the initial work a while a go).

rouault commented 5 months ago

I'm a bit confused by your explanation. To clarify: what is the formula to apply to transform from ETRS 89 height (H) to DKMSL/DKLAT depth (d) given the grid value (g) ?

kbevers commented 5 months ago

From the documentation of DKLAT (only available in Danish, unfortunately):

image

D is depth, L is the LAT-value from the grid and h is the ellipsoidal height. DKMSL is similar.

Here's a relevant section from my gie tests:

# DKLAT2022 (Ellipsoid height -> DKLAT depth)
# -------------------------------------------------------------------------------
operation   +proj=pipeline
                +step +proj=vgridshift +grids=./DKLAT/dklat_2022.tif +multiplier=-1
                +step +proj=axisswap +order=1,2,-3,4
tolerance   0.1 mm

# Aarhus Bugt
accept 10.3052 56.1434 24.0789
expect 10.3052 56.1434 14.3140

# Venø Bugt
accept  8.6826 56.5277 33.2506
expect  8.6826 56.5277  6.0901

# Drogden Rende
accept 12.7036 55.5887 27.0668
expect 12.7036 55.5887  8.7467

# Horns Rev
accept  7.9870 55.5135 36.1733
expect  7.9870 55.5135  3.5649

PROJ performs the transformation as I expect it to (although the pipeline is more complicated than it could be):

PROJ_DATA=data ./bin/projinfo -s EPSG:4937 -t EPSG:10550 -o PROJ
Candidate operations found: 1
-------------------------------------
Operation No. 1:

DERIVED_FROM(EPSG):10563, ETRS89 to DKLAT(2023) depth (1), 0.1 m, Denmark - offshore.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=vgridshift +grids=dk_sdfi_dklat_2023.tif +multiplier=1
  +step +proj=axisswap +order=1,2,-3
  +step +proj=unitconvert +xy_in=rad +

# same results as the first gie test case at Aarhus Bugt
echo 56.1434 10.3052 24.0789 | PROJ_DATA=data ./bin/cs2cs EPSG:4937 EPSG:10550    
56.14   10.31 14.31

I believe everything is in order but in case it isn't I'm happy to try to clarify. Note though, I'm going on vacation tomorrow morning so I might not answer this for a while (will be back at the latest the 20th). Feel free to merge if there isn't any issues.

rouault commented 5 months ago

ok, thanks for the clarifications. so indeed the value in the grid behaves similarly to the geoid undulation. We'd probably need to clarify the GTG spec with that use case, but this PR looks good by itself