---
title: "Concentration de CO2 dans l'atmosphère depuis 1958"
author: "Pierre LACOSTE"
date: "13/12/2021"
output: pdf_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## 1) Importation et verification des données

Les données on été importé le 13/12/2021. 

```{r}
path = "./weekly_in_situ_co2_mlo.csv"
data_url = "https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv"
```

On verifie si les données ont été téléchargé, si ce n'est pas le cas on les téléchage
```{r}
if (file.exists(path)==F) {
  download.file(data_url,path)
}

```
On importe les données, en regardant le jeux de données, on se rend compte que les 44 première lignes expliquent le tableau, il faut donc les enlever lors de l'importation. 
```{r}
data = read.csv(path, skip = 44, header = F, dec = ".", sep = ",")
```

Le tableau contient 2 colones:

- La date 
- La concentration en C0² en micromol/mol

Les colonnes n'ont pas de titres donc on en ajoute un pour plus de clarté dans les analyses.
```{r}
names(data)=c("Date","C_ppm")
names(data)
```

On verifie que les tableau vas bien de 1958 à 2021.

```{r}
head(data)
tail(data)
```
Les mesures vont bien de Mars 1958 à septembre 2021.

Maintenant il faut verifier si des données sont manquantes.

```{r}
na_records = apply(data, 1, function (x) any(is.na(x)))
data[na_records,]
```
Il n'y en a pas, bonne nouvelle pour notre anlayse.

Il faut maintenant verifier la classe de chaque colones.

```{r}
str(data)
```
La colonne `date` est un vecteur chaine de caractère et la colone concentration (`C_ppm`) est un vecteur de nombre.

Pour pouvoir travailler sur les données il faut transformer les chaines de caractères en date.

```{r}
library(parsedate)
convert_week = function(iso) {
	as.character(parse_iso_8601(iso))
}

data$Date = as.Date(convert_week(data$Date))
```

```{r}
class(data$Date)
```
Les poitns sont dans l'ordre chronoligique donc pas besoin de les trier.

On verifie que les semaines font bien 7 jours.

```{r}
all(diff(data$date) == 7)
```
## Représentation graphique:

On Fait une première représentaion graphique pour visualiser les données.

```{r}
plot(data, type = "l")
```
On observe une tendence positive avec une oscilation.

Zoom sur 2ans

```{r}
plot(tail(data,100), type = "l")
```


Pour visuliser l'evolution globale on vas se concentrer sur les mesure anuelles.

Les oscilation sont anuelles et le pic semble au printemps.
Pour eviter de demarer une année dans le pic on va commencer l'année au 1er aout.

```{r}
pic_annuel = function(annee) {
      debut = paste0(annee-1,"-08-01")
      fin = paste0(annee,"-08-01")
      semaines = data$Date > debut & data$Date <= fin
      return(c(debut,fin,mean(data$C_ppm[semaines]),sd(data$C_ppm[semaines])))
}

annees = 1959:2021

annee = NULL
dates = NULL
mean_CO2 = NULL
sd_CO2 = NULL

for (i in annees) {
  annee = c(annee, rep(i,2))
  CO2 = pic_annuel(i)
  dates = c(dates,CO2[1],CO2[2])
  mean_CO2 = c(mean_CO2,rep(CO2[3],2))
  sd_CO2 = c(sd_CO2,rep(CO2[4],2))
}
inc_annuelle = data.frame(annee = annee, dates, mean_CO2 ,sd_CO2)
str(inc_annuelle)
inc_annuelle$dates = as.Date(convert_week(inc_annuelle$dates))
inc_annuelle$mean_CO2 = as.numeric(inc_annuelle$mean_CO2)
inc_annuelle$sd_CO2 = as.numeric(inc_annuelle$sd_CO2)
```


On refait le graphique en ajoutant par dessus la tendance.

```{r}
plot(data, type = "l")

points(inc_annuelle$mean_CO2~inc_annuelle$dates, col = "red", type = "l")
#points((inc_annuelle$mean_CO2-inc_annuelle$sd_CO2)~inc_annuelle$dates, col = "blue", type = "l")
#points((inc_annuelle$mean_CO2+inc_annuelle$sd_CO2)~inc_annuelle$dates, col = "blue", type = "l")
```

## Etude de la tendance

Pour estimer la tendance on peut utiliser une régression linéaire.

Il faut d'abord créer le modèle.
```{r}
m = lm(inc_annuelle$mean_CO2~inc_annuelle$dates)
```

```{r}
summary(m)

```
Notre modèle est bon, il explique 98% des données (R²=0.98).
La droite à pour intercept 3.244e+02 et pour pente 4.369e-3


```{r}
plot(data, type = "l", xlim = c(as.Date("1958-01-01"),as.Date("2030-12-31")), ylim = c(300,450))

abline(a = 3.244e+02, b = 4.369e-03, col = "blue")

abline(v=as.Date("2025-01-01"), lty=2)
abline (412, 0, lty = 2)
```
Maintenant on vas tenter d'extrapoler avec notre modèle à l'horizon 2025.

On peut extrapoler grâce au graphique qu'a l'horizon 2025 la quantité de CO2 sera d'environ 412 mircomol par mol.
On cette valeur est déjà dépassé en 2021 ce qui montre que notre modèle devrait être affiné pour ameliorer nos predictions.