#2011, Eulalie Sauthier, Master Thesis


source("IDW_.r") #charge les fonctions utiles
source("RL.r")
source("temps_retour.r")
source("maxmonthtoyear.r")
source("plot.ch.2.r")
source("plot.ch.3.r")
source("mae.r")
source("mpe.r")
source("rmse.r")
source("biais.r")
source("plot.ch.r") 
library(MASS) #appelle la librarie

long_lat_alt_t<-read.csv("Glaciers/résumé/Long_lat_ss_toutes.csv")# charge les détails des stations connues
long_lat_alt<-cbind(long_lat_alt_t[1], long_lat_alt_t[3:5])
prec_bcp_stations<-read.csv("Glaciers/precipitations_plus_stations_.csv")#charge les précipitations des stations connues
prec_bcp_stations<-prec_bcp_stations[3:length(colnames(prec_bcp_stations))]#on enlève l'année et le mois 
nbstation<-length(colnames(prec_bcp_stations)) #nombre de stations

score_mae_u<- vector(length= nbstation)+0 #déclare les vecteurs pour le calcule d'erreur
score_mpe_u<- vector(length= nbstation)+0
score_rmse_u<- vector(length= nbstation)+0
score_biais_u<- vector(length= nbstation)+0

rownames(long_lat_alt)<- long_lat_alt[ ,1] #renomme les lignes
params_obs<-matrix(nrow= nbstation,ncol=3); #déclare la matrice des GEV paramètres
for(i in  1: nbstation) #boucle sur toutes les stations
params_obs[i, 1:3]<-temps_retour(prec_bcp_stations,i) #calcule les GEV paramètres pour les station connues
rownames(params_obs)<-long_lat_alt[ ,1]#garde les noms de lignes

for ( u in 1: nbstation){ #on enlève une station à la fois pour la validation
	if (u == 1 ) { #gère le cas où la ligne enlevée est la première
		known_data<-cbind(params_obs[2: nbstation, ],long_lat_alt[2: nbstation,2:4]) 	}else if (u == nbstation){ #gère le cas où la ligne enlevée est la dernière 
		known_data<-cbind(params_obs[1: nbstation-1, ],long_lat_alt[1: nbstation-1,2:4]) 
	}else{ #gère le cas où la ligne enlevée est au milieu 
		temp <- rbind(params_obs[1:(u-1), 1:3],params_obs[(u+1):nbstation,1:3]) 
		rownames(temp)[1]<-"Engelberg" #renomme la ligne
		temp_lla <- rbind(long_lat_alt[1:(u-1),2:4], long_lat_alt[(u+1):nbstation,2:4]) 
		known_data<-cbind(temp,temp_lla)
		}
	
interpolating_data<-long_lat_alt[u,2:4] #lie les stations avec leurs données
parameter_for_RLplot<-IDW_(known_data,interpolating_data) #calcule leurs GEV paramètres

inject<-parameter_for_RLplot[nbstation, ] #ajoute la station pour l'interpoler

if (u == 1 ) { #selon le cas de figure on remet la station à interpoler
	parameter_for_RLplot<-rbind(inject ,known_data)
}else if (u == nbstation){
	parameter_for_RLplot<-rbind(known_data, inject)
}else{
	parameter_for_RLplot <-rbind(known_data[1:(u-1),], inject,known_data[(u+1): nbstation-1,])
}

return_level<-vector(length=(length(rownames(parameter_for_RLplot)))) #déclare le vecteur de niveau retour

for(i in 1:length(rownames(parameter_for_RLplot))) 
return_level[i]<- RL(parameter_for_RLplot[i,],50) #calcul les niveaux retour pour temps de retour de 50ans

}
fake_Data_all <-cbind(long_lat_alt[2:4], unlist(return_level)) #met les données et les GEV paramètres ensemble
name<-paste("../Image/Carte_Suisse_TR_toutes_uneaune_",rownames(params_obs)[u],".pdf" ) #nom du fichier
pdf(name) #sauve en pdf
plot.ch(fake_Data_all) #plot les niveaux retours en interpolant à la station enlevée 
dev.off() #ferme l'image

##
#Calcul des scores
#####
size_<-length((params_obs[ ,1])) #déclaration des vecteurs pour les calcul des scores
score_mae<- vector(length= size_)+0
score_mpe<- vector(length= size_)+0
score_rmse<- vector(length= size_)+0
score_biais<- vector(length= size_)+0

for (i in 1: size_){
#calcule l'erreur en comparant les quantiles des stations connues et de la station inconnue avec 4 différentes méthodes
out<-maxmonthtoyear(prec_bcp_stations,i) 
quantiles_obs<-sort(out)
p_k<-ppoints(sort(out))
quantiles_estimes <-qgev(p_k, parameter_for_RLplot[i,1], parameter_for_RLplot[i,2], parameter_for_RLplot[i,3])
score_mae[i]<-mae(quantiles_obs, quantiles_estimes) #sauve un scores pour chaque station retirée
score_mpe[i]<-MPE(quantiles_obs, quantiles_estimes)
score_rmse[i]<-RMSE(quantiles_obs, quantiles_estimes)
score_biais[i]<-biais(quantiles_obs, quantiles_estimes)
#sauve un plot des quantiles connus et obervés
pdf(paste("../Image/QQplot",colnames(prec_bcp_stations)[i],".pdf"))
qqplot(quantiles_obs, quantiles_estimes,log=c("x","y"),main=paste("QQplot",rownames(params_obs)[u]))
lines(seq(1,350),seq(1,350))
dev.off()

out<-maxmonthtoyear(prec_bcp_stations,u)
quantiles_obs<-sort(out)
p_k<-ppoints(sort(out))
#remplit un vecteur pour sauver les scores
quantiles_estimes <-qgev(p_k, parameter_for_RLplot[u,1], parameter_for_RLplot[u,2], parameter_for_RLplot[u,3])
score_mae_u[u]<-mae(quantiles_obs, quantiles_estimes)
score_mpe_u[u]<-MPE(quantiles_obs, quantiles_estimes)
score_rmse_u[u]<-RMSE(quantiles_obs, quantiles_estimes)
score_biais_u[u]<-biais(quantiles_obs, quantiles_estimes)
scores<-cbind(score_mae,score_mpe,score_rmse,score_biais)#sauve en une matrice
write.matrix(file=paste("Scores_toutes",rownames(params_obs)[u]), scores) #sauve la matrice
}

##
# Génère les carte selon les 4 moyens de calcul d'erreur avec un seuil
####
name<-paste("../Image/Carte_Suisse_erreur_MAE",".pdf" )
pdf(name)
plot.ch.2(fake_Data_all,scores[,1],10,"MAE")
dev.off()

name<-paste("../Image/Carte_Suisse_erreur_MPE",".pdf" )
pdf(name)
plot.ch.2(fake_Data_all,scores[,2],40,"MPE")
dev.off()

name<-paste("../Image/Carte_Suisse_erreur_RMSE",".pdf" )
pdf(name)
plot.ch.2(fake_Data_all,scores[,3],10,"RMSE")
dev.off()

name<-paste("../Image/Carte_Suisse_erreur_Biais",".pdf" )
pdf(name)
plot.ch.2(fake_Data_all,scores[,4],10,"Biais")
dev.off()

##
# Génère les carte selon les 4 moyens de calcul d'erreur avec gradient d'erreur
####

pdf(paste("../Image/Error_plot","MAE",".pdf"));source("plot.ch.3.r");plot.ch.3(fake_Data_all,scores[,1],"MAE");dev.off()
pdf(paste("../Image/Error_plot","MPE",".pdf"));source("plot.ch.3.r");plot.ch.3(fake_Data_all,scores[,2],"MPE");dev.off()
pdf(paste("../Image/Error_plot","RMSE",".pdf"));source("plot.ch.3.r");plot.ch.3(fake_Data_all,scores[,3],"RMSE");dev.off()
pdf(paste("../Image/Error_plot","Biais",".pdf"));source("plot.ch.3.r");plot.ch.3(fake_Data_all,scores[,4],"biais");dev.off()
