[gdal-dev] GTiff format does not preserve WKT EXTENSION nodes?

[Etienne Tourigny]

Hi all,

I have been trying to override the MODIS WKT definition (derived from
the source HDF4 files using gdal_translate), which apparently is not
sufficient for correct transformations.
The correct srs is  '+proj=sinu +R=6371007.181 +nadgrids=@null
+wktext', whereas the one stored in GTiff files is '+proj=sinu
+lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs
'

I tried to override the projection ref with gdal_translate, but the
resulting file (in GTiff format) does not keep the EXTENSION nodes.

Examples below.

Is this a limitation of GTiff format that cannot be fixed?  Is there
another format (except VRT) that can keep EXTENSION nodes?
The workaround is to override the SRS with -s_srs='...' when a
transformation is needed, but it would be better if the file itself
could store that information.

thanks,
Etienne

$ gdal_translate  -a_srs '+proj=sinu +R=6371007.181 +nadgrids=@null
+wktext' MCD12Q1.IGBP.h11v07.2002.tif tmp2.tif

$ gdalsrsinfo tmp2.tif

PROJ.4 : '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs '

OGC WKT :
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6371007.181,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Sinusoidal"],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

which is the same(almost) as the input file:

$ gdalsrsinfo MCD12Q1.IGBP.h11v07.2002.tif

PROJ.4 : '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs '

OGC WKT :
PROJCS["unnamed",
    GEOGCS["Unknown datum based upon the custom spheroid",
        DATUM["Not_specified_based_on_custom_spheroid",
            SPHEROID["Custom spheroid",6371007.181,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Sinusoidal"],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]



However, using vrt format shows that gdal_translate does set the
output srs correctly:

$ gdal_translate -of vrt -co COMPRESS=DEFLATE -a_nodata 255 -a_srs
'+proj=sinu +R=6371007.181 +nadgrids=@null +wktext'
MCD12Q1.IGBP.h11v07.2002.tif tmp2.vrt
$ gdalsrsinfo tmp2.vrt

PROJ.4 : '+proj=sinu +R=6371007.181 +nadgrids=@null +wktext'

OGC WKT :
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6371007.181,0],
            EXTENSION["PROJ4_GRIDS","@null"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Sinusoidal"],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    EXTENSION["PROJ4","+proj=sinu +R=6371007.181 +nadgrids=@null +wktext"]]


For reference here is how I generate the GTiff from HDF4 files (similar to