[gdal-dev] Setting vertical units in a raster

Even Rouault even.rouault at spatialys.com
Thu Oct 29 13:03:20 PDT 2015


Jack,

> Hey,
> So I am exporting a raster and would like to set the VerticalUnitsGeoKey to
> feet. I have tried to use the SetUnitType( const char * ) member function
> of GDALRasterBand without success.

SetUnitType() will write a GDAL metadata in the TIFFTAG_GDAL_METADATA XML tag

> I have looked at the internal code the
> driver uses to set this flag, but would like to avoid having to reopen the
> file to do this after I have finished the export. Is there any standard
> wayto set GeoTIFF tags?

See the following Python script for an example :

from osgeo import gdal,osr

sr = osr.SpatialReference()
sr.SetFromUserInput("""COMPD_CS["unknown",
PROJCS["WGS 84 / UTM zone 11N",
    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",-117],
    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","32611"]],
VERT_CS["Unnamed",
    VERT_DATUM["North American Vertical Datum 1988",2005,
        AUTHORITY["EPSG","5103"]],
    UNIT["foot",0.3048,AUTHORITY["EPSG","9002"]],
    AXIS["Up",UP]]]""")

ds = gdal.GetDriverByName('GTiff').Create('/vsimem/tiff_srs_compd_cs.tif',1,1)
ds.SetProjection(sr.ExportToWkt())
ds = None

# Check back
gdal.SetConfigOption('GTIFF_REPORT_COMPD_CS','YES')
ds = gdal.Open('/vsimem/tiff_srs_compd_cs.tif')
wkt = ds.GetProjectionRef()
gdal.SetConfigOption('GTIFF_REPORT_COMPD_CS',None)
print(wkt)

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com


More information about the gdal-dev mailing list