# 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. 
#
# Les données fournies ont été obtenues à partir d'un
# plus gros fichier, extraction du site CMIP6,
# les tasmax régionales ont été calculées par moyennisation
# de tasmax par point de grille, proportionnellement à la
# surface de la région couverte par chaque point de grille. 
#
# Pour la séance de TD, il faut attribuer un jeu de données
# par groupe d'étudiant·e·s (binômes par exemple). 
# Dans le fichier csv, le séparateur est le point-virgule,
# et la décimale est la virgule.

library(ggplot2)
# commande setwd() adéquate à insérer ici
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)) 
                # données factuelles, de 2000 à 2030

# Analyse pour une région donnée

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)
thetaX <- fit_gev(donX,getParams = TRUE)
names(thetaX)
thetaX$mle
thetaX$se_mle
thetaY <- fit_gev(donY,getParams = TRUE)
thetaY$mle
thetaY$se_mle
muX <- thetaX$mle[1]
sigX <- thetaX$mle[2]
gamX <- thetaX$mle[3]
muY <- thetaX$mle[1]
sigY <- thetaY$mle[2]
gamY <- thetaY$mle[3]

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

ggplot(don,aes(x=tasmax,col=monde)) +
     geom_histogram(fill="white",alpha=0.5,position="identity") +
     geom_function(fun=function(x){10*dgev(x,loc=muX,scale=sigX,shape=gamX)},col="blue") +
     geom_function(fun=function(x){10*dgev(x,loc=muY,scale=sigY,shape=gamY)},col="red")