Opened 9 years ago
Closed 9 years ago
#5889 closed defect (duplicate)
Inconsistent read/write of geotiff spatial reference system
Reported by: | bwallin | Owned by: | warmerdam |
---|---|---|---|
Priority: | normal | Milestone: | |
Component: | default | Version: | 1.9.2 |
Severity: | normal | Keywords: | geotiff srs |
Cc: |
Description
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.
Note:
See TracTickets
for help on using tickets.
Duplicate of #3220 that has been fixed into libgeotiff upstream (and GDAL 2.0 internal libgeotiff)