<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="display:none"><!--P{margin-top:0;margin-bottom:0;} --></style>
</head>
<body dir="ltr" style="font-size:14pt;color:#000000;background-color:#FFFFFF;font-family:Calibri,Arial,Helvetica,sans-serif;">
<p></p>
<div class="_rp_R3 rpHighlightAllClass rpHighlightBodyClass allowTextSelection" style="">
<div style="">
<div class="_rp_S3 ms-font-weight-regular ms-font-color-neutralDark" role="button" tabindex="0">
<div>
<div dir="ltr" style="background-color:white;"><font size="4" face="Calibri,Arial,Helvetica,sans-serif" color="black"><span style="font-size:14pt;background-color:white;" dir="ltr">
<div style="margin-top:0;margin-bottom:0;"> </div>
<div><font face="Calibri,Arial,Helvetica,sans-serif">Hello, I'm having trouble using gdalwarp to adjust the heights in a DEM from EGM96 heights to ellipsoidal heights.</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif"><br>
</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif">
<div>There is an example in the GDAL documentation (<a href="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=" target="_blank">https://gdal.org/programs/gdalwarp.html#examples</a>
, near the end of the page):</div>
<div style="margin:14pt 30pt;">
<ul style="margin-top:14pt;margin-bottom:14pt;">
<li><em>To transform a DEM from geoid elevations (using EGM96) to WGS84 ellipsoidal heights:</em></li></ul>
<div style="margin:14pt 30pt;"><em>New in version 2.2.</em></div>
<pre style="margin-top:14pt;margin-bottom:14pt;"><em>gdalwarp -overwrite in_dem.tif out_dem.tif -s_srs EPSG:4326+5773 -t_srs EPSG:4979</em></pre>
</div>
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.<br>
</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif"><br>
</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif">If I set CPL_DEBUG=ON, I see the message:</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif"><br>
</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif"><font size="-1" face="Courier New">GDALWARP: Source SRS is a compound CRS but lacks +geoidgrids</font></font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif"><br>
</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif">Using gdalsrsinfo to examine the SRS "EPSG:4326+5773" shows a WKT that notably does not refer to +geoidgrids in the PROJ string:
<pre style="margin-top:14pt;margin-bottom:14pt;">$ 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]]]
</pre>
A colleague using GDAL v2.4.1 had a different result. Note the reference to +geoidgrids=egm95_15.gtx and the PROJ4_GRIDS extension:<br>
<pre style="margin-top:14pt;margin-bottom:14pt;"> </pre>
</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif">
<pre style="margin-top:14pt;margin-bottom:14pt;">$ 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"]]]
</pre>
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:<br>
<br>
<font face="Courier New,monospace">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</font></font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif"><br>
It seems that gdalwarp is having trouble parsing the <font size="-1" face="Courier New">
+5773</font> (vertical) part of the EPSG code into the required <font size="-1" face="Courier New">
+geoidgrids=egm96_15.gtx</font> in the PROJ string.</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif"><br>
</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif">I wondered whether GDAL was having trouble finding the PROJ files, but GDAL seems to find
<font size="-1" face="Courier New">proj.db</font> okay (this is my output with PROJ_DEBUG=3 and CPL_DEBUG=ON)...<br>
<div>
<div>
<div>
<div>
<div>
<div>
<pre style="margin-top:14pt;margin-bottom:14pt;">$ 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]]]</pre>
</div>
</div>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>
<div>I expected that we should be seeing a <font size="-1" face="Courier New">+geoidgrids=egm96_15.gtx</font>
</div>
<div>at the end of: <font size="-1" face="Courier New">PROJ.4 : +proj=longlat +datum=WGS84 +vunits=m +no_defs</font></div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
<br>
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 ?</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif">I would appreciate any insight anyone may have.
</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif">Thanks,</font></div>
<div><font face="Calibri,Arial,Helvetica,sans-serif">Alex</font></div>
</span></font></div>
</div>
</div>
</div>
</div>
</body>
</html>