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",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,
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",500000], 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 [email protected] https://lists.osgeo.org/mailman/listinfo/gdal-dev
