1

I apply the following Spatial Reference to a file:

spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(3411)    
outraster.SetProjection(spatialRef.ExportToWkt() )

geotransform = (-3850000.0, 25000.0 ,0.0 ,5850000.0, 0.0, -25000.0)
outraster.SetGeoTransform(geotransform)

If I use "icechart.GetGeoTransform()" and "icechart.GetProjection()" in Python, the output looks correct.

>>>icechart = gdal.Open(r'C:\Users\max\Desktop\test\icechart_persistencemap.tif')
>>>icechart.GetProjection()
'PROJCS["NSIDC_Sea_Ice_Polar_Stereographic_North",GEOGCS["GCS_Unspecified datum based upon 
the Hughes 1980 ellipsoid",DATUM["D_",SPHEROID["Hughes_1980",6378273,298.279411123064]],
PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Stereographic_North_P
ole"],PARAMETER["standard_parallel_1",70],PARAMETER["central_meridian",-45],PARAMETER
["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'

>>>icechart.GetGeoTransform()
(-3850000.0, 25000.0, 0.0, 5850000.0, 0.0, -25000.0)

However, gdalinfo tells me that there is no Spatial Reference defined, neither does QuantumGIS.

When I get the Projection and Geotransform from another file, rather than defining it manually, gdalinfo finds the Spatial Reference:

outraster.SetGeoTransform(icechart.GetGeoTransform())
outraster.SetProjection(icechart.GetProjection())

Also when I reproject the raster above with gdalwarp, gdalinfo and QGIS recognize the Spatial Reference:

os.system('gdalwarp -s_srs EPSG:3411 -tr 25000 -25000 -t_srs EPSG:3575 -of GTiff ' 
+ infile + ' ' + outfile) 

So what is wrong with the first lines of code in this post? Any advice?

0

It seems I found an answer, even though I do not know why.

This does NOT work:

spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(3411)
outraster.SetProjection(spatialRef.ExportToWkt() )
outband.FlushCache()

However THIS does work:

spatialRef = osr.SpatialReference()
spatialRef.ImportFromProj4('+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs')
outraster.SetProjection(spatialRef.ExportToWkt() )
outband.FlushCache()

Could it be that gdal does not like the Wkt representation of EPSG:3411 but is happy with the one created from the proj4 version? ( Proj4 and WKT are taken from here or from gdal library.)

| improve this answer | |

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.