#2011, Eulalie Sauthier, Master Thesis

#setwd("/Users/Lilosme/Documents/EPFL/pdmlilo/PDM/ACP/data") #chemin du workspace
rm(list=ls()) #nettoie les variables
graphics.off() #ferme les graphe ouverts

source("contour/plot.ch.borders.R") #charge les fonctions nécessaires
source("contour/rotate4image.R")
source("contour/addscale.R")
source("IDW_.r")
source("RL.r")
source("temps_retour.r")
source("IDW_vec.r")
require(splancs) #charge la librarie

x11(width=6.9,heigh=4.8);par(mar=c(0,0,0,0)); #fixe la taille de l'image

x<-(0:385)*1000+479500 #gridded points (1x1 km)
y<-(0:228)*1000+73500
x.mat <- matrix(rep(x,times=length(y)),ncol=length(x),byrow=T) #créer la grille de CH
y.mat <- matrix(rep(rev(y),each=length(x)),ncol=length(x),byrow=T)
z.mat <- as.matrix(read.table("contour/altitude_grid.txt"))*1000 #altitude de chaque points de la grille CH
x.vec<- as.vector(x.mat) #passe en vecteur
y.vec<- as.vector(y.mat)
z.vec<- as.vector(z.mat)
all.coord<-cbind(x.vec,y.vec,z.vec) #merge les vecteurs
all.coord<-data.frame(all.coord)
names(all.coord)<-c("x","y","z") #renomme les colonnes

###
# Calcule les paramètres GEV pour chaque points de la grille à partir de l'altitude
######

long_lat_alt_t<-read.csv("Glaciers/résumé/Long_lat_ss_toutes.csv") #charge les paramètres 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
prec_bcp_stations<-prec_bcp_stations[3:length(colnames(prec_bcp_stations))]#on supprime l'année et le mois

nbstation<-length(colnames(prec_bcp_stations)) #fixe le nombre de stations
rownames(long_lat_alt)<-long_lat_alt[ ,1] #redonne les noms aux lignes
params_obs<-matrix(nrow= nbstation,ncol=3); # créer la matrice pour les GEv paramètres des stations connues

for(i in  1: nbstation)
params_obs[i, 1:3]<-temps_retour(prec_bcp_stations,i) #génère les GEV paramètres pour chaque station connue
rownames(params_obs)<-long_lat_alt[ ,1] #rend le nom des lignes
write.matrix(file="params_obs",params_obs) #sauve les GEV paramètres
params_obs <-read.table("params_obs") #charge les GEV paramètres
known_data<-cbind(params_obs,long_lat_alt[,2:4]) #merge les GEV paramètres et les détails sur les stations
interpolating_data<-all.coord

IDW_vec(known_data,interpolating_data)# interpolation les GEV paramètres pour les 90'000 points de la grille

parameter_for_RLplot <-read.table("sauve_IDW_VEC") #ouvre les GEV paramètres
parameter_for_RLplot<-parameter_for_RLplot[2:(length(rownames(parameter_for_RLplot))),] #enlève la première ligne

return_level<-vector(length=(length(rownames(parameter_for_RLplot))))#créer vecteur pour niveaux retour
for(i in 1:length(rownames(parameter_for_RLplot)))
return_level[i]<- RL(parameter_for_RLplot[i,],50) #calcul les niveaux retour pour les 90'000 points grille
write.matrix(file="RL_toutes_50", return_level) #sauve les valeurs
return_level <-read.table("RL_toutes_50") #charge les valeurs sauvées
return_level<-as.matrix(return_level) #charge en matrix

all.coord<-cbind(x.vec,y.vec, return_level) # matrix à ploter
all.coord<-data.frame(all.coord) # passe en data.frame
names(all.coord)<-c("x","y","RL") #renomme colonnes

borders<-plot.ch.borders(unit="m")
inpol<-inpip(all.coord[,c("x","y")],borders) #On garde les points en CH

matimage<-matrix(return_level,nrow=nrow(z.mat),ncol=ncol(z.mat)) # crées matrice de niveaux retour
matimage[!1:length(return_level) %in% inpol]<-NA #remplit les valeurs hors de CH de NA

plot(range(borders[,1])+c(-1,1), range(borders[,2])+c(-40,1), type="n", axes=F)#plot contour
col<-terrain.colors(20)#palette de couleurs
image(x,y,rotate4image(matimage),col=col,lwd=1.5,axes=FALSE,xlab="",ylab="",add=T) #plot les niveaux retour
polygon(borders,border="black",lwd=1.5)
addscale(return_level[inpol],col ,unit="km",nbdigits=0,spacey=0) #ajoute l'échelle
