GDAL 3.0.1: coordinate transformation with vertical shift

[and-marsh]

Expected behavior and actual behavior.

Hello, i'm trying to understand how i can properly work with geoids in GDAL 3. I have a next transformation:
From EPSG:4979(WGS84 3D) transform to EPSG2162(Sierra Leone 1968 / UTM zone 29N) with vertical shift in accordance with EPSG:3855(EGM2008 heigh). I expected, that first, coordinates will be projected from (WGS84) to (Sierra Leone 1968 / UTM zone 29N) and then vertical shift will be applied.
It woks this way in GDAL 2.4.0 with PROJ 5:

/gdaltransform --debug on -s_srs EPSG:4979 -t_srs EPSG:2162+3855
OGRCT: Source: +proj=longlat +datum=WGS84 +no_defs
OGRCT: Target: +proj=utm +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +units=m +geoidgrids=egm08_25.gtx +vunits=m +no_defs
OGRCT: Source: +proj=utm +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +units=m +geoidgrids=egm08_25.gtx +vunits=m +no_defs
OGRCT: Target: +proj=longlat +datum=WGS84 +no_defs
-11 7
279062.796759759 774121.565647936 -63.4349140766573

gdaltransform --debug on -s_srs EPSG:4979 -t_srs EPSG:2162                                  130 ↵
OGRCT: Source: +proj=longlat +datum=WGS84 +no_defs
OGRCT: Target: +proj=utm +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +units=m +no_defs
OGRCT: Source: +proj=utm +zone=29 +ellps=clrk80 +towgs84=-88,4,101,0,0,0,0 +units=m +no_defs
OGRCT: Target: +proj=longlat +datum=WGS84 +no_defs
-11 7
279062.796759759 774121.565647936 -32.7834358848631

gdaltransform --debug on -s_srs EPSG:4979 -t_srs EPSG:4326+3855                             130 ↵
OGRCT: Source: +proj=longlat +datum=WGS84 +no_defs
OGRCT: Target: +proj=longlat +datum=WGS84 +geoidgrids=egm08_25.gtx +vunits=m +no_defs
OGRCT: Source: +proj=longlat +datum=WGS84 +geoidgrids=egm08_25.gtx +vunits=m +no_defs
OGRCT: Target: +proj=longlat +datum=WGS84 +no_defs
-11 7
-11 7 -30.652172088623

It looks like vertical shift was applied to the coordinates obtained after projection.
But in GDAL 3 i have next results:

gdaltransform --debug on -s_srs EPSG:4979 -t_srs EPSG:2162+3855                            130 ↵
-11 7
OGRCT: Selecting transformation +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=egm08_25.gtx +multiplier=1 +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=-88 +y=4 +z=101 +step +inv +proj=cart +a=6378249.145 +rf=293.465 +step +proj=pop +v_3 +step +proj=longlat +a=6378249.145 +rf=293.465 +step +proj=utm +zone=29 +a=6378249.145 +rf=293.465 (WGS 84 to EGM2008 height (1) + Inverse of Sierra Leone 1968 to WGS 84 (1) + UTM zone 29N)
279062.796869158 774121.565475025 -30.652172088623

gdaltransform --debug on -s_srs EPSG:4979 -t_srs EPSG:2162                                 130 ↵
-11 7
279062.796807871 774121.565603866 -32.7820108355954

gdaltransform --debug on -s_srs EPSG:4979 -t_srs EPSG:4326+3855                            130 ↵
-11 7
-11 7 -30.652172088623

Here, if i understand correctly, first was applied vertical shift, then coordinates was projected. It described in piplene: WGS 84 to EGM2008 height (1) + Inverse of Sierra Leone 1968 to WGS 84 (1) + UTM zone 29N
I found out that GDAL 3 always applies vertical shift to source datum, not to the target datum.

Can you explain please is it correct to use such transformation? What result is more correct?

Operating system

macOS Mojave 10.14.6

GDAL version and provenance

GDAL 3.0.1, PROJ 6.2.0

[rouault]

Regarding if horizontal transformation then vertical or vertical then horizontal must be applied, I haven't found yet a definitive answer on this. For that particular case thoug, EPSG:4979 to EPSG:2162+3855, applying the vertical transformation first using egm08_25.gtx grid makes sense since this grid is to transform between EGM2008 height and WGS 84 ellipsoidal height, and I believe that the grid assumes horizontal coordinates in WGS 84, so I think the order used here, vertical shift then horizontal shift, makes sense.

What result is more correct?

Here the main difference between GDAL 2.4 (-63.4m) and GDAL 3.0 (-32.7m) mostly comes from the fact that PROJ4/5 transformation was also applying a transformation on Z when doing conversion between the WGS84 and Sierra Leone 1968 datum going through geocentric coordinates, and during PROJ 6 development, we deemed that this transformation didn't make sense when the CRS was a compound made of a horizontal + vertical CRS. So GDAL 3 + PROJ6 result is certainly closer to the "truth" to GDAL 2.x + PROJ4/5

So all in all I think GDAL 3 results are OK

[and-marsh]

@rouault Thank you for your answer!

Some more tests for GDAL 3 vertical shift:
I tried to transform from EPSG:4979(WGS84 3D) to EPSG:3040(ETRS89 / UTM zone 28N (N-E))+5732( Belfast heigh)
EPSG:5732 is related to ETRS89 datum and EPSG:3040 too.

testepsg "EPSG:3040+5732"
Validate Succeeds.
WKT[EPSG:3040+5732] =
COMPD_CS["ETRS89 / UTM zone 28N (N-E) + Belfast height",
    PROJCS["ETRS89 / UTM zone 28N (N-E)",
        GEOGCS["ETRS89",
            DATUM["European_Terrestrial_Reference_System_1989",
                SPHEROID["GRS 1980",6378137,298.257222101,
                    AUTHORITY["EPSG","7019"]],
                TOWGS84[0,0,0,0,0,0,0],
                AUTHORITY["EPSG","6258"]],
            PRIMEM["Greenwich",0,
                AUTHORITY["EPSG","8901"]],
            UNIT["degree",0.0174532925199433,
                AUTHORITY["EPSG","9122"]],
            AUTHORITY["EPSG","4258"]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",0],
        PARAMETER["central_meridian",-15],
        PARAMETER["scale_factor",0.9996],
        PARAMETER["false_easting",500000],
        PARAMETER["false_northing",0],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Northing",NORTH],
        AXIS["Easting",EAST],
        AUTHORITY["EPSG","3040"]],
    VERT_CS["Belfast height",
        VERT_DATUM["Belfast Lough",2005,
            AUTHORITY["EPSG","5131"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Gravity-related height",UP],
        AUTHORITY["EPSG","5732"]]]

...

PROJ.4 rendering of [EPSG:3040+5732] = +proj=utm +zone=28 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs

I understand that WGS84 and ETRS89 are almost the same and i expect, that vertical shift will be applied to this transformation but it's not :

gdaltransform --debug on -s_srs EPSG:4979 -t_srs EPSG:3040+5732                            130 ↵
-7 54
1023866.07679486 6013191.39274665 0

So, i tried to replace EPSG:4979(WGS84(3D)) with EPSG:4937(ETRS89(3D)) and it worked fine:

gdaltransform --debug on -s_srs EPSG:4937 -t_srs EPSG:3040+5732
-7 54
OGRCT: Selecting transformation +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=OSGM15_Belfast.gtx +multiplier=1 +step +proj=utm +zone=28 +ellps=GRS80 +step +proj=axisswap +order=2,1 (ETRS89 to Belfast height (2) + ETRS89 to WGS 84 (1) + Inverse of ETRS89 to WGS 84 (1) + UTM zone 28N)
1023866.07679486 6013191.39274665 -57.0370325025414

Also i tried to replace target SRS with proj4 string, that i got from testepsg EPSG:3040+5732 and added geoid file +geoidgrids=OSGM15_Belfast.gtx:

gdaltransform --debug on -s_srs EPSG:4937 -t_srs "+proj=utm +zone=28 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs +geoidgrids=OSGM15_Belfast.gtx"
-7 54
1023866.07679486 6013191.39274665 -57.0370325025414

./gdaltransform --debug on -s_srs EPSG:4979 -t_srs "+proj=utm +zone=28 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +vunits=m +no_defs +geoidgrids=OSGM15_Belfast.gtx"
-7 54
1023866.07679486 6013191.39274665 -57.0370325025414

And this way it work fine for both: EPSG:4979 and EPSG:4973.
So, i found during debug, that vertical shift is applied only if source datum and vertical shift datum are exactly the same or if we force it with proj4 +geoidgrids=OSGM15_Belfast.gtx! In first case related file OSGM15_Belfast.gtx was found for target vertical SRS with DatabaseContext::lookForGridAlternative (.../proj-6.2.0/src/iso19111/factory.cpp:956) in grid_alternatives table. But as i said it works only if source datum and vertical shift datum are the same.

Is it bug or not?
I'm just trying to understand, how can i use it right way. I wanted to use it to transform from EPSG:4979(always) to any other system.

Thank you in advance for your answer!

[rouault]

Regarding EPSG:4979 to EPSG:3040+5732, you've exactly hit the work I've done this week in OSGeo/PROJ#1608

With that work applied, then it works out of the box:

$ src/projinfo -s EPSG:4979 -t EPSG:3040+5732 --spatial-test intersects
Candidate operations found: 4
-------------------------------------
Operation n°1:

unknown id, ETRS89 to Belfast height (2) + Inverse of ETRS89 to WGS 84 (1) + UTM zone 28N, 1.014 m, UK - Northern Ireland - onshore, at least one grid missing

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=OSGM15_Belfast.gtx +multiplier=1 +step +proj=utm +zone=28 +ellps=GRS80 +step +proj=axisswap +order=2,1
[...]

[and-marsh]

That sounds great! Thank you! I will wait for new version

[rouault]

Closing this as it is only PROJ business.