# Fichier permettant d'analyser les données haute-résolution
# de tasmax pour les régions de France métropolitaine.
# Le modèle est CERFACS-CNRM-CM6-1-HR, càd une version 
# haute-résolution (càd presque du modèle régional) du modèle 
# océan-atmosphère de Météo-France, avec l'aide du CERFACS.
# 
# Le fichier "yearlymax_tasmax.csv" contient simplement 
# les tasmax (max annuels) par région (les colonnes) 
# pour chaque année entre 1850 et 2100. 
#
# Ces tasmax régionales sont obtenues par moyennisation
# de tasmax par point de grille, proportionnellement à la
# surface de la région couverte par chaque point de grille 
# (voir fichier Prepare_CMIP6_France.R).
#
# 

#setwd("/Users/worms/Documents/Enseignement/TEDS/Données/Données CNRM-CM6-1-HR régions françaises")
library(ggplot2)
data=read.csv("yearlymax_tasmax.csv",sep=";",dec=",")
dataY <- subset(data,(year<=1879)) 
                # proxy des données contre-factuelles, de 1850 à 1880
dataX <- subset(data,(year>=2000)&(year<=2029)) 
# mX <- as.vector(apply(dataX[,-1],2,mean)) # temp moyennes par région sur la période 2000-2030
# mY <- as.vector(apply(dataY[,-1],2,mean)) # temp moyennes par région sur la période 1850-1880
# plot(mX,mY); abline(a=0,b=1,color="grey") 

# Analyse pour une région donnée
if(T){
numreg=1 # 1 = Île de France (2ème colonne de dataX), 
donX <- dataX[,numreg+1] # maintenant dataX est une colonne, pour l'IdF seulement
donY <- dataY[,numreg+1] 
n <- nrow(dataX)
don <- data.frame(c(donX,donY),rep(c("F","CF"),c(n,n)))
names(don) <- c("tasmax","monde")
ggplot(don,aes(x=tasmax,col=monde)) +
  geom_histogram(fill="white",alpha=0.5,position="identity")
#rug(donX,col="blue"); rug(donY,col="red")

summary(dataX)
estimX <- fit_gev(donX,getParams = TRUE)
names(estimX)
muX=estimX$mle[1]
sigX=estimX$mle[2]
gamX=estimX$mle[3]
seX=estimX$se_mle

summary(dataY)
estimY <- fit_gev(donY,getParams = TRUE)
muY=estimY$mle[1]
sigY=estimY$mle[2]
gamY=estimY$mle[3]
seY=estimY$se_mle


t=36.5
resRR <- calc_riskRatio_gev(returnValue=t,y1=donX,y2=donY)
resRR

x=seq(32,39,by=0.1)
plot(x,dgev(x,loc=muY,scale=sigY,shape=gamY),col="red",type="l")
lines(x,dgev(x,loc=muX,scale=sigX,shape=gamX),col="blue",type="l")
}