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