Tuesday, 23 August 2016

Python - update georeferencing information in multiple files

You need to install GDAL Python bindings first.

https://pythongisandstuff.wordpress.com/2011/07/07/installing-gdal-and-ogr-for-python-on-windows/

Here is a basic script as template:

 import os  
 import fnmatch  
 from datetime import date  
 from osgeo import gdal  
 from gdalconst import *  
 # logfile-path  
 logfile = "c:/temp/geotifheader_update.log"  
 #path to GeoTIFs  
 tif_filefolder = "S:/MG_DATEN/Rasterdaten"  
 #shift for all images - we just add a fixed shift to all images   
 shift_x = 2000000.950  
 shift_y = 999999.800  
 # get all TIF files in folder, including subfolders  
 def list_files(dir):  
   r = []  
   for root, dirs, files in os.walk(dir):  
     #for name in files:  
     for name in fnmatch.filter(files, "*.tif"):  
       r.append(os.path.join(root, name))  
   return r  
 #write message to log file  
 def writelog(message):  
   with open(logfile, 'a') as myfile:  
     myfile.write(message)  
 #init log file  
 def initlog():  
   with open(logfile, 'a') as myfile:  
     myfile.write(str(date.today())+"\n")  
 #update georeferencing information for each image  
 def setHeader(filepath):  
   dataset = gdal.OpenEx(filepath, GA_Update)  
   geotransform = dataset.GetGeoTransform()  
   if not geotransform is None:  
     print "**file**: ", filepath  
     print 'Origin = (', geotransform[0], ',', geotransform[3], ')'  
     print 'Pixel Size = (', geotransform[1], ',', geotransform[5], ')'  
     ulx = geotransform[0];  
     uly = geotransform[3];  
     #simple check whether image should be updated or not  
     if (ulx > 1000000 or ulx == 0 or uly > 1000000 or uly == 0):  
       m = "current insertion point out of bounds: " + str(ulx) + "," + str(uly)+"\n"  
       print m  
       writelog(m)  
     else:  
       ulx = ulx + shift_x  
       uly = uly + shift_y  
       lrx = ulx + (dataset.RasterXSize * geotransform[1])  
       lry = uly + (dataset.RasterYSize * geotransform[5])  
       gt = [ulx, (lrx - ulx) / dataset.RasterXSize, 0, uly, 0, (lry - uly) / dataset.RasterYSize]  
       dataset.SetGeoTransform(gt)  
       m = "image has been updated.\n"  
       print m  
       writelog(m)  
   else:  
     m = "not a geotiff.\n "  
     print m  
     writelog(m)  
 files = list_files(tif_filefolder)  
 initlog()  
 for file in files:  
   writelog("-----------------------\n")  
   writelog(file+"\n")  
   setHeader(file)  
   writelog("-----------------------\n")  
 '''  
 for single file:  
 setHeader("S:/GIS/DTMAV/Relief/relief_dtmav30.tif")  
 '''  
 # https://gis.uster.ch/dokumentation/datenkonvertierung/gdal-tools  
 # set path=C:\Program Files\GDAL;C:\Python27  

No comments:

Post a Comment