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