Krovak projection not properly stored by GTiff

[Mike Taves]

The GTiff driver does not properly store Krovak projections, namely EPSG codes 2065, 5221, 5513 and 5514. While the EPSG codes can be imported by OSR, and exported to PROJ.4 and WKT outputs, the data is lost when written to GeoTIFF.

from osgeo import osr, gdal
osr.UseExceptions()
gdal.UseExceptions()
drv = gdal.GetDriverByName('GTiff')
for code in [2065, 5221, 5513, 5514]:
    print('EPSG:' + str(code))
    fname = '/vsimem/tmp' + str(code)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(2065)
    print('Input PROJ.4: ' + srs.ExportToProj4())
    ds = drv.Create(fname, 2, 3)
    ds.SetProjection(srs.ExportToWkt())
    ds = None
    ds = gdal.Open(fname)
    wkt_out = ds.GetProjection()
    ds = None
    gdal.Unlink(fname)
    srs_out = osr.SpatialReference()
    srs_out.ImportFromWkt(wkt_out)
    print('Output PROJ.4: ' + srs_out.ExportToProj4())

shows

EPSG:2065
Input PROJ.4: +proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +towgs84=589,76,480,0,0,0,0 +pm=ferro +units=m +no_defs 
Output PROJ.4: 
EPSG:5221
Input PROJ.4: +proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +towgs84=589,76,480,0,0,0,0 +pm=ferro +units=m +no_defs 
Output PROJ.4: 
EPSG:5513
Input PROJ.4: +proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +towgs84=589,76,480,0,0,0,0 +pm=ferro +units=m +no_defs 
Output PROJ.4: 
EPSG:5514
Input PROJ.4: +proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +towgs84=589,76,480,0,0,0,0 +pm=ferro +units=m +no_defs 
Output PROJ.4: 

Note the none of these have an output PROJ.4, just an empty string. Further note that the first line of WKT is "LOCAL_CS["S-JTSK (Ferro) / Krovak",", rather than PROJCS["S-JTSK (Ferro) / Krovak", from the input srs.