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

[Even Rouault]

On mardi 16 juin 2020 20:22:57 CEST Comer, Alex wrote:
> 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.proofp
> oint.com/v2/url?u=https-3A__slack-2Dredir.net_link-3Furl-3Dhttps-253A-252F-2
> 52Fgdal.org-252Fprograms-252Fgdalwarp.html-2523examples-26v-3D3&d=DwMFAw&c=q
> 
v9gmtOiXdJCj9gMlMKtaw&r=c5qG9R_C_Uk15wO7YCz0Kmxn1feCJcKpkdWjm6OrnHA&m=L
BJGDE
> isk_mehbQp1YF7U2vnpaqKAyB1WlXNejSAiHc&s=SPJdv5gp9MY8S77hKdINBORe-
qDXH6QvCJa4
> mrKxl9c&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

Yes, since GDAL 3 / PROJ 6, this is the expected way to proceed, even if admitedly not ideal 
and could benefit from enhancements, since EPSG:5773 is no longer automatically expanded 
(for "good" reasons) to +geoidgrids=egm96_15.gtx , and gdalwarp in height adjustment 
mode requires an explicit mention of geoidgrids.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200616/3bc6b01e/attachment-0001.html>