[gdal-dev] Source SRS is a compound CRS but lacks +geoidgrids

[Comer, Alex]

Hello, I'm having trouble using gdalwarp to adjust the heights in a DEM from EGM96 heights to ellipsoidal heights.

There is an example in the GDAL documentation (https://gdal.org/programs/gdalwarp.html#examples<https://urldefense.proofpoint.com/v2/url?u=https-3A__slack-2Dredir.net_link-3Furl-3Dhttps-253A-252F-252Fgdal.org-252Fprograms-252Fgdalwarp.html-2523examples-26v-3D3&d=DwMFAw&c=qv9gmtOiXdJCj9gMlMKtaw&r=c5qG9R_C_Uk15wO7YCz0Kmxn1feCJcKpkdWjm6OrnHA&m=LBJGDEisk_mehbQp1YF7U2vnpaqKAyB1WlXNejSAiHc&s=SPJdv5gp9MY8S77hKdINBORe-qDXH6QvCJa4mrKxl9c&e=> , near the end of the page):

  *   To transform a DEM from geoid elevations (using EGM96) to WGS84 ellipsoidal heights:

New in version 2.2.

gdalwarp -overwrite in_dem.tif out_dem.tif -s_srs EPSG:4326+5773 -t_srs EPSG:4979

When I run this command using GDAL 3.1.0 and PROJ 7.0.1, the output is unchanged from the input. My expectation is that the heights would be adjusted in the output file.

If I set CPL_DEBUG=ON, I see the message:

GDALWARP: Source SRS is a compound CRS but lacks +geoidgrids

Using gdalsrsinfo to examine the SRS "EPSG:4326+5773" shows a WKT that notably does not refer to +geoidgrids in the PROJ string:

$ gdalsrsinfo EPSG:4326+5773
PROJ.4 : +proj=longlat +datum=WGS84 +vunits=m +no_defsOGC WKT2:2018 :
COMPOUNDCRS["WGS 84 + EGM96 height",
    GEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        CS[ellipsoidal,2],
            AXIS["geodetic latitude (Lat)",north,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433]],
            AXIS["geodetic longitude (Lon)",east,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433]],
        USAGE[
            SCOPE["unknown"],
            AREA["World"],
            BBOX[-90,-180,90,180]],
        ID["EPSG",4326]],
    VERTCRS["EGM96 height",
        VDATUM["EGM96 geoid"],
        CS[vertical,1],
            AXIS["gravity-related height (H)",up,
                LENGTHUNIT["metre",1]],
        USAGE[
            SCOPE["unknown"],
            AREA["World"],
            BBOX[-90,-180,90,180]],
        ID["EPSG",5773]]]



A colleague using GDAL v2.4.1 had a different result. Note the reference to +geoidgrids=egm95_15.gtx and the PROJ4_GRIDS extension:



$ gdalsrsinfo EPSG:4326+5773
PROJ.4 : +proj=longlat +datum=WGS84 +geoidgrids=egm96_15.gtx +vunits=m +no_defsOGC WKT :
COMPD_CS["WGS 84 + EGM96 geoid height",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    VERT_CS["EGM96 geoid height",
        VERT_DATUM["EGM96 geoid",2005,
            EXTENSION["PROJ4_GRIDS","egm96_15.gtx"],
            AUTHORITY["EPSG","5171"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Up",UP],
        AUTHORITY["EPSG","5773"]]]



Indeed, with GDAL 3.1.0 & PROJ 7.0.1 I am able to achieve the desired result if I use a PROJ string rather than the EPSG:horiz+vertical code. For example something like:

gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=egm96_15.gtx" -t_srs "+proj=longlat +datum=WGS84 +no_def" input.tif output.tif

It seems that gdalwarp is having trouble parsing the +5773 (vertical) part of the EPSG code into the required +geoidgrids=egm96_15.gtx in the PROJ string.

I wondered whether GDAL was having trouble finding the PROJ files, but GDAL seems to find proj.db okay (this is my output with PROJ_DEBUG=3 and CPL_DEBUG=ON)...

$ gdalsrsinfo EPSG:4326+5773
gdalsrsinfo: got arg #1 : [EPSG:4326+5773]
gdalsrsinfo: trying to open with GDAL
gdalsrsinfo: trying to get SRS from user input [EPSG:4326+5773]
PROJ: pj_open_lib(proj.db): call fopen(/usr/share/proj/proj.db) - succeeded
gdalsrsinfo: got SRS from user input
gdalsrsinfo: bGotSRS: 1 bValidate: 0 pszOutputType: default bPretty: 1gdalsrsinfo: PrintSRS( oSRS, proj4, 1, 1 )PROJ.4 : +proj=longlat +datum=WGS84 +vunits=m +no_defsgdalsrsinfo: PrintSRS( oSRS, wkt2, 1, 1 )OGC WKT2:2018 :
COMPOUNDCRS["WGS 84 + EGM96 height",
    GEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        CS[ellipsoidal,2],
            AXIS["geodetic latitude (Lat)",north,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433]],
            AXIS["geodetic longitude (Lon)",east,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433]],
        USAGE[
            SCOPE["unknown"],
            AREA["World"],
            BBOX[-90,-180,90,180]],
        ID["EPSG",4326]],
    VERTCRS["EGM96 height",
        VDATUM["EGM96 geoid"],
        CS[vertical,1],
            AXIS["gravity-related height (H)",up,
                LENGTHUNIT["metre",1]],
        USAGE[
            SCOPE["unknown"],
            AREA["World"],
            BBOX[-90,-180,90,180]],
        ID["EPSG",5773]]]

I expected that we should be seeing a +geoidgrids=egm96_15.gtx
at the end of: PROJ.4 : +proj=longlat +datum=WGS84 +vunits=m +no_defs

Am I doing something wrong, or misunderstanding how things work now? Did something change with GDAL 3.1.0 or perhaps it should not be used with PROJ 7.0.1 ?
I would appreciate any insight anyone may have.
Thanks,
Alex
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200616/4355171b/attachment.html>