#!/usr/bin/python

import sys, math
from osgeo import ogr
from osgeo import osr

from osgeo import gdal

## file format of the input file with GCPs
#SOURCE;TARGET
#POINT(4431417.47 5650061.50);POINT(4431525.63 5650015.02)
#POINT(4431330.69 5650080.67);POINT(4431435.25 5650033.86)
#POINT(4431399.45 5650094.34);POINT(4431506.66 5650047.09)

if len(sys.argv) != 4:
	print "Usage: {0} <CSV file with SOURCE TARGET GCPs in WKT format> <EPSG number of the src coordinate system> <order of polynomial used for warping>".format(sys.argv[0]) 
	quit()

fname = sys.argv[1]
srcEPSG = int(sys.argv[2])
order = int(sys.argv[3])

transform = None

if srcEPSG != 25832:
	#31468:
	source = osr.SpatialReference()
	target = osr.SpatialReference()
	source.ImportFromEPSG(srcEPSG)
	target.ImportFromEPSG(25832)
	transform = osr.CoordinateTransformation(source, target)

opts = '-order {0} -f Memory'.format(order)

destPoints = []

with open(fname) as f:
	for line in f:
		src, dst = line.split(';', 2)
		if src == 'SOURCE':
			continue
		u = ogr.CreateGeometryFromWkt(src)
		v = ogr.CreateGeometryFromWkt(dst)
		opts += " -gcp {0:.12f} {1:.12f} {2:.12f} {3:.12f}".format(u.GetX(), u.GetY(), v.GetX(), v.GetY())
		destPoints.append(v)


srcDS = gdal.OpenEx(fname, open_options = ['GEOM_POSSIBLE_NAMES=SOURCE'])
ds = gdal.VectorTranslate( "", srcDS, options = opts)

layer = ds.GetLayer(0)
n = 0
delta = 0
for feature in layer:
	s = destPoints.pop(0)
	t = feature.GetGeometryRef()
	if transform:
		s.Transform(transform)
		t.Transform(transform)
#	print "{0} {1}".format(s, t)
	delta += (s.GetX() -  t.GetX()) ** 2 + (s.GetY() -  t.GetY()) ** 2
	n += 1

delta /= n
delta = math.sqrt(delta)
print("mean deviation is {0}".format(delta))