# TD Fondamentaux physiques des Echanges biosphere Atmoshere
# author: "| Raia Silvia Massad, Sebastien St Jean, Benjamin Loubet, Pauline Buysse\n|
#   UMR AgroParisTech-INRA EcoSys \n|\n| 20 janvier 2022\n"
# subtitle: M1 AETPF et M1 bio gÃ©ochimie appliquÃ©e


# Attendus du TD ===============================================================
# - Savoir realiser une analyse critique de donnees experimentales de Bioclimatologie.
# - Savoir evaluer un bilan radiatif et d'energie, ses differents termes et en comprendre leurs
# dynamiques.
# - evaluer les connaissances acquises pendant le module.
# - Savoir utiliser quelques fonctions de R 
# 
# ![Mat meteorologique a Grignon (48.844224 N, 1.951929 E )](ICOS-2015.jpg)
# 
# <https://ecosys.versailles-grignon.inra.fr/icos/index-fr.html>


# Le langage R et l'environnent Rstudio ========================================
# R (<https://www.r-project.org/>) est un environnement intégré de manipulation de données, de calcul et de preparation de graphiques. Toutefois, ce n'est pas seulement un "autre" environnement statistique (comme SPSS ou SAS, par exemple), mais aussi un langage de programmation complet et autonome.
# 
# Comme tous ces langages, le R est interprété, c'est-a-dire qu'il requiert un autre programme - interpréteur - pour que ses commandes soient exécutées. Par opposition, les programmes de langages compiles, comme le C ou le C++, sont d'abord convertis en code machine par le compilateur puis directement executes par l'ordinateur. Cela signifie donc que lorsque l'on programme en R, il n'est pas possible de plaider l'attente de la fin de la phase de compilation pour perdre son temps au travail.
# 
# Parmi ses caractéristiques particulièrement intéressantes, on note :
# 
#   - langage base sur la notion de vecteur, ce qui simplifie les calculs mathématiques et réduit considérablement le recours aux structures itératives (boucles for, while, etc.) ;
#   - pas de typage ni de déclaration obligatoire des variables ;
#   - programmes courts, en général quelques lignes de code seulement ;
#   - temps de développement très court.
# 
# Plus d'informations
# 
# <https://cran.r-project.org/doc/contrib/Goulet_introduction_programmation_R.pdf>
# 
# RStudio (<https://www.rstudio.com/>) est un environnement de developpement multiplateforme gratuit, libre et open source pour R.
# 
# 
# 
# On a 4 fenetres :
# 
#   - Editeur et visualisation
#   - Console
#   - Environnement et Historique
#   - Fichier, Graphs, Help


## Scripts et librairies  ======================================================
# 
# Vous pouvez sauvegarder vos script dans un fichier text qui se termine par ".R".
# Vous pourrez ainsi re-executer les scripts en R `source("myscript.R")`.
# 
### Installer et charger des "packages" ========================================

#install.packages(chron)
#install.packages(openair)
#install.packages(solaR)
#install.packages(lattice)
#install.packages(ggplot2)
#install.packages(dplyr)

install.packages("RAtmosphere_1.1.tar.gz", repos = NULL, type = "source")


library(RAtmosphere)
library(chron) # to handle time 
library(openair) # for Wind rose and corPlot
library(solaR) # for projection of irradiance

# Not necessary: but could be useful for extending graphical possibilities  
library(lattice) # to plot with lattice package
library(ggplot2) # to plot with ggplot package
library(dplyr) # to aggregate 


### Préciser le dossier de travail ==============================================
# Données d'entrées, scripts et sorties seront charges ou sauvegardes dans ce repertoire
# 
# De manière général, il faut éviter les espaces et les caractères accentues dans les nom et les chemins de vos fichiers_

#setwd("C:/Users/.../TD ICOS")


### Charger les données du TD ==================================================
#
#Données du TD sont sur la clé USB
#
# * Data
#     - Données météo et flux
#         + "FR-Gri_meteo_2013.csv"
#         + "FR-Gri_flux_2013.csv"
#     - Description des variables 
#         + "NomDesVariables.pdf"
#     - Information sur les pratiques agricoles
#         + "Itiniraire technique.xlsx"
# 
# Charger les donnees et reconnaitre le type de separateur et le format des donnees (field separator, decimal points, "NA", ...) ainsi que les dates et les heures (utiliser la fonction ci-dessous). 

#lecture des données
meteo = read.table("FR-Gri_meteo_2013.csv",header=TRUE,sep=";",dec=".",na.strings="NA")

# changement du format de date (complexe mais nécessaire)
meteo$Date.Time = as.POSIXct(strptime(meteo$Date.Time, format="%d/%m/%Y %H:%M"))
# extraction de l'heure
meteo$Time = strftime(meteo$Date.Time, format="%H:%M")
meteo$hournb = hour(meteo$Date.Time)
# date de début
startDate = meteo$Date.Time[1]
# date de fin
endDate = meteo$Date.Time[nrow(meteo)]
#jour julien 
meteo$DOY = as.integer(julian(meteo$Date.Time,startDate)) +1
#mois
meteo$month = months(meteo$Date.Time)
#numéro du mois
meteo$monthnb = as.numeric(format(meteo$Date.Time,"%m"))


#On utilisera la fonction `LoadICOS` dans le TD


### Explorer la structure des donnees ==========================================
# 
# meteo est un objet de type `data.frame`. Vous pouvez utiliser "help" et executer l'exemple du help.

# help of a known function
?data.frame
# help on an unknown name function 
??data
# run the example of the help
example(data.frame)


### Obtenir des informations sur la structure, le nom des colonnes etc =========

str(meteo) ## structure de l'object
summary(meteo) ## informations statistiques sur les donnees
names(meteo) ## nom des colonnes
head(meteo) ## premiere ligne du tableau

## Manipuler un `data.frame`====================================================

### Operations simples =========================================================

# Acceder aux donnees et aux colonnes 

# accès par index
meteo[1,]
meteo[,2]

meteo[1:3,]
meteo[c(1,3),]

#accès conditionnel
meteo[meteo$Ta>0, ]

# accès par nom de colonne
meteo$Ta

meteo[,c("Ta","Date.Time")]

#affectation d'une nouvelle colonne
meteo$NewCol = meteo$Ta + 1 


### Construire vos propres fonctions ===========================================
# 
# Exemple : fonction pour charger les données du TD
#     
    
loadICOS = function(filename)
{
 data = read.table(filename,header=TRUE,sep=";",dec=".",na.strings="NA")
 data$Date.Time = as.POSIXct(strptime(data$Date.Time, format="%d/%m/%Y %H:%M"))
 data$Time = strftime(data$Date.Time, format="%H:%M")
 data$hournb = hour(data$Date.Time)
 startDate = data$Date.Time[1]
 endDate = data$Date.Time[nrow(data)]
 data$DOY = as.integer(julian(data$Date.Time,startDate)) +1
 data$month = months(data$Date.Time)
 data$monthnb = as.numeric(format(data$Date.Time,"%m"))
 # valeur de retour
 return (data)
}

# voici comment appeler cette fonction
meteo = loadICOS("FR-Gri_meteo_2013.csv")


### Sélectionner, organiser et transformer vos donnees =========================
# 
# Exemple utiliser des accolades, `subset` et `aggregate`.
# 
#extraction des données du mois de janvier pour lesquels la température de l'air est > 0
meteo[meteo$monthnb == 1 & meteo$Ta > 0, ]

# sous-ensembels de données avec la fonction subset
meteo.janvier = subset(meteo, monthnb==1, select=c(Ta,RH)) # on garde 2 colonnes Ta, RH 
meteo.hiver = subset(meteo, monthnb %in% c(12,1,2))

# agrégation de données: aggregate
# calcul de valeurs moyennes journalières
meteo.dailymean = aggregate(meteo,by=list(meteo$DOY),FUN=mean)

# Pour aller plus loin :  
# Les package `dplyr`,`tidyr`  et l'operateur `%>%`


## Faire des graphiques ===========================================================
#
### Fonction basique "plot" dans R ----
#
plot(Ta~Date.Time,meteo,type="l")

### une fonction plus haut niveau : Lattice package xyplot ----
xyplot(SWC_5cm+SWC_20cm+SWC_30cm~Date.Time,meteo,auto.key=TRUE,type="l")

### Une fonction très haut niveau : ggplot2 package  ----
# la fonction qplot
qplot(x = Date.Time,y = Ta,data=meteo,geom=c("point","smooth"))

# la syntaxe plus générale ggplot
ggplot(data=meteo,aes(x = Date.Time,y = Ta)) + # on dit quelle variable est concernée
  geom_point() + #on précise le type de graphique ensuite 
  geom_smooth() #on peut ajouter plusiurs types de graphiques ou autre fonctions avec un "+"

### Graphique de rose des vents ----

windRose(meteo, ws="WS",wd="WD")


## Quelques fonctions ==========================================================

#manipulation de données 
cbind(meteo$Ta,meteo$SWin)

#statistiques simples 
summary(meteo$Ta)

#moyenne
mean(meteo$Ta)

#ecart-type
sd(meteo$Ta)

# Cumuls 
cumsum(c(1:5))

# agrégations de données 
aggregate(data=meteo,Ta~hournb,FUN = mean)

aggregate(data=meteo,Ta~monthnb,FUN = sd)

# Corrélations 
x = subset(meteo,select=c(Ta,SWin,RH,WS))
corPlot(x)

# L'alternative si la librairie openair ne peut être utilisé
library(corrplot)
C=cor(x,use="na.or.complete")
corrplot(C,method='square',type='upper')

# Decomposition en série temporelle
x = subset(meteo,select=SWin)[meteo$month=="juillet",]
Serie_temp=ts(x,frequenc=24)
dec=decompose(Serie_temp)
plot(dec)

# azimuth du soleil
lat=48.84424
lon=-1.951929
SZA(meteo$Date.Time[c(1,100,1000)],lat,lon)

# Régression linéaire avec la fonction lm()
x = c(1:100)
y = 2*x + rnorm(100,0,0.1*x)
df = data.frame(x,y)
lmres = lm(y~x,df)
summary(lmres)

plot(y~x)
abline(lmres,col="red")

plot(lmres)


