depth to depth or two CRS with orientation down issue?

[barry-gallagher]

Problem description

In my custom proj.db, attached, I am trying to connect three compound CRS depth objects with a concatenated operation of two grid_transformations. My concatenated_operation of the heights seems to work. I then tried just calling the grid_transformation and got unexpected results. I may have made the depth operations incorrectly but swapping the source and target has no effect which seems wrong even if I had something set wrong.

Example of problem

proj.zip

// Reversing the source and target has no effect, both CRS are coordinate system EPSG:6498 (down)
echo 48.3 -123 10.0 | cs2cs --only-best --no-ballpark NOAA:8506 NOAA:8507
48.30   -123.00 -9.17

echo 48.3 -123 10.0 | cs2cs --only-best --no-ballpark NOAA:8507 NOAA:8506
48.30   -123.00 -9.17

// Reversing the source and target on the height CRS looks ok
echo 48.3 -123 10.0 | cs2cs --only-best --no-ballpark NOAA:8511 NOAA:8512
48.30   -123.00 10.83

echo 48.3 -123 10.0 | cs2cs --only-best --no-ballpark NOAA:8512 NOAA:8511
48.30   -123.00 9.17

// this is the concatenated height,   2.08 is the expected result of adding 0.83 and 1.24
echo 48.3 -123 0.0 | cs2cs --only-best --no-ballpark NOAA:8550 NOAA:8512
48.30   -123.00 2.08

// this is the concatenated depth  -0.41  the signs are getting switched in the transform
echo 48.3 -123 0.0 | cs2cs --only-best --no-ballpark NOAA:8546 NOAA:8507
48.30   -123.00 -0.41

projinfo --hide-ballpark -s NOAA:8546 -t NOAA:8507 --spatial-test intersects | more
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of 'NAD83(FBN) + LMSL (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth to NAD83(FBN) + NAVD88 depth' + Inverse of 'NAD83(FBN) + MHHW (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth to NAD83(FBN) + LMSL (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth', 0.346410161513776 m, 17203 Washington - Strait of Juan de Fuca Local Mean Sea Level 1983-2001

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=axisswap +order=1,2,-3
  +step +proj=vgridshift
        +grids=us_noaa_nos_NAD83(FBN)_LMSL_NAVD88_(WAjuandefuca01_vdatum_3.4_20141224_1983-2001).tif
        +multiplier=1
  +step +proj=axisswap +order=1,2,-3
  +step +proj=vgridshift
        +grids=us_noaa_nos_NAD83(FBN)_MHHW_LMSL_(WAjuandefuca01_vdatum_3.4_20141224_1983-2001).tif
        +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

Expected Output

Depth and height operations would hopefully be equal but opposite sign.

// Reversing the source and target has no effect, both CRS are coordinate system EPSG:6498 (down)
echo 48.3 -123 10.0 | cs2cs --only-best --no-ballpark NOAA:8506 NOAA:8507
48.30   -123.00 9.17

echo 48.3 -123 10.0 | cs2cs --only-best --no-ballpark NOAA:8507 NOAA:8506
48.30   -123.00 10.83

// Reversing the source and target on the height CRS looks ok
echo 48.3 -123 10.0 | cs2cs --only-best --no-ballpark NOAA:8511 NOAA:8512
48.30   -123.00 10.83

echo 48.3 -123 10.0 | cs2cs --only-best --no-ballpark NOAA:8512 NOAA:8511
48.30   -123.00 9.17

// this is the concatenated height,   2.08 is the expected result of adding 0.83 and 1.24
echo 48.3 -123 0.0 | cs2cs --only-best --no-ballpark NOAA:8550 NOAA:8512
48.30   -123.00 2.08

// this is the concatenated depth  -0.41  the signs are getting switched in the transform
echo 48.3 -123 0.0 | cs2cs --only-best --no-ballpark NOAA:8546 NOAA:8507
48.30   -123.00 -2.08

Environment Information

• PROJ version (9.4)
• Operation System Information (Windows 11)

Installation method

• conda

[rouault]

I've investigated that, and I believe the root cause if that operations such as EPSG:11488 are abusing the "Geog3D to Geog2D+GravityRelatedHeight (gtx)" method, in a way it doesn't anticipate. This method should only be used for a Geographic3D CRS to a compound CRS. Here you use 2 compound CRS. The code in PROJ that deals with "Geog3D to Geog2D+GravityRelatedHeight" and "Geog3D to Geog2D+Depth are actually the same, with the only difference that it does a vertical axis inversion prior to applying the geoid model if the target vertical CRS has a down orientation.

I'm not sure how that should be resolved. We are really a in very advanced use case. I believe it would be cleaner for those depth to depth transformation to use EPSG:1084 "Vertical Offset by Grid Interpolation (gtx)" to handle the vertical transformation with vertical axis=UP, and add an explicit height/depth reversal operation before and after (which will likely require extra tricks in fixStepsDirection()...). Or you might want to write custom PROJ based operations with an explicit pipeline

[barry-gallagher]

Thanks @rouault, I actually started with EPSG:1084 doing vertical only transforms but ran into issues with the horizontal so we thought we would enforce that the correct horizontal was used with the tidal datum models we produce. It may be that I didn't have interpolation_crs_code set when I did my initial testing -- I'll retry using EPSG:1084.

You had suggested the pipelines in #3819 when I was asking about doing extra transforms to 3D and I actually had those working. I recently went back to concatenated (seems more object oriented/modular) when we started leaning away from trying to get all our tidal datums to have transforms to a singular 3D pivot. I meant to ask on the mailing list about what the right compound to compound would be since I only saw the EPSG:1084 (vertical only), EPSG:1088 (3d to compound) or EPSG:1115 (hydroid) as the closest options.

I'll go back to the drawing board and test a bit more. Thanks again.

[Rjialceky]

Of particular concern to us is to retain "3D" CRS attribution (whether Geog3D or compound) in source data, through to coordinate expression in certain target CRS connected by reversible transformations.
The "Geog3D to Geog2D+GravityRelatedHeight" (H = h - ζ) and "Geog3D to Geog2D+Depth" (D = ζ - h) [reversible] methods (per EPSG Guidance Note 7-2) will give the correct results using Geog3D or compound referencing if the sense of the parameters in the defined equations are treated consistently, including through concatenation:
(1) 'H', 'h', and 'ζ' are always up-orientation; 'D' is always down-orientation.
(2) 'h' is the vertical component of the "Geog3D" (or compound) source CRS/datum; 'H' and 'D' are associated with the vertical component of the "Geog2D+[GravityRelatedHeight, Depth]" target CRS/datum, resp. This convention remains intact when reversing operations. 'ζ' is the height correction model file.
(3) Concatenated operations involve cross-parameter & height- or depth-orientation mapping between h and that of H and D.

A straightforward example is a (reversible) compound to compound via concatenated operations pair, through a pivot in the Geog3D; e.g., EPSG:10496 "ETRS89 + DVR90(2013) height to ETRS89 + DVR90(2023) height (1)":
In operation 1: "ETRS89 to ETRS89 + DVR90(2013) height (1)", the compound input vertical CRS DVR90(2013) height is assigned as 'H' and the "Geog3D" pivot is h = H + ζ (dvr90_2013.tif).
In operation 2: "ETRS89 to ETRS89 + DVR90(2023) height (1)", the ETRS89 value of h from operation 1 is used to calculate the compound output vertical CRS DVR90(2023) H = h - ζ (dvr90_2023.tif)

Correction of prior NOAA examples (last two use concatenated operations):
Aside: We are also standardizing on ζ in Geodetic TIFF grid (GTG) format, so the "Geog3D" to Geog2D+GravityRelatedHeight and Geog2D+Depth coordinate operation method codes should be EPSG:1124 and EPSG:1128, resp.

NOAA:8506 = NAD83(FBN) + LMSL (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth
to
NOAA:8507 = NAD83(FBN) + MHHW (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth
This is a "Geog3D to Geog2D+Depth" method, but the reverse of the one registered below.
(NOAA:11488)
So the reverse transform is h = ζ - D, but will negate the value of h for depth orientation result.
Only one grid operation is involved with this: ζ = -0.83
Input LMSL depth = +10, means input D = +10
Output MHHW depth means = -h = -(-0.83 - (+10) = +10.83

NOAA:8507 = NAD83(FBN) + MHHW (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth
to
NOAA:8506 = NAD83(FBN) + LMSL (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth
This should be registered as a "Geog3D to Geog2D+Depth" method, so transform is D = ζ - h
(NOAA:11488)
Only one grid operation is involved with this: ζ = -0.83
Input MHHW depth = +10, means input h = -10
Output LMSL depth = D = -0.83 - (-10) = +9.17

NOAA:8511 = NAD83(FBN) + LMSL (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) height
to
NOAA:8512 = NAD83(FBN) + MHHW (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) height
This is a "Geog3D to Geog2D+GravityRelatedHeight" method, but the reverse of the one registered below.
(NOAA:11493)
So the reverse transform is h = H + ζ
Only one grid operation is involved with this: ζ = -0.83
Input LMSL height = +10, means input H = +10
Output MHHW height = h = -0.83 + (+10) = +9.17

NOAA:8512 = NAD83(FBN) + MHHW (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) height
to
NOAA:8511 = NAD83(FBN) + LMSL (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) height
This should be registered as a "Geog3D to Geog2D+GravityRelatedHeight" method, so transform is H = h - ζ
(NOAA:11493)
Only one grid operation is involved with this: ζ = -0.83
Input MHHW height = +10, means input h = +10
Output LMSL height = H = +10 - (-0.83) = +10.83

NOAA:8550 = NAD83(FBN) + NAVD88 height
to
NOAA 8512 = NAD83(FBN) + MHHW (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) height
Two concatenated grid operations of "Geog3D to Geog2D+GravityRelatedHeight" method are involved, reverse of the two registered; so the transform is h = H + ζ
1) NAD83(FBN) + LMSL (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) height to NAD83(FBN) + NAVD88 height
ζ_NOAA:11497 = -1.24
Input NAVD88 height = 0.0 means H = 0.0
Output LMSL height h = 0.0 + (-1.24) = -1.24
2) NAD83(FBN) + MHHW (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) height to NAD83(FBN) + LMSL (ibid.) height
ζ_NOAA:11493 = -0.83
Input LMSL height = -1.24 means H = -1.24
Output MHHW height h = -1.24 + (-0.83) = -2.07

NOAA:8546 = NAD83(FBN) + NAVD88 depth
to
NOAA:8507 = NAD83(FBN) + MHHW (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth
Two concatenated grid operations of "Geog3D to Geog2D+Depth" method are involved, reverse of the two registered; so the transform is h = ζ - D, but will negate the value of h for depth orientation result.
1) NAD83(FBN) + LMSL (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth to NAD83(FBN) + NAVD88 depth
ζ_NOAA:11492 = -1.24
Input NAVD88 depth = 0.0 means D = 0.0
Output LMSL depth means -h = -(-1.24 - 0.0) = +1.24
2) NAD83(FBN) + MHHW (WAjuandefuca01/vdatum_3.4_20141224/1983-2001) depth to NAD83(FBN) + LMSL (ibid.) depth
ζ_NOAA:11488 = -0.83
Input LMSL depth = +1.24 means input D = +1.24
Output MHHW depth means -h = -(-0.83 - (+1.24)) = +2.07