I am writing a simple geotiff using the GDAL python bindings, but the projection WKT reported by gdalinfo
differs from the source. The script
from osgeo import gdal, osr
import numpy
# Create synthetic data
raster = numpy.zeros((512, 512), dtype=numpy.uint8)
raster[256, 256] = 1
top_left = (100000, 200000)
pixel_size = 30
output_fn = 'test.tif'
driver = gdal.GetDriverByName('GTiff')
geotiff = driver.Create(output_fn, 512, 512, 1, gdal.GDT_Byte)
src_geotran = [top_left[0], pixel_size, 0,
top_left[1], 0, -pixel_size]
geotiff.SetGeoTransform(src_geotran)
srs = osr.SpatialReference()
srs.ImportFromEPSG(3413)
srs_wkt = srs.ExportToPrettyWkt()
print(srs_wkt)
geotiff.SetProjection(srs_wkt)
geotiff.GetRasterBand(1).WriteArray(raster)
geotiff = None
produces the 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,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin",70],
PARAMETER["central_meridian",-45],
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"]]
Whereas the WKT from gdalinfo is
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],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]],
PROJECTION["Polar_Stereographic"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","3413"]]
Notably, the latitude_of_origin and central_meridian show as 0 in gdalinfo, when they should be 70 and -45 respectively.
Tested with Python 2.7.8, GDAL 1.11.1, Numpy 1.9.1
Best Answer
Looks like a bug to me, since the projection parameters were not correctly stored or read when reopening (I can't tell which is correct). No errors were raised, and Set / Get / Import methods return a 0 (= good) status, so there are no signals of any errors.
You can add this bug to here. If you need a User ID, it is from here. It's also good practice to see if this bug is already listed, but it can be difficult to search for this kind of bug.
Here is a boiled down version of what you are seeing, tested with GDAL 1.11.1:
The differences in the PROJ.4 output are more critical in my opinion than the WKT.