#2011, Eulalie Sauthier, Master Thesis

#to run the scipt: 
#setwd("/Users/Lilosme/Documents/EPFL/pdmlilo/PDM/ACP/data")
#source("PCA_main.r")

#clean the variables
rm(list=ls())

#library and call of others function
source("plot.band.r")
library(dyn)
library(latticeExtra)
library(RColorBrewer)
library(MASS)
library(tseries)
source("cercle_corr_saison.r")
source("cercle_corr_zones.r")
source("hiver.r")
source("sprintemps.r")
source("ete.r")
source("automne.r")
source("slope.r")
source("trend_pz.r")
source("seuil.r")
source("statio.r")
source("temps_retour.r")
source("test_statio.r")
source("monthtoyear.r")
source("maxmonthtoyear.r")
source("cross_corr.r")
source("fitexp.r")
source("IDW_.r")
source("plot.ch.r")
source("scatter_p_vs_r.r")
source("eta_chap.r")
source("RL.r")
source("mae.r")
source("MPE.r")
source("RMSE.r")
source("biais.r")
source("pot.r")
source("glac_prec_Davos.r")
source("glac_prec_engel.r")
source("glac_prec_santis.r")
source("glac_prec_sion.r")
source("forcast.r")
source("statio_mensuel.r")
source("scatter_p_vs_gl_decal.r")

####
#Import de données
####

#Valais
all_all<-read.csv("VS_METEO_DATASET/2006/data_76_10.csv")
all_vs<-all_all[,3:12]
nb_sta_vs<-10
#plateau romand
all_all_ro<-read.csv("Climap/plateau_romand.csv")
all_romand<-all_all_ro[2:6]
#zwch
all_all_zwch<-read.csv("Climap/ZWCH.csv")
all_zwch<-all_all_zwch[1:7]
toutes<-read.csv("Climap/Toutes.csv")
#engadine
all_all_enga<-read.csv("Climap/engadine.csv")
all_enga<-all_all_enga[2:4]
#tessin
all_all_tessin<-read.csv("Climap/tessin.csv")
all_tessin<-all_all_tessin[1:3]
nbstation<-10

##Glisement de terrain
GT=c(1978.250,1981.250,1986.250,1987.250,1995.333,1999.333,2000.250)
long=length(GT)
zero_gt=vector(length=long)+0
write.matrix(file="GT", GT)

#Inondations sur la Suisse
Inondations=c(1978.000,1987.167,1987.333,1987.5,1989.000,1990,1992.667,1998.000,1999.333,2000.750,2002,2002.917,2005.5,2008.917)
Inondations_prin<-c(1987.167,1987.333,1999.333)
Inondations_ete<-c(1987.5,2005.5)
Inondations_automne<-c(1992.667,2000.750)
Inondations_hiver<-c(1978.000,1989.000,1990.000,1998.000,2002,2002.917,2008.917)

long_in=length(Inondations)
zero_in=vector(length=long_in)-1

long_in_prin=length(Inondations_prin)
zero_in_prin=vector(length=long_in_prin)-1
long_in_ete=length(Inondations_ete)
zero_in_ete=vector(length=long_in_ete)
long_in_automne=length(Inondations_automne)
zero_in_automne=vector(length=long_in_automne)
long_in_hiver=length(Inondations_hiver)
zero_in_hiver=vector(length=long_in_hiver)

#NAO
NAO<-read.csv("NAO/NAO_76_10.csv")

#piézométrie
LM_import<-read.csv("PIEZO_DATASET/LM_2010.csv")
LM_2010<-LM_import[1:5]

LM_import_09<-read.csv("PIEZO_DATASET/LM_PCA_76-09.csv")
LM_2009<-LM_import_09[3:12]
nsta_pz<-length(colnames(LM_2009))

#Val d'AOSTE
VAD_70_07 <- read.csv("VAD/VAD1970_07.csv")
VAD_70_03 <- read.csv("VAD/VAD_1970_2003.csv")

ylab_<-"Cumul sur l'axe F1"

#Glaciers
prec_lp<-read.csv("Glaciers/P_long_period.csv")
prec_lp<-prec_lp[3:6]
prec_lp_ttes<-read.csv("Glaciers/P_long_period_ttes_.csv")

#température
temp_lp<-read.csv("Glaciers/P_long_period.csv")
temp_lp<-temp_lp[8:11]

#Glaciers -data
glaciers_gr<-read.csv("Glaciers/résumé/glaciers_gr.csv")
glaciers_ow<-read.csv("Glaciers/résumé/glaciers_engel.csv")
glaciers_gl<-read.csv("Glaciers/résumé/glaciers_Santis.csv")
glaciers_vs<-read.csv("Glaciers/résumé/glaciers_vs.csv")

glaciers_gr<-glaciers_gr[2:4]
glaciers_ow<-glaciers_ow[2:3]
glaciers_gl<-glaciers_gl[2:5]
glaciers_vs<-glaciers_vs[2:7]

long_lat_alt<-read.csv("Glaciers/résumé/Long_lat.csv")

############################################################################

##
#PCA
####

## Simple plot des données
#################################################

pdf("../Image/VAD_70_07.pdf")
plot(window(ts(VAD_70_07[2:9], start = c(1970, 1),end=c(2007,12) ,frequency = 12)), main="Précipitations VAD 07")
dev.off()
pdf("../Image/VAD_70_03_1.pdf")
plot(window(ts(VAD_70_03[2:8], start = c(1970, 1),end=c(2003,12) ,frequency = 12)),main="Précipitations VAD 03")
dev.off()
pdf("../Image/VAD_70_03_2.pdf")
plot(window(ts(VAD_70_03[9:12], start = c(1970, 1),end=c(2003,12) ,frequency = 12)),main="Précipitations VAD 03")
dev.off()

pdf("../Image/piezo_general_09.pdf")
plot(window(ts(LM_2009, start = c(1976, 5),end=c(2009,12) ,frequency = 12)),main="Précipitations pz 09")
dev.off()
pdf("../Image/piezo_general_10.pdf")
 plot(window(ts(LM_2010, start = c(1976, 5),end=c(2009,12) ,frequency = 12)),main="Précipitations pz 10")
dev.off()
pdf("../Image/prec_VS.pdf")
 plot(window(ts(all_vs, start = c(1976, 1),end=c(2010,12) ,frequency = 12)),main="Précipitations VS")
dev.off()
pdf("../Image/prec_ZWCH.pdf")
 plot(window(ts(all_zwch, start = c(1976, 1),end=c(2010,12) ,frequency = 12)),main="Précipitations ZWCH")
dev.off()
pdf("../Image/prec_tessin.pdf")
 plot(window(ts(all_tessin, start = c(1976, 1),end=c(2010,12) ,frequency = 12)),main="Précipitations TI")
dev.off()
pdf("../Image/prec_enga.pdf")
 plot(window(ts(all_enga, start = c(1976, 1),end=c(2010,12) ,frequency = 12)),main="Précipitations GR")
dev.off()
pdf("../Image/prec_romand.pdf")
 plot(window(ts(all_romand, start = c(1976, 1),end=c(2010,12) ,frequency = 12)),main="Précipitations Romandie")
dev.off()

##Test de stationarité des données
##################################################
for (i in 1:length(colnames(all_vs))) statio_mensuel(all_vs,i,1976,1,2010,12,12);
for (i in 1:length(colnames(all_zwch))) statio_mensuel(all_zwch,i,1976,1,2010,12,12);
for (i in 1:length(colnames(all_tessin))) statio_mensuel(all_tessin,i,1976,1,2010,12,12);
for (i in 1:length(colnames(all_enga))) statio_mensuel(all_enga,i,1976,1,2010,12,12);
for (i in 1:length(colnames(all_romand))) statio_mensuel(all_romand,i,1976,1,2010,12,12);

for (i in 1:length(colnames(VAD_70_07[2:9]))) statio_mensuel(VAD_70_07[2:9],i,1970,1,2007,12,12);
for (i in 1:length(colnames(VAD_70_03[2:13]))) statio_mensuel(VAD_70_03[2:13],i,1970,1,2003,12,12);

for (i in 1:length(colnames(LM_2010))) statio_mensuel(LM_2010,i,1976,1,2010,12,12);
for (i in 1:length(colnames(LM_2010))) statio_mensuel(LM_2009,i,1976,1,2009,12,12);


## FUNCTION calls:génération des cercle de corrélation
##################################################
z<-cercle_corr_zones(all_vs,nb_sta_vs,"Valais",1976,1,2010,12,12)
z_zwch<-cercle_corr_zones(all_zwch,7,"Suisse Centrale Ouest",1976,1,2010,12,12)
z_tessin<-cercle_corr_zones(all_tessin,3,"Tessin",1976,1,2010,12,12)
z_enga<-cercle_corr_zones(all_enga,3,"Engadine",1976,1,2010,12,12)
z_romand<-cercle_corr_zones(all_romand,5,"Plateau Romand",1976,1,2010,12,12)
z_toutes<-cercle_corr_zones(toutes,28,"Toutes les stations",1976,1,2010,12,12)
#VAD
z_vad_7007<-cercle_corr_zones(VAD_70_07[2:9],8,"Val D'Aoste 70-07",1970,1,2007,12,12)
z_vad_7003<-cercle_corr_zones(VAD_70_03[2:13],12,"Val D'Aoste 70_03",1970,1,2003,12,12)


## Génération des trends de ACP
#################################################################
print("graphe des tendances")
pdf("../Image/tendances_annuelles.pdf")
layout(matrix(c(1,2), nrow = 2), widths = c(0.9, 0.6))
par(mar = c(4, 4, 3, 2) + 0.1)
plot(z, plot.type="single", lty=1:2, main="Graphe des tendances annuelles sur la Suisse", col="red",ylab=ylab_, ylim=c(-25,60))
lines(z_zwch,col="green")
lines(z_enga,col="cyan")
lines(z_romand,col="blue")
lines(z_tessin,col="orange")
lines(z_vad_7003,col="magenta")
abline(h=0,col="black")
points(x=Inondations,y=zero_in, pch=21, col="red", bg="red")
points(x=GT,y=zero_gt, pch=23, col="black", bg="black")

par(mar = c(1, 0, 1, 5) + 0.1)
 plot(1:3, rnorm(3), pch = 1, lty = 1, ylim=c(-2,2), type = "n", axes = FALSE, ann = FALSE)
 legend("topright",c(paste("Valais", "(",round(slope(z),3)*100,"%)"),paste("Suisse Z W", "(",round(slope(z_zwch),3)*100,"%)"),paste("Engadine", "(",round(slope(z_enga),3)*100,"%)"),paste("Plateau Ro", "(",round(slope(z_romand),3)*100,"%)"),paste("Tessin", "(",round(slope(z_tessin),3)*100,"%)"),paste("VAD 70_03", "(",round(slope(z_vad_7007),3)*100,"%)"), "GT Peillettes","Inondations"), pch=c(NA,NA,NA,NA,NA,NA,23,21),pt.bg=c("red","black"), lty=1, col = c("red", "green","cyan","blue","orange","magenta", "red","black"),ncol=2) 
dev.off()

####
# PCA saisons
##############################

#valais
Valais_printemps_test<-sprintemps(all_vs, 111, nb_sta_vs)
Valais_ete_test<-ete(all_vs, 111, nb_sta_vs)
Valais_automne_test<-automne(all_vs, 111, nb_sta_vs)
Valais_hiver_test<-hiver(all_vs, 111, nb_sta_vs)

#Valais excel  
printemps_vs<-read.csv("VS_METEO_DATASET/2006/printemps.csv")
Valais_printemps<-printemps_vs[2:11]
automne_vs<-read.csv("VS_METEO_DATASET/2006/automne.csv")
Valais_automne <-(automne_vs[2:11])
ete_vs<-read.csv("VS_METEO_DATASET/2006/ete.csv")
Valais_ete<-(ete_vs[2:11])
hiver_vs<-read.csv("VS_METEO_DATASET/2006/hiver.csv")
Valais_hiver <-(hiver_vs[2:11])

#tessin
Tessin_printemps<-sprintemps(all_tessin, 111,3)
Tessin_ete<-ete(all_tessin, 111,3)
Tessin_automne<-automne(all_tessin, 111,3)
Tessin_hiver<-hiver(all_tessin, 111,3)

#Plateau romand
PR_printemps<-sprintemps(all_romand, 111,5)
PR_ete<-ete(all_romand, 111,5)
PR_automne<-automne(all_romand, 111,5)
PR_hiver<-hiver(all_romand, 111,5)

#Engadine
#Enga_printemps<-sprintemps(all_enga, 111,3)
#Enga_ete<-ete(all_enga, 111,3)
#Enga_automne<-automne(all_enga, 111,3)
#Enga_hiver<-hiver(all_enga, 111,3)

# Engadine 
printemps_enga<-read.csv("Climap/print_enga.csv")
Enga_printemps <-printemps_enga[1:3]
automne_enga<-read.csv("Climap/aut_enga.csv")
Enga_automne <-automne_enga[1:3]
ete_enga<-read.csv("Climap/ete_enga.csv")
Enga_ete <-ete_enga[1:3]
hiver_enga<-read.csv("Climap/hiv_enga.csv")
Enga_hiver <-hiver_enga[1:3]

#ZWCH
ZWCH_printemps<-sprintemps(all_zwch, 111,7)
ZWCH_ete<-ete(all_zwch, 111,7)
ZWCH_automne<-automne(all_zwch, 111,7)
ZWCH_hiver<-hiver(all_zwch, 111,7)

#VAD
nbst<-8 #Nb de stations
data_vad <-VAD_70_07[2:9]

VAD_printemps<-sprintemps(data_vad, 140, nbst)
VAD_ete<-ete(data_vad, 140, nbst)
VAD_automne<-automne(data_vad, 140, nbst)
VAD_hiver<-hiver(data_vad, 140, nbst)
filename_vad<-"VAD_70_07"

###SPRING###############################################################

zed_valais_spring<-cercle_corr_zones(Valais_printemps, nb_sta_vs,"Valais_printemps",1976,3,2010,5,3)
zed_zwch_spring<-cercle_corr_zones(ZWCH_printemps,7,"ZWCH_printemps",1976,3,2010,5,3)
zed_tessin_spring<-cercle_corr_zones(Tessin_printemps,3,"Tessin_printemps",1976,3,2010,5,3)
zed_enga_spring<-cercle_corr_zones(Enga_printemps,3,"Engadine_printemps",1976,3,2010,5,3)
zed_romand_spring<-cercle_corr_zones(PR_printemps,5,"Plateau romand_printemps",1976,3,2010,5,3)
vad_prin<-cercle_corr_zones(VAD_printemps, nbst,paste(filename_vad,"Printemps"),1976,3,2010,5,3)


pdf("../Image/tendances_spring.pdf")
layout(matrix(c(1,2), nrow = 2), widths = c(0.9, 0.6))
par(mar = c(4, 4, 4, 2) + 0.1)
plot(zed_valais_spring, main="Printemps",col="red",ylab="Cumul sur l'axe F1",xlab="Temps",ylim=c(-20,25))

lines(zed_zwch_spring,col="green")
lines(zed_enga_spring,col="cyan")
lines(zed_romand_spring,col="blue")
lines(zed_tessin_spring,col="orange")
lines(vad_prin,col="magenta")
abline(h=0,col="black")
points(x=GT,y=zero_gt, pch=23, col="red", bg="red")

par(mar = c(1, 2, 1, 5) + 0.1)
 plot(1:3, rnorm(3), pch = 1, lty = 1, ylim=c(-2,2), type = "n", axes = FALSE, ann = FALSE)
 legend("topright",c(paste("Valais", "(",round(slope(zed_valais_spring),3)*100,"%)"),paste("Suisse Z-W", "(",round(slope(zed_zwch_spring),3)*100,"%)"),paste("Engadine", "(",round(slope(zed_enga_spring),3)*100,"%)"),paste("Plateau romand", "(",round(slope(zed_romand_spring),3)*100,"%)"),paste("Tessin", "(",round(slope(zed_tessin_spring),3)*100,"%)"), paste("VAD", "(",round(slope(vad_prin),3)*100,"%)"),"Glissements"), pch=c(NA,NA,NA,NA,NA,NA,23),pt.bg=c("red"), lty=1, col = c("red", "green","cyan","blue","orange","magenta","red"),ncol=2)
dev.off()
#####ETE###############################################################

print("ete")
zed_ete<-cercle_corr_zones(Valais_ete, nb_sta_vs,"Valais_ete",1976,6,2010,8,3)
zed_zwch_ete<-cercle_corr_zones(ZWCH_ete,7,"ZWCH_ete",1976,6,2010,8,3)
zed_tessin_ete<-cercle_corr_zones(Tessin_ete,3,"Tessin_ete",1976,6,2010,8,3)
zed_enga_ete<-cercle_corr_zones(Enga_ete,3,"Engadine_ete",1976,6,2010,8,3)
zed_romand_ete<-cercle_corr_zones(PR_ete,5,"Plateau_romand_ete",1976,6,2010,8,3)
vad_summer<-cercle_corr_zones(VAD_ete, nbst,paste(filename_vad,"Eté"),1976,6,2010,8,3)

pdf("../Image/tendances_été.pdf")
layout(matrix(c(1,2), nrow = 2), widths = c(0.9, 0.6))
par(mar = c(4, 5, 3, 1) + 0.1)
plot(zed_ete, main="Eté",col="red",ylab="Cumul sur l'axe F1",ylim=c(-30,25))
lines(zed_zwch_ete,col="green")
lines(zed_enga_ete,col="cyan")
lines(zed_romand_ete,col="blue")
lines(zed_tessin_ete,col="orange")
lines(vad_summer,col="magenta")
abline(h=0,col="black")
points(x=Inondations_ete,y=zero_in_ete, pch=21, col="black",bg="black")

par(mar = c(7, 7, 3, 0) + 0.1)
 plot(1:3, rnorm(3), pch = 1, lty = 1, ylim=c(-2,2), type = "n", axes = FALSE, ann = FALSE)
legend("bottomleft",c(paste("Valais", "(",round(slope(zed_ete),3)*100,"%)"),paste("Suisse Z-W", "(",round(slope(zed_zwch_ete),3)*100,"%)"),paste("Engadine", "(",round(slope(zed_enga_ete),3)*100,"%)"),paste("Plateau romand", "(",round(slope(zed_romand_ete),3)*100,"%)"),paste("Tessin", "(",round(slope(zed_tessin_ete),3)*100,"%)"),paste("VAD", "(",round(slope(vad_summer),3)*100,"%)"),"Inondations"), pch=c(NA,NA,NA,NA,NA,NA,21),pt.bg=c("black"), lty=1, col = c("red", "green","cyan","blue","orange","magenta","black"),ncol=2)
dev.off()

######AUTOMNE#################################################################

zed_valais_aut<-cercle_corr_zones(Valais_automne, nb_sta_vs,"Valais_aut",1976,9,2010,11,3)
zed_zwch_aut<-cercle_corr_zones(ZWCH_automne,7,"ZWCH_aut",1976,9,2010,11,3)
zed_tessin_aut<-cercle_corr_zones(Tessin_automne,3,"Tessin_aut",1976,9,2010,11,3)
zed_enga_aut<-cercle_corr_zones(Enga_automne,3,"Engadine_aut",1976,9,2010,11,3)
zed_romand_aut<-cercle_corr_zones(PR_automne,5,"Plateau romand_aut",1976,9,2010,11,3)
vad_fall<-cercle_corr_zones(VAD_automne, nbst,paste(filename_vad,"Automne"),1976,9,2010,11,3)
pdf("../Image/tendances_automne.pdf")
layout(matrix(c(1,2), nrow = 2), widths = c(0.9, 0.6))
par(mar = c(4, 3, 3, 2) + 0.1)
plot(zed_valais_aut, main="Automne",col="red",ylab="Cumul sur l'axe F1",ylim=c(-20,30))
lines(zed_zwch_aut,col="green")
lines(zed_enga_aut,col="cyan")
lines(zed_romand_aut,col="blue")
lines(zed_tessin_aut,col="orange")
lines(vad_fall,col="magenta")
points(x=Inondations_automne,y=zero_in_automne, pch=21, col="black",bg="black")
abline(h=0,col="black")
par(mar = c(1, 4, 1, 5) + 0.1)
plot(1:3, rnorm(3), pch = 1, lty = 1, ylim=c(-2,2), type = "n", axes = FALSE, ann = FALSE)
legend("topleft",c(paste("Valais", "(",round(slope(zed_valais_aut),3)*100,"%)"),paste("Suisse Z-W", "(",round(slope(zed_zwch_aut),3)*100,"%)"),paste("Engadine", "(",round(slope(zed_enga_aut),3)*100,"%)"),paste("Plateau romand", "(",round(slope(zed_romand_aut),3)*100,"%)"),paste("Tessin", "(",round(slope(zed_tessin_aut),3)*100,"%)"),paste("VAD", "(",round(slope(vad_fall),3)*100,"%)"),"Inondations"), pch=c(NA,NA,NA,NA,NA,NA,21),pt.bg=c("black"),lty=1, col = c("red", "green","cyan","blue","orange","magenta","black"),ncol=2)
dev.off()

########Winter###########################################

print("hiver")
zed_valais_h<-cercle_corr_zones(Valais_hiver, nb_sta_vs,"Valais_hiver",1976,1,2010,12,3)
write.matrix(file="z_valais_hiver", zed_valais_h)

zed_zwch_h<-cercle_corr_zones(ZWCH_hiver,7,"ZWCH_hiver",1976,1,2010,12,3)
zed_tessin_h<-cercle_corr_zones(Tessin_hiver,3,"Tessin_hiver",1976,1,2010,12,3)
zed_enga_h<-cercle_corr_zones(Enga_hiver,3,"Engadine_hiver",1976,1,2010,12,3)
zed_romand_h<-cercle_corr_zones(PR_hiver,5,"Plateau romand_hiver",1976,1,2010,12,3)
vad_winter<-cercle_corr_zones(VAD_hiver, nbst,paste(filename_vad,"Hiver"),1976,1,2010,12,3)
pdf("../Image/tendances_hiver.pdf")
layout(matrix(c(1,2), nrow = 2), widths = c(0.9, 0.6))
par(mar = c(5, 4, 3, 2) + 0.1)
plot(zed_valais_h, main="Hiver",col="red",ylab="Cumul sur l'axe F1",ylim=c(-28,35))
lines(zed_zwch_h,col="green")
lines(zed_enga_h,col="cyan")
lines(zed_romand_h,col="blue")
lines(zed_tessin_h,col="orange")
lines(vad_winter,col="magenta")
abline(h=0,col="black")
points(x=Inondations_hiver,y=zero_in_hiver, pch=21, col="black",bg="black")
par(mar = c(9, 6, 3, 3) + 0.1)
plot(1:3, rnorm(3), pch = 1, lty = 1, ylim=c(-2,2), type = "n", axes = FALSE, ann = FALSE)
legend("bottomleft",c(paste("Valais", "(",round(slope(zed_valais_h),3)*100,"%)"),paste("Suisse Z-W", "(",round(slope(zed_zwch_h),3)*100,"%)"),paste("Engadine", "(",round(slope(zed_enga_h),3)*100,"%)"),paste("Plateau romand", "(",round(slope(zed_romand_h),3)*100,"%)"),paste("Tessin", "(",round(slope(zed_tessin_h),3)*100,"%)"),paste("VAD", "(",round(slope(vad_winter),3)*100,"%)"),"Inondations"), pch=c(NA,NA,NA,NA,NA,NA,21),pt.bg=c("black"), lty=1, col = c("red", "green","cyan","blue","orange","magenta","black"),ncol=2)
dev.off()


###
#PIEZOMETRIE
#################################################################################

print("piézométrie")
region_10<-"Tendances Valaisannes et piézométriques annuelles"
LM_10<-cercle_corr_zones(LM_2010,5,"Piézomètres 2010",1976,1,2010,12,12)
pdf("../Image/tendances_lm_all_2010.pdf")

layout(matrix(c(1,2), nrow = 2), widths = c(0.9, 0.6))
par(mar = c(4, 5, 3, 3) + 0.1)
plot(LM_10, plot.type="single", lty=1:2, main=region_10, col="red",ylab="Cumul sur l'axe PC1",xlab="Temps")
lines(z, col="blue")
abline(h=0,col="black")

par(mar = c(1, 3, 2, 6) + 0.1)
plot(1:3, rnorm(3), pch = 1, lty = 1, ylim=c(-2,2), type = "n", axes = FALSE, ann = FALSE)
legend("topright",c("Pfiézomètres","Stations valaisannes"), lty=1, col = c("red","blue"),ncol=2)
dev.off()

#piezometrie 76-2009
######################

LM_09<-cercle_corr_zones(LM_2009,nsta_pz,"Piézomètres 2009",1976,1,2009,12,12)

pdf("../Image/tendances_lm_all_09.pdf")
layout(matrix(c(1,2), nrow = 2), widths = c(0.9, 0.6))
par(mar = c(1, 3, 3, 2) + 0.1)

plot(LM_09, plot.type="single", lty=1:2, main=region_10, col="red",ylab="Cumul sur l'axe F1")
lines(z, col="blue")
abline(h=0,col="black")#points(x=GT,y=zero_gt, pch=23, col="red")
par(mar = c(1, 3, 2, 8) + 0.1)
plot(1:3, rnorm(3), pch = 1, lty = 1, ylim=c(-2,2), type = "n", axes = FALSE, ann = FALSE)
legend("topright",c("Piézométrie","Pluies"), lty=1, col = c("red","blue"),ncol=2)
dev.off()

#seasonal PZ cirlce and trends
##############################

#2009
#####

PZ_printemps<-sprintemps(LM_2009, 140,nsta_pz)
PZ_ete<-ete(LM_2009, 140,nsta_pz)
PZ_automne<-automne(LM_2009, 140, nsta_pz)
PZ_hiver<-hiver(LM_2009, 140, nsta_pz)

pz_1<-cercle_corr_zones(PZ_printemps, nsta_pz,"PZ 09 Printemps",1976,3,2010,5,3)
pz_2<-cercle_corr_zones(PZ_ete, nsta_pz,"PZ 09 Eté",1976,6,2010,8,3)
pz_3<-cercle_corr_zones(PZ_automne, nsta_pz,"PZ 09 Automne",1976,9,2010,11,3)
pz_4<-cercle_corr_zones(PZ_hiver, nsta_pz,"PZ 09 Hiver",1976,1,2010,12,3)
write.matrix(file="pz_valais_hiver", pz_4)

trend_pz(pz_1, nsta_pz,"printemps_09_pz",z,zero_09)
trend_pz(pz_2, nsta_pz,"ete_09_pz",z,zero_09)
trend_pz(pz_3, nsta_pz,"automne_09_pz",z,zero_09)
trend_pz(pz_4, nsta_pz,"hiver_09_pz",z,zero_09)

#2010
#####

PZ_printemps_10<-sprintemps(LM_2010, 140,5)
PZ_ete_10<-ete(LM_2010, 140,5)
PZ_automne_10<-automne(LM_2010, 140,5)
PZ_hiver_10<-hiver(LM_2010, 140,5)

pz_1_10<-cercle_corr_zones(PZ_printemps_10,5,"PZ 10 Printemps",1976,3,2010,5,3)
pz_2_10<-cercle_corr_zones(PZ_ete_10,5,"PZ 10 Eté",1976,6,2010,8,3)
pz_3_10<-cercle_corr_zones(PZ_automne_10,5,"PZ 10 Automne",1976,9,2010,11,3)
pz_4_10<-cercle_corr_zones(PZ_hiver_10,5,"PZ 10 Hiver",1976,1,2010,12,3)

trend_pz(pz_1_10,5,"printemps_10_pz",z,zero_10)
trend_pz(pz_2_10,5,"ete_10_pz",z,zero_10)
trend_pz(pz_3_10,5,"automne_10_pz",z,zero_10)
trend_pz(pz_4_10,5,"hiver_10_pz",z,zero_10)

###
#cross-correlation
#####################

cross_corr(LM_10,z) #valais
prewhiten(LM_10,z)

#####
#NAO
#################################################


#mettre sous la bonne forme pour utiliser ts()
data_m<- vector(length=420)
k<-1
for (i in 1:35){
for (j in 1:12){
	data_m[k] <- NAO[i,j] 
	k<-k+1
	}}
#transorme en série temporelle	
data_ts<-ts(data_m,start=c(1976,1),end=c(2010,12),frequency=12)

pdf("../Image/nao.pdf")
layout(matrix(c(1,2), nrow = 2), widths = c(0.9, 0.6))
par(mar = c(4, 5, 3, 3) + 0.1)
plot(smooth(data_ts),type="h", col="blue",ylab="Valeurs NAO") 
par(mar = c(4, 5, 3, 3) + 0.1)
box(which = "figure")
plot(z,col="red", ylab="Cumul sur PC1",ylim=c(-25,60))
lines(z_zwch,col="green")
lines(z_romand,col="blue")
lines(z_tessin, col="orange")
lines(z_vad_7003,col="magenta")
legend("topright",c("Valais", "ZWCH", "Romandie","Tessin", "VAD"), col=c("red","green", "blue", "orange", "magenta"),lty=1, cex=0.6)
dev.off()


##
#VAD
###################################

pdf("../Image/VAD0703_vs.pdf")
plot(z_vad_7007, plot.type="single", lty=1:2, main="Val d'Aoste", col="red",ylab="Cumul sur l'axe F1")
legend("topright",c(paste("VAD 70-07", "(",round(slope(z_vad_7007),3)*100,"%)"),paste("VAD 70-03", "(",round(slope(z_vad_7003),3)*100,"%)")), pch=21, pt.bg="white", lty=1, col = c("red", "green"))
lines(z_vad_7003,col="green")
abline(h=0,col="black")
lines(z,"black",type="l")
dev.off()

#####################################################################################################
# PART 2 : donnéees période plus longue +glaciers
#####################################################################################################

#####
#Precipitation - glaciers
#####

glaciers_gr_ts<-ts(glaciers_gr[ ,1], start=c(1893,1))
pdf("../Image/Précipitation_lp.pdf")
plot(window(ts(prec_lp, start = c(1864, 1),end=c(2011,3) ,frequency = 12)))
dev.off()

#test la stationarité des données générale
for (i in 1:length(colnames(prec_lp))){
statio_mensuel(prec_lp,i,1864,1,2011,3,12)
statio_mensuel(temp_lp,i,1864,1,2011,3,12);
}


####
# Glaciers et Prec
####
print("4st")
#a lancer individuellement
#glac_prec_Davos()
#glac_prec_engel()
#glac_prec_santis()
#glac_prec_sion()

#####
#Stationarité, temps retour et autocorrélation des Précipitations plus Températures
#####

residu<-matrix(nrow=500, ncol=4)
i<-0
for (i in 1:length(colnames(prec_lp))){
seuille_lp<-150
seuille_lp_T<-15
statio_mensuel(prec_lp,i,1864,1,2011,3,12);
statio_mensuel(temp_lp,i,1864,1,2011,3,12);
temps_retour(prec_lp,i)
temps_retour(temp_lp,i)
}

graphics.off()

#####
#Peak-over-Threshold (POT)
#####
size_fenetre<-10
fen<-matrix(nrow=47,ncol=15)

for(i in 1:4){
buf<-pot(prec_lp,i,50,1864,2011,3,12)
buf_T<-pot(temp_lp,i, 50,1864,2011,3,12)
}

#####
#Etat données - glaciers
#####
seuille_l<-(-300);

for (i in 1:length(colnames(glaciers_gr)))
statio_mensuel(glaciers_gr,i,1893,1,2009,1,1);
for (i in 1:length(colnames(glaciers_ow)))
statio_mensuel(glaciers_ow,i,1894,1,2009,1,1);

for (i in 1:length(colnames(glaciers_gl)))
statio_mensuel(glaciers_gl,i,1912,1,2009,1,1);

for (i in 1:length(colnames(glaciers_vs)))
statio_mensuel(glaciers_vs,i,1847,1,2009,1,1);
graphics.off()

####
#FORCASTING THE DATA
####
forcast()

###########################################################################################
#PART 3
############################################################################################
print("beging of GEV: shit")
#stop()
#####
#Théorie des valeurs extrêmes
#####

#cette partie est là pour tester la stationarité des stations. GEV_RL.r permet de généraler la carte de niveaux retour et GEV_une_toutes.r permet d générer la carte des erreurs sur les stations.

long_lat_alt_t<-read.csv("Glaciers/résumé/Long_lat_ss_toutes.csv")#
long_lat_alt<-cbind(long_lat_alt_t[1], long_lat_alt_t[3:5])
print("2.74")
prec_bcp_stations<-read.csv("Glaciers/precipitations_plus_stations__.csv")prec_bcp_stations<-prec_bcp_stations[3:length(colnames(prec_bcp_stations))]#on enlève l'année et le mois 

######
# on test la stationarité des milles stations pour supprimer celles qui ne sont pas stationaires
######

kpps<-matrix(nrow=length(colnames(prec_bcp_stations)),ncol=3)
for (i in 1:length(colnames(prec_bcp_stations))){
		print(i)

	print("sort de seuil")
	#test stationarité graphiquement
	statio_mensuel(prec_bcp_stations,i,1864,1,2011,3,12);
	#test stationarité numériquement
	kpps[i,1]<-colnames(prec_bcp_stations)[i]
	buff<-kpss.test(na.omit(prec_bcp_stations[,i]))
	kpps[i,2]<-buff$statistic
	kpps[i,3]<-buff$p.value
	}
