5

I have 2 DTM datasets, the first has a vertical datum of WGS84, while the second has an orthometric vertical datum (maybe EGM96). gdalinfo seems to only provide the horizontal datum info:

Size is 5551, 3250
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 16N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 16N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-87,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["EPSG",32616]]

The output above is the exact same for both datasets.

I want to convert the second dataset's vertical datum to WGS84 using gdalwarp, but I've tried many variations of the following commands to no avail (egm96_15.gtx is in /usr/share/proj):

gdalwarp in_dtm.tif out_wgs84_dtm.tif -s_srs EPSG:4326+5773 -t_srs EPSG:4979    
gdalwarp  in_dtm.tif out_wgs84_dtm.tif -s_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=/path/to/egm96_15.gtx" -t_srs "+proj=longlat +datum=WGS84 +no_def"

I've also tried converting the WGS84 DTM to EGM96, but no matter what I try, the output file's elevations aren't any different from the input file. Multiple people have asked how to do this before, but none of the answers to other posts are working for me.

Why won't it work?

Possible issues:

  • I'm using a DTM and not a DEM
  • Versions: GDAL=3.3.2, PROJ=7.2
  • Both datasets have the same DATUM tag, but store data in different reference frames
Improve this question
6
  • WGS84 (EPSG:4326) does not care about the vertical dimension. How would you expect a transformation to occur on the z-dimension then? You can also have a look there.
    – swiss_knight
    Commented Oct 5, 2021 at 21:16
  • Hmmm...do you know why the first gdalwarp command in my post didn't work?
    – Daniel C Jacobs
    Commented Oct 5, 2021 at 21:30
  • Welcome to GIS SE. As a new user, please take the Tour, which emphasizes the importance of asking One question per Question. Questions with bullet lists of questions are likely to be closed as needing more focus.
    – Vince
    Commented Oct 5, 2021 at 21:54
  • Proj supports compound CRS with syntax 2Dcode + verticalcode EPSG:4326+5773. I think you are doing the right thing. Perhaps your Proj version cannot handle the conversion.
    – user30184
    Commented Oct 6, 2021 at 6:48
  • @Vince the bulleted questions are all in support of "Why won't it work?"
    – Daniel C Jacobs
    Commented Oct 6, 2021 at 14:32

2 Answers 2

Reset to default
4

Your syntax works for me with Proj version 8.1.1, at least the Z coordinates are changing.

gdaltransform -s_srs EPSG:4326+5773 -t_srs EPSG:4979
Enter X Y [Z [T]] values separated by space, and press Return.
10 20 0
10 20 26.504430770874

What operation Proj is using can be checked with projinfo

projinfo -s EPSG:4326+5773 -t EPSG:4979

Operation No. 1:

unknown id, Inverse of WGS 84 to EGM96 height (1), 1 m, World.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

WKT2:2019 string:
COORDINATEOPERATION["Inverse of WGS 84 to EGM96 height (1)",
    SOURCECRS[
        COMPOUNDCRS["WGS 84 + EGM96 height",
            GEOGCRS["WGS 84",
                ENSEMBLE["World Geodetic System 1984 ensemble",
                    MEMBER["World Geodetic System 1984 (Transit)"],
                    MEMBER["World Geodetic System 1984 (G730)"],
                    MEMBER["World Geodetic System 1984 (G873)"],
                    MEMBER["World Geodetic System 1984 (G1150)"],
                    MEMBER["World Geodetic System 1984 (G1674)"],
                    MEMBER["World Geodetic System 1984 (G1762)"],
                    ELLIPSOID["WGS 84",6378137,298.257223563,
                        LENGTHUNIT["metre",1]],
                    ENSEMBLEACCURACY[2.0]],
                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]]],
            VERTCRS["EGM96 height",
                VDATUM["EGM96 geoid"],
                CS[vertical,1],
                    AXIS["gravity-related height (H)",up,
                        LENGTHUNIT["metre",1]]],
            ID["EPSG",9707]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            ENSEMBLE["World Geodetic System 1984 ensemble",
                MEMBER["World Geodetic System 1984 (Transit)"],
                MEMBER["World Geodetic System 1984 (G730)"],
                MEMBER["World Geodetic System 1984 (G873)"],
                MEMBER["World Geodetic System 1984 (G1150)"],
                MEMBER["World Geodetic System 1984 (G1674)"],
                MEMBER["World Geodetic System 1984 (G1762)"],
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]],
                ENSEMBLEACCURACY[2.0]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,3],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["ellipsoidal height (h)",up,
                    ORDER[3],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",4979]]],
    METHOD["PROJ-based operation method (approximate): +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=us_nga_egm96_15.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1"],
    OPERATIONACCURACY[1],
    USAGE[
        SCOPE["unknown"],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    REMARK["For WGS 84 to EGM96 height (1) (EPSG:10084): Replaces WGS 84 to EGM84 height (1) (CT 15781). Replaced by WGS 84 to EGM2008 height (1) and (2) (CTs 3858-59).  For reversible alternative see WGS 84 to WGS 84 + EGM96 height (1) (CT 9708). An executable using spherical harmonics is also available."]]
Improve this answer
9
  • Got missing grid error when running projinfo: Grid us_nga_egm96_15.tif needed but not found on the system. This file exists in /usr/share/proj, and projinfo works when run from this directory. Do I need to add this to my PATH or something?
    – Daniel C Jacobs
    Commented Oct 6, 2021 at 14:38
  • I have a /usr/share/proj containing lots of grids, and then a /usr/local/share/proj containing almost nothing. When I copy us_nga_egm96_15.tif to /usr/local/share/proj, the command works. How can I tell linux to look in /usr/local/share/proj?
    – Daniel C Jacobs
    Commented Oct 6, 2021 at 15:13
  • 1
    What do you have as PROJ_LIB now proj.org/usage/environmentvars.html?
    – user30184
    Commented Oct 6, 2021 at 15:23
  • It isn't set to anything. I tried export PROJ_LIB=/usr/share/proj, but then projinfo command results in ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db. Seems like proj was installed to the local share folder, but the resource files installed to the root share folder
    – Daniel C Jacobs
    Commented Oct 6, 2021 at 16:00
  • 1
    It seems that gdalwarp should be able to do that gdal.org/programs/gdalwarp.html Starting with GDAL 2.2, if the SRS has an explicit vertical datum that points to a PROJ.4 geoidgrids, and the input dataset is a single band dataset, a vertical correction will be applied to the values of the dataset. Ask help from the gdal-dev mailing list if you do not manage to make it work. See also lists.osgeo.org/pipermail/gdal-dev/2017-November/047762.html.
    – user30184
    Commented Oct 6, 2021 at 17:17
1

Use vDatum to convert between vertical datums. vDatum is FOSS4G and multi-platform. There is a command line version for Windows users. See this link.

Improve this answer
4
  • Tried to use on my file and it wouldn't process. Didn't get any errors, it was just stuck on 0% complete when trying to convert.
    – Daniel C Jacobs
    Commented Oct 13, 2021 at 19:12
  • Can you make your data available
    – GBG
    Commented Oct 13, 2021 at 21:25
  • Unfortunately, I cannot
    – Daniel C Jacobs
    Commented Oct 14, 2021 at 16:51
  • I know VDATUM does points and lists of points. Does it do raster and vector files too? ... Yes, the link says: Supported file format: text(ASCII), LiDAR(.LAS, *.LAZ) version 1.0, 1.1, 1.2 and 1.4 with Classification, ESRI ASCII Raster(.ASC), ESRI 3D shapefile, GeoTIFF (on Windows OS), and ERDAS IMAGINE (on Windows OS)
    – Dave X
    Commented Oct 14, 2021 at 20:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.