ENH : New fct an script python to add into tiff files, GCPs

python_src/utils/addGCP.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+""" 
+    addGCP.py
+    ==========
+
+    Python script to add into a tiff file, GCPs information
+ 
+"""
+
+import argparse
+import gdal
+import osr
+import re
+
+
+# addGCP function
+def addGCP(dictKWL, inTiff):
+    """ 
+    Function to add into a tiff file, GCPs information
+    """
+
+    ###### Create GCPs list from input dict ######
+    # Select keys with gcp in it
+    filtered_dictKWL = {k:v for k,v in dictKWL.items() if "gcp" in k}
+
+    # Create 5 lists : for colunms (pixel), lines, lon, lat and hgt  
+    col_list = [float(filtered_dictKWL[k])
+                for k in filtered_dictKWL.keys() if "im_pt.x" in k]
+    line_list = [float(filtered_dictKWL[k])
+                 for k in filtered_dictKWL.keys() if "im_pt.y" in k]
+    lon_list = [float(filtered_dictKWL[k])
+                for k in filtered_dictKWL.keys() if "world_pt.lon" in k]
+    lat_list = [float(filtered_dictKWL[k])
+                for k in filtered_dictKWL.keys() if "world_pt.lat" in k]
+    hgt_list = [float(filtered_dictKWL[k])
+                for k in filtered_dictKWL.keys() if "world_pt.hgt" in k]
+    
+    gcp_number = int(filtered_dictKWL['support_data.geom.gcp.number'])
+
+    # Check list size
+    if len(col_list) != gcp_number or len(line_list) != gcp_number or \
+       len(lon_list) != gcp_number or len(lat_list) != gcp_number or \
+       len(hgt_list) != gcp_number :
+        print("Wrong size for gcp lists ")
+        quit()
+
+    gcp_list = []
+    for i in range(0, gcp_number):
+        gcp = gdal.GCP(lon_list[i], lat_list[i], hgt_list[i],
+                       col_list[i], line_list[i]) # lon, lat, hgt, col, line
+        gcp_list.append(gcp)
+
+    ###### Add GCPs into the input tiff ######
+    ds = gdal.Open(inTiff, gdal.GA_Update)
+
+    # First : Adapt Projetion for WGS84
+    wkt = ds.GetProjection()
+
+    # if no projection defined => projection by default (ESPG 4326)
+    if not wkt :
+        sr = osr.SpatialReference()
+        sr.ImportFromEPSG(4326) # ESPG : 4326 = WGS84
+        wkt = sr.ExportToWkt()
+
+    # Set GCPs into tiff image
+    ds.SetGCPs(gcp_list, wkt)
+        
+    ds = None
+    
+# Main
+if __name__ == "__main__":
+
+    ###### Get the main argument : geom file and input tiff ######
+    # Check arguments
+    parser = argparse.ArgumentParser()
+    parser.add_argument("geomfile", help="input geom file with correct GCPs")
+    parser.add_argument("intiff", help="input tiff image to include GCPs in it")
+    args = parser.parse_args()
+    
+    ###### Read geom file ######
+    f_geom =open(args.geomfile, "r")
+    keywordlist = f_geom.read().split("\n") 
+    keywordlist = filter(None, keywordlist)
+    # Transfrom file content to dict
+    dictKWL = { i.split(':')[0] : re.sub(r"[\n\t\s]*", "", i.split(':')[1]) for i in keywordlist }
+
+    addGCP(dictKWL, args.intiff)