[gdal-dev] showing vertical datum in wgs84 or its variants
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",50], 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",50], 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",50], 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 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] showing vertical datum in wgs84 or its variants
Michael, > 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",50], > 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, EPSG:32651+7662 is an invalid CRS. The vertical part should be a vertical CRS, and not a Geocentric CRS as here. GDAL / VRT is forgiving, but there's no way to encode this in GeoTIFF. There's no way in WKT 1 to model a ProjectedCRS with a ellipsoidal height (WKT 2:2018 will allow it with the concept of Projected 3D CRS), so using plain EPSG:32651 is your best option here. ~ Actually I lied... There would be a way, but it is definitely not recommended. The historical GeoTIFF spec allows for "ellipsoid Vertical CS" ( http://geotiff.maptools.org/spec/geotiff6.html#6.3.4 ), but ISO 19111 does not allow this, and this will be likely no longer available as it in a future revision of the GeoTIFF spec. So if using listgeo to dump the TIFF you want to fix, and you insert the following in the Keyed_Information section of the listgeo dump GTCitationGeoKey (Ascii,45): "WGS 84 / UTM zone 31N + WGS84 ellipsoid height" VerticalCSTypeGeoKey (Short,1): VertCS_WGS_84_ellipsoid VerticalUnitsGeoKey (Short,1): Linear_Meter and then using geotifcp -g to forge a new GeoTIFF file, you will get COMPD_CS["WGS 84 / UTM zone 31N + WGS 84 ellipsoid 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",50], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AXIS["Easting",EAST], AXIS["Northing",NORTH], AUTHORITY["EPSG","32651"]], VERT_CS["unknown", VERT_DATUM["Not specified (based on WGS 84 ellipsoid)",2002, AUTHORITY["EPSG","6030"]], UNIT["metre",1.0, AUTHORITY["EPSG","9001"]], AXIS["Up",UP]]] But this would be considered as invalid as well by today's standards. (EPSG:6030 is not a vertical datum, but a geodetic datum) Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
