[gdal-dev] showing vertical datum in wgs84 or its variants

[Michael Smith]

All,

If I have a DEM with a vertical datum that is not wgs84, it shows properly in gdal. For example, 

COMPD_CS["WGS 84 / UTM zone 51N + EGM2008 geoid height",
    PROJCS["WGS 84 / UTM zone 51N",
        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"]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",0],
        PARAMETER["central_meridian",123],
        PARAMETER["scale_factor",0.9996],
        PARAMETER["false_easting",500000],
        PARAMETER["false_northing",0],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Easting",EAST],
        AXIS["Northing",NORTH],
        AUTHORITY["EPSG","32651"]],
    VERT_CS["EGM2008 geoid height",
        VERT_DATUM["EGM2008 geoid",2005,
            EXTENSION["PROJ4_GRIDS","egm08_25.gtx"],
            AUTHORITY["EPSG","1027"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Up",UP],
        AUTHORITY["EPSG","3855"]]]

If I then project this to WGS84 vertical, it does display as a vrt.

gdalwarp -of vrt -t_srs epsg:32651+7662 myfile.tif myfile_wgs.vrt

Coordinate System is:
COMPD_CS["WGS 84 / UTM zone 51N + WGS 84 (G1674)",
    PROJCS["WGS 84 / UTM zone 51N",
        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"]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",0],
        PARAMETER["central_meridian",123],
        PARAMETER["scale_factor",0.9996],
        PARAMETER["false_easting",500000],
        PARAMETER["false_northing",0],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Easting",EAST],
        AXIS["Northing",NORTH],
        AUTHORITY["EPSG","32651"]],
    GEOCCS["WGS 84 (G1674)",
        DATUM["World_Geodetic_System_1984_G1674",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","1155"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Geocentric X",OTHER],
        AXIS["Geocentric Y",OTHER],
        AXIS["Geocentric Z",NORTH],
        AUTHORITY["EPSG","7662"]]]

(and the z values do get properly projected). But when I convert from vrt to geotiff, 

gdal_translate -of gtiff myfile_wgs.vrt myfile_wgs.tif

It no longer shows the compd_cs

PROJCS["WGS 84 / UTM zone 51N",
    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"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",123],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32651"]]

I do have the environmental variable set to always show a compd_cs, GTIFF_REPORT_COMPD_CS=TRUE

Is this a special wgs simplification that gdal does on the vertical datum?


-- 
Michael Smith

Remote Sensing/GIS Center
US Army Corps of Engineers