3

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

1

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:

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)

The differences in the PROJ.4 output are more critical in my opinion than the WKT.

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"]]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.