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
_______________________________________________ gdal-dev mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
