search

OSGeo / PROJ

PROJ - Cartographic Projections and Coordinate Transformations Library
https://proj.org
Other
1.74k stars 787 forks source link

depth to depth or two CRS with orientation down issue? #4118

Open barry-gallagher opened 6 months ago

barry-gallagher commented 6 months ago

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

Installation method

rouault commented 6 months ago

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 commented 6 months ago

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 commented 6 months ago

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