plotKilyosProfilesInclLastDiv.py~

from osgeo import ogr
from scipy import integrate
import matplotlib.pyplot as plt

source = ogr.Open('/home/nlibassi/Geodesy/Thesis/Project/Vector/ITRF96TM30/ProfilePoints/2001_2017_ProfilePtsWEllHts.shp')
layer = source.GetLayer()

outFile = '/home/nlibassi/Geodesy/Thesis/Project/ProfileData/ProfileScratch2.txt'

profNames = ['pr0000106', 'pr4000106', 'pr3000106', 'pr7700106', 'pr7000106', 'pr2000106', 'pr6000106', 'pr1000106', 'pr5000106']

def sameSign(x, y):
	return (x<0 and y<0) or (x>0 and y>0)

def calcArea(x1, x2, y1, y2, y3, y4):
	"""
	given 2 x-coordinates and 4 y-coordinates, 
	calculates the equation of two lines and returns 
	the difference of the two integrals.
	"""
	m1 = (y1-y2)/(x1-x2)
	m2 = (y3-y4)/(x1-x2)
	b1 = y1-m1*x1
	b2 = y3-m2*x1
	f1 = lambda x: m1*x + b1
	f2 = lambda x: m2*x + b2
	i1 = integrate.quad(f1, x1, x2)[0] #returns integral as 1st item in tuple
	i2 = integrate.quad(f2, x1, x2)[0]
	#return abs(i1-i2)
	return i1-i2


for pname in profNames:
	distList = []
	elevOrigList = []
	elevNewList = []
	#diff11mList = []
	layer.SetAttributeFilter("Filename = " + "'" + pname + "'")

	for feature in layer: 
		distList.append(feature.GetField('DistanceTo'))
		elevOrig = feature.GetField('EllHt2001')
		elevNew = feature.GetField('EllHt2017')
		elevOrigList.append(elevOrig)
		elevNewList.append(elevNew)

	plt.figure()
	plt.plot(distList, elevOrigList)
	plt.plot(distList, elevNewList)
	if pname == 'pr0000106':

		plt.plot(33.07, 36.71, 'b+')
		plt.plot(29.51, 36.67, 'g+')
		plt.axhline(y=36.71, color='b', linestyle='-')
		plt.axhline(y=36.67, color='g', linestyle='-')
		plt.fill_between(distList, elevOrigList, elevNewList)
		#will do for each profile but single test here

		hZipped = zip(elevOrigList, elevNewList)
		deltaList = []

		delta = 0
		for h in hZipped:
			if None not in h:
				delta = h[1] - h[0]
				deltaList.append(delta)
			else:
				deltaList.append(None)

		allZipped = zip(distList, elevOrigList, elevNewList, \
		deltaList)
		#print allZipped[:10]
		allZippedComp = zip(allZipped, allZipped[1:])
		#print allZippedComp[:10]
		area = 0
		areaList = []
		count = 0
		countList = []

		with open(outFile, 'a') as out:
			for a in allZippedComp:
				if None not in a[0] and None not in a[1] and sameSign(a[0][3], a[1][3]):
					localArea = calcArea(a[0][0], a[1][0], a[0][1], a[1][1], a[0][2], a[1][2])
					area += localArea
					out.write(str(localArea) + '\t' + str(area) + '\n')
					#count += 1
				elif None not in a[0] and None not in a[1] and not sameSign(a[0][3], a[1][3]):
					localArea = calcArea(a[0][0], a[1][0], a[0][1], a[1][1], a[0][2], a[1][2])
					area += localArea
					areaList.append(area)
					out.write('\n')
					#countList.append(count)
					#count = 0
					area = 0
			print areaList 		
				
	plt.suptitle('Profile ' + pname[2:5], fontsize=14)
	plt.xlabel('Distance from Benchmark (m)')
	plt.ylabel('Ellipsoidal Height (m)')
	plt.legend(['June 2001', 'Mar 2017'], loc='upper right')
	#plt.show()
	plt.savefig('/home/nlibassi/Geodesy/Thesis/Project/ProfileData/Plots/Prof' + pname[2:5] + '_200106_201703_EllHtTest2.png') 

'''
Sample a from allZippedComp:
((-43.705, 40.1947, 38.64408, -1.550619999999995), (-37.435, 39.3647, 38.61456, -0.7501400000000018))
'''