#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on September 2022

@author: Phil Vernant philippe.vernant@umontpellier.fr
"""

# import libraries
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from pyproj import Transformer

#- Centipede -------------------------------------------------------------------------------

gps = pd.read_csv('Centripede220923.csv',names=['id','date','point','geom','lat','lon','hauteur','Z','b1','b2','b3','b4','b5','b6'],index_col=False) # lecture du fichier
#gps = gps[gps['point'].str.contains('base')==False] # on elimine la ligne contenant la base
gps['sta'] =  gps.point.str[0:2].apply(str.lower) # on creer une nouvelle colonne contenant juste les 2 premières lettres
stations = gps.sta.unique() # on liste les stations
print ('points mesurées :', stations)
#conversion en UTM31
# on transforme les coordonnées WGS84 en UTM31N
transformer = Transformer.from_crs("EPSG:4326", "epsg:32631", always_xy=True)
gps['E'],gps['N'] = (transformer.transform(gps.lon,gps.lat))

#creation d'un tableau de donnees pour les valeurs moyennes  
gps_m = pd.DataFrame(columns=['sta','n_obs','E','N','Z','sE','sN','sZ'])
gps_m.sta = stations

for i, sta in enumerate(stations):   # calcul des moyennes et écarts types
#  print ('processing :',sta)
  gps_m.n_obs[i] = len(gps.E[gps.sta==sta])
  gps_m.E[i] = np.mean(gps.E[gps.sta==sta])
  gps_m.N[i] = np.mean(gps.N[gps.sta==sta])
  gps_m.Z[i] = np.mean(gps.Z[gps.sta==sta])
  gps_m.sE[i] = np.std(gps.E[gps.sta==sta])
  gps_m.sN[i] = np.std(gps.N[gps.sta==sta])
  gps_m.sZ[i] = np.std(gps.Z[gps.sta==sta])

print('Moyenne et ecrat type pour centipede:')
print(gps_m)

#- Trimble R8 ------------------------------------------------------------------------------------

gps_r8 = pd.read_csv('bat22_2009.csv',names=['point','E','N','Z'],index_col=False) # lecture du fichier
gps_r8 = gps_r8[gps_r8['point'].str.contains('base')==False] # on elimine la ligne contenant la base
gps_r8['sta'] =  gps_r8.point.str[0:2].apply(str.lower) # on creer une nouvelle colonne contenant juste les 2 premières lettres
stations_r8 = gps_r8.sta.unique() # on liste les stations
print ('points mesurées :', stations_r8)

#creation d'un tableau de donnees pour les valeurs moyennes  
gps_m_r8 = pd.DataFrame(columns=['sta','n_obs','E','N','Z','sE','sN','sZ']) # creation d'un dataframe pour mettre en memoire les resultats
gps_m_r8.sta = stations_r8 # on rempli avec le nom des stations

for i, sta in enumerate(stations_r8):   # calcul des moyennes et écarts types
#  print ('processing :',sta)
  gps_m_r8.n_obs[i] = len(gps_r8.E[gps_r8.sta==sta])
  gps_m_r8.E[i] = np.mean(gps_r8.E[gps_r8.sta==sta])
  gps_m_r8.N[i] = np.mean(gps_r8.N[gps_r8.sta==sta])
  gps_m_r8.Z[i] = np.mean(gps_r8.Z[gps_r8.sta==sta])
  gps_m_r8.sE[i] = np.std(gps_r8.E[gps_r8.sta==sta])
  gps_m_r8.sN[i] = np.std(gps_r8.N[gps_r8.sta==sta])
  gps_m_r8.sZ[i] = np.std(gps_r8.Z[gps_r8.sta==sta])

print('Moyenne et ecrat type pour le R8:')
print(gps_m_r8)


#- calcul des lignes de base pour Centipede et R8 et comparaison ---------------------------------------
lb = pd.DataFrame(columns=['sta','LB','sLB','LB_r8','sLB_r8',]) # creation d'un dataframe pour mettre en memoire les resultats
lb.sta = np.arange(0,len(stations)**2) # on rempli pour créer le dataframe, mais ça sera modifié après
k = 0
for i, station1 in enumerate(stations):
	for j in np.arange(i+1,len(stations)):
		station2 = stations[j]
		print('ligne de base :',station1, station2)
		lb.sta[k] = str(station1)+'-'+str(station2)
		DE = gps_m.loc[gps_m['sta'] == station1, 'E'].iloc[0] - gps_m.loc[gps_m['sta'] == station2, 'E'].iloc[0]
		sDE= np.sqrt((gps_m.loc[gps_m['sta'] == station1, 'sE'].iloc[0])**2 + (gps_m.loc[gps_m['sta'] == station2, 'sE'].iloc[0])**2) # propagation des incertitudes
		DN = gps_m.loc[gps_m['sta'] == station1, 'N'].iloc[0] - gps_m.loc[gps_m['sta'] == station2, 'N'].iloc[0]
		sDN= np.sqrt((gps_m.loc[gps_m['sta'] == station1, 'sN'].iloc[0])**2 + (gps_m.loc[gps_m['sta'] == station2, 'sN'].iloc[0])**2) # propagation des incertitudes
		LB = np.sqrt(DE**2+DN**2)
		sLB = np.sqrt(((DE/LB)**2 * sDE**2) + ((DN/LB)**2 * sDN**2)) # propagation des incertitudes
		lb.LB[k]  = LB
		lb.sLB[k] = sLB
		DE_r8 = gps_m_r8.loc[gps_m_r8['sta'] == station1, 'E'].iloc[0] - gps_m_r8.loc[gps_m_r8['sta'] == station2, 'E'].iloc[0]
		sDE_r8= np.sqrt((gps_m_r8.loc[gps_m_r8['sta'] == station1, 'sE'].iloc[0])**2 + (gps_m_r8.loc[gps_m_r8['sta'] == station2, 'sE'].iloc[0])**2) # propagation des incertitudes
		DN_r8 = gps_m_r8.loc[gps_m_r8['sta'] == station1, 'N'].iloc[0] - gps_m_r8.loc[gps_m_r8['sta'] == station2, 'N'].iloc[0]
		sDN_r8= np.sqrt((gps_m_r8.loc[gps_m_r8['sta'] == station1, 'sN'].iloc[0])**2 + (gps_m_r8.loc[gps_m_r8['sta'] == station2, 'sN'].iloc[0])**2) # propagation des incertitudes
		LB_r8 = np.sqrt(DE_r8**2+DN_r8**2)
		sLB_r8= np.sqrt(((DE_r8/LB_r8)**2 * sDE_r8**2) + ((DN_r8/LB_r8)**2 * sDN_r8**2)) # propagation des incertitudes
		lb.LB_r8[k]  = LB_r8
		lb.sLB_r8[k] = sLB_r8
#		print('   longueur ligne de base Centipede = ',np.round(LB,3),u"\u00B1",np.round(sLB,3),'m')
#		print('   longueur ligne de base R8        = ',np.round(LB_r8,3),u"\u00B1",np.round(sLB_r8,3),'m')
		k = k + 1

lb = lb.dropna()

print(lb)

#- calcul de la pente et prise en compte des incertitudes -------------------------------

from scipy.odr import *

x = np.array(lb['LB'], dtype=float)
inds = np.argsort(x)
x = x[inds]
y = np.array(lb['LB_r8'], dtype=float)[inds]
x_err = np.array(lb['sLB'], dtype=float)[inds]
y_err = np.array(lb['sLB_r8'], dtype=float)[inds]
print('number of samples:',len(x))

# Define a function linear through the origin
def linearo_func(P, x):
     return P[0]*x
linearo_model = Model(linearo_func) # Create a model for fitting.
datas = RealData(x, y, sx=x_err, sy=y_err) # Create a RealData object using our initiated data from above.
odr = ODR(datas, linearo_model, beta0=[1.]) # Set up ODR with the model and data.
out = odr.run() # Run the regression.
# Print the results
out.pprint()
slope = np.round(out.beta[0],0)
sd_sl = np.round(out.sd_beta[0],0)
print(slope,sd_sl)


#- 
plt.figure()
x_fit = np.linspace(0, x[-1], 10)
y_fit = linearo_func(out.beta, x_fit)
plt.errorbar(x, y, xerr=x_err, yerr=y_err, linestyle='None', marker='o', markersize=2, color='blue', zorder=10)
plt.plot(x_fit, y_fit,linestyle='-',color='grey')

plt.show()