Inconsistent read/write of geotiff spatial reference system

[bwallin]

Writing a geotiff with the SRS defined by EPSG:3413 and then reading gives significantly different SRS parameters.

import difflib
from osgeo import gdal, osr
gdal.UseExceptions()
osr.UseExceptions()

srs = osr.SpatialReference()
srs.ImportFromEPSG(3413)

# Use a virtual file system
fname = '/vsimem/tmp.tif'
drv = gdal.GetDriverByName('GTiff')
ds = drv.Create(fname, 1, 1)
assert ds.SetProjection(srs.ExportToWkt()) == 0

# Close and re-open
ds = None
ds = gdal.Open(fname)
srs_out = osr.SpatialReference()
assert srs_out.ImportFromWkt(ds.GetProjection()) == 0

# Show differences
print('Differences in PROJ.4')
print(''.join(difflib.ndiff(
    [srs.ExportToProj4() + '\n'],
    [srs_out.ExportToProj4() + '\n'])))

print('Differences in WKT')
print(''.join(difflib.ndiff(
    srs.ExportToPrettyWkt().splitlines(True),
    srs_out.ExportToPrettyWkt().splitlines(True))))

ds = None
gdal.Unlink(fname)

gives the diff output

Differences in PROJ.4
- +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
?                               -         ^^^
+ +proj=stere +lat_0=90 +lat_ts=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
?                                        ^

Differences in WKT
  PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",
      GEOGCS["WGS 84",
          DATUM["WGS_1984",
              SPHEROID["WGS 84",6378137,298.257223563,
                  AUTHORITY["EPSG","7030"]],
              AUTHORITY["EPSG","6326"]],
-         PRIMEM["Greenwich",0,
+         PRIMEM["Greenwich",0],
?                             +
-             AUTHORITY["EPSG","8901"]],
-         UNIT["degree",0.0174532925199433,
+         UNIT["degree",0.0174532925199433],
?                                         +
-             AUTHORITY["EPSG","9122"]],
          AUTHORITY["EPSG","4326"]],
      PROJECTION["Polar_Stereographic"],
-     PARAMETER["latitude_of_origin",70],
?                                    -
+     PARAMETER["latitude_of_origin",0],
-     PARAMETER["central_meridian",-45],
?                                  ^^^
+     PARAMETER["central_meridian",0],
?                                  ^
      PARAMETER["scale_factor",1],
      PARAMETER["false_easting",0],
      PARAMETER["false_northing",0],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]],
-     AXIS["X",EAST],
-     AXIS["Y",NORTH],
      AUTHORITY["EPSG","3413"]]

Note the different latitude of origin and central meridian.

[Even Rouault]

Duplicate of #3220 that has been fixed into libgeotiff upstream (and GDAL 2.0 internal libgeotiff)