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

[Patrick Young]

Is this a special case for EPSG:5773, or is it a bad idea to express the
vertical datum via EPSG codes in general?

Patrick

On Tue, Jun 16, 2020 at 2:34 PM Even Rouault <even.rouault at spatialys.com>
wrote:

> 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=LBJGDE
>
> >
> 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
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200616/38f68a8a/attachment.html>