<div dir="ltr">Is this a special case for EPSG:5773, or is it a bad idea to express the vertical datum via EPSG codes in general?<div><br></div><div>Patrick</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Jun 16, 2020 at 2:34 PM Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><u></u>
<div style="font-family:monospace;font-size:9pt;font-weight:400;font-style:normal">
<p style="margin:0px;text-indent:0px">On mardi 16 juin 2020 20:22:57 CEST Comer, Alex wrote:</p>
<p style="margin:0px;text-indent:0px">> Hello, I'm having trouble using gdalwarp to adjust the heights in a DEM from</p>
<p style="margin:0px;text-indent:0px">> EGM96 heights to ellipsoidal heights.</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> There is an example in the GDAL documentation</p>
<p style="margin:0px;text-indent:0px">> (<a href="https://gdal.org/programs/gdalwarp.html#examples" target="_blank">https://gdal.org/programs/gdalwarp.html#examples</a><<a href="https://urldefense.proofp" target="_blank">https://urldefense.proofp</a></p>
<p style="margin:0px;text-indent:0px">> <a href="http://oint.com/v2/url?u=https-3A__slack-2Dredir.net_link-3Furl-3Dhttps-253A-252F-2" target="_blank">oint.com/v2/url?u=https-3A__slack-2Dredir.net_link-3Furl-3Dhttps-253A-252F-2</a></p>
<p style="margin:0px;text-indent:0px">> 52Fgdal.org-252Fprograms-252Fgdalwarp.html-2523examples-26v-3D3&d=DwMFAw&c=q</p>
<p style="margin:0px;text-indent:0px">> v9gmtOiXdJCj9gMlMKtaw&r=c5qG9R_C_Uk15wO7YCz0Kmxn1feCJcKpkdWjm6OrnHA&m=LBJGDE</p>
<p style="margin:0px;text-indent:0px">> isk_mehbQp1YF7U2vnpaqKAyB1WlXNejSAiHc&s=SPJdv5gp9MY8S77hKdINBORe-qDXH6QvCJa4</p>
<p style="margin:0px;text-indent:0px">> mrKxl9c&e=> , near the end of the page):</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> * To transform a DEM from geoid elevations (using EGM96) to WGS84</p>
<p style="margin:0px;text-indent:0px">> ellipsoidal heights:</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> New in version 2.2.</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> gdalwarp -overwrite in_dem.tif out_dem.tif -s_srs EPSG:4326+5773 -t_srs</p>
<p style="margin:0px;text-indent:0px">> EPSG:4979</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> When I run this command using GDAL 3.1.0 and PROJ 7.0.1, the output is</p>
<p style="margin:0px;text-indent:0px">> unchanged from the input. My expectation is that the heights would be</p>
<p style="margin:0px;text-indent:0px">> adjusted in the output file.</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> If I set CPL_DEBUG=ON, I see the message:</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> GDALWARP: Source SRS is a compound CRS but lacks +geoidgrids</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> Using gdalsrsinfo to examine the SRS "EPSG:4326+5773" shows a WKT that</p>
<p style="margin:0px;text-indent:0px">> notably does not refer to +geoidgrids in the PROJ string:</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> $ gdalsrsinfo EPSG:4326+5773</p>
<p style="margin:0px;text-indent:0px">> PROJ.4 : +proj=longlat +datum=WGS84 +vunits=m +no_defsOGC WKT2:2018 :</p>
<p style="margin:0px;text-indent:0px">> COMPOUNDCRS["WGS 84 + EGM96 height",</p>
<p style="margin:0px;text-indent:0px">> GEOGCRS["WGS 84",</p>
<p style="margin:0px;text-indent:0px">> DATUM["World Geodetic System 1984",</p>
<p style="margin:0px;text-indent:0px">> ELLIPSOID["WGS 84",6378137,298.257223563,</p>
<p style="margin:0px;text-indent:0px">> LENGTHUNIT["metre",1]]],</p>
<p style="margin:0px;text-indent:0px">> PRIMEM["Greenwich",0,</p>
<p style="margin:0px;text-indent:0px">> ANGLEUNIT["degree",0.0174532925199433]],</p>
<p style="margin:0px;text-indent:0px">> CS[ellipsoidal,2],</p>
<p style="margin:0px;text-indent:0px">> AXIS["geodetic latitude (Lat)",north,</p>
<p style="margin:0px;text-indent:0px">> ORDER[1],</p>
<p style="margin:0px;text-indent:0px">> ANGLEUNIT["degree",0.0174532925199433]],</p>
<p style="margin:0px;text-indent:0px">> AXIS["geodetic longitude (Lon)",east,</p>
<p style="margin:0px;text-indent:0px">> ORDER[2],</p>
<p style="margin:0px;text-indent:0px">> ANGLEUNIT["degree",0.0174532925199433]],</p>
<p style="margin:0px;text-indent:0px">> USAGE[</p>
<p style="margin:0px;text-indent:0px">> SCOPE["unknown"],</p>
<p style="margin:0px;text-indent:0px">> AREA["World"],</p>
<p style="margin:0px;text-indent:0px">> BBOX[-90,-180,90,180]],</p>
<p style="margin:0px;text-indent:0px">> ID["EPSG",4326]],</p>
<p style="margin:0px;text-indent:0px">> VERTCRS["EGM96 height",</p>
<p style="margin:0px;text-indent:0px">> VDATUM["EGM96 geoid"],</p>
<p style="margin:0px;text-indent:0px">> CS[vertical,1],</p>
<p style="margin:0px;text-indent:0px">> AXIS["gravity-related height (H)",up,</p>
<p style="margin:0px;text-indent:0px">> LENGTHUNIT["metre",1]],</p>
<p style="margin:0px;text-indent:0px">> USAGE[</p>
<p style="margin:0px;text-indent:0px">> SCOPE["unknown"],</p>
<p style="margin:0px;text-indent:0px">> AREA["World"],</p>
<p style="margin:0px;text-indent:0px">> BBOX[-90,-180,90,180]],</p>
<p style="margin:0px;text-indent:0px">> ID["EPSG",5773]]]</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> A colleague using GDAL v2.4.1 had a different result. Note the reference to</p>
<p style="margin:0px;text-indent:0px">> +geoidgrids=egm95_15.gtx and the PROJ4_GRIDS extension:</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> $ gdalsrsinfo EPSG:4326+5773</p>
<p style="margin:0px;text-indent:0px">> PROJ.4 : +proj=longlat +datum=WGS84 +geoidgrids=egm96_15.gtx +vunits=m</p>
<p style="margin:0px;text-indent:0px">> +no_defsOGC WKT : COMPD_CS["WGS 84 + EGM96 geoid height",</p>
<p style="margin:0px;text-indent:0px">> GEOGCS["WGS 84",</p>
<p style="margin:0px;text-indent:0px">> DATUM["WGS_1984",</p>
<p style="margin:0px;text-indent:0px">> SPHEROID["WGS 84",6378137,298.257223563,</p>
<p style="margin:0px;text-indent:0px">> AUTHORITY["EPSG","7030"]],</p>
<p style="margin:0px;text-indent:0px">> AUTHORITY["EPSG","6326"]],</p>
<p style="margin:0px;text-indent:0px">> PRIMEM["Greenwich",0,</p>
<p style="margin:0px;text-indent:0px">> AUTHORITY["EPSG","8901"]],</p>
<p style="margin:0px;text-indent:0px">> UNIT["degree",0.0174532925199433,</p>
<p style="margin:0px;text-indent:0px">> AUTHORITY["EPSG","9122"]],</p>
<p style="margin:0px;text-indent:0px">> AUTHORITY["EPSG","4326"]],</p>
<p style="margin:0px;text-indent:0px">> VERT_CS["EGM96 geoid height",</p>
<p style="margin:0px;text-indent:0px">> VERT_DATUM["EGM96 geoid",2005,</p>
<p style="margin:0px;text-indent:0px">> EXTENSION["PROJ4_GRIDS","egm96_15.gtx"],</p>
<p style="margin:0px;text-indent:0px">> AUTHORITY["EPSG","5171"]],</p>
<p style="margin:0px;text-indent:0px">> UNIT["metre",1,</p>
<p style="margin:0px;text-indent:0px">> AUTHORITY["EPSG","9001"]],</p>
<p style="margin:0px;text-indent:0px">> AXIS["Up",UP],</p>
<p style="margin:0px;text-indent:0px">> AUTHORITY["EPSG","5773"]]]</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> Indeed, with GDAL 3.1.0 & PROJ 7.0.1 I am able to achieve the desired result</p>
<p style="margin:0px;text-indent:0px">> if I use a PROJ string rather than the EPSG:horiz+vertical code. For</p>
<p style="margin:0px;text-indent:0px">> example something like:</p>
<p style="margin:0px;text-indent:0px">> </p>
<p style="margin:0px;text-indent:0px">> gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_defs</p>
<p style="margin:0px;text-indent:0px">> +geoidgrids=egm96_15.gtx" -t_srs "+proj=longlat +datum=WGS84 +no_def"</p>
<p style="margin:0px;text-indent:0px">> input.tif output.tif</p>
<p style="margin:0px;text-indent:0px"> </p>
<p style="margin:0px;text-indent:0px">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.</p>
<p style="margin:0px;text-indent:0px"> </p>
<p style="margin:0px;text-indent:0px">Even</p>
<p style="margin:0px;text-indent:0px"> </p>
<p style="margin:0px;text-indent:0px">-- </p>
<p style="margin:0px;text-indent:0px">Spatialys - Geospatial professional services</p>
<p style="margin:0px;text-indent:0px"><a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a></p></div>_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/gdal-dev" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a></blockquote></div>