[GIS] gdal python set projection of a raster not working

coordinate systemgdalraster

Gdal 1.9.2, ECW driver: ECW/ERDAS Compressed Wavelets (SDK 3.x)

I am having problems with the script that supposed to set projection of the ECW to EPSG:27700, i.e. the British_National_Grid.

Please could you have a look at the script and suggest what could be wrong?
I created two versions of teh script and neither of them work 🙁

I have checked and it does not seem to be the ECW file format/ECW sdk issue as the script does not updte the projection for Geotiffs either.

script 1:

import sys, os, os.path, string, array, math, shutil, glob, re, locale, math, traceback, time, datetime, csv
from osgeo import gdal, gdalconst, osr
from osgeo.gdalconst import *
import subprocess

driver = gdal.GetDriverByName("ECW")
driver = gdal.GetDriverByName("GTiff")

# set OSGB96
sr = osr.SpatialReference()
sr.SetProjection ("EPSG:27700")
sr_wkt = sr.ExportToWkt()
print sr_wkt

inDir = str(sys.argv[1])
for file in os.listdir(inDir):
    ds = gdal.Open(inDir + "\\" + file, gdal.GA_Update)
    if ds:
        print('Updating projection for ' + file)
        res = ds.SetProjection(sr_wkt)
        if res != 0:
            print('Setting projection failed ' + str(res))
        ds = None # save, close
    else:
        print('Could not open with GDAL: ' + file)

and script 2:

import sys, os, os.path, string, array, math, shutil, glob, re, locale, math, traceback, time, datetime, csv
from osgeo import gdal, gdalconst, osr
from osgeo.gdalconst import *
import subprocess

driver = gdal.GetDriverByName("ECW")
driver = gdal.GetDriverByName("GTiff")

ds_reference = gdal.Open('D:\\reference_geotif\\reference.tif', gdal.GA_ReadOnly)
prj_wkt = ds_reference.GetProjectionRef()
print prj_wkt

inDir = str(sys.argv[1])
for file in os.listdir(inDir):
    ds = gdal.Open(inDir + "\\" + file, gdal.GA_Update)
    if ds:
        print('Updating projection for ' + file)
        res = ds.SetProjection(prj_wkt)
        if res != 0:
            print('Setting projection failed ' + str(res))
        ds = None # save, close
    else:
        print('Could not open with GDAL: ' + file)

Best Answer

In the first script, you should use srs.ImportFromEPSG rather than srs.SetProjection.

i.e.

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

...

ds.SetProjection(srs.ExportToWkt())

It is hard to tell where the second script is going wrong without access to the reference image you are using.

Related Question