Commit 6168f415 authored by David Pinaud's avatar David Pinaud

final version subject #1 from David Pinaud

parent 821f83d9
--- ---
title: "Votre titre" title: 'Traitement sujet 1 : Concentration de CO~2~ dans l''atmosphère depuis 1958'
author: "Votre nom" author: "David Pinaud"
date: "La date du jour" date: "2023-09-12"
output: html_document output:
pdf_document:
df_print: paged
--- ---
...@@ -10,24 +12,134 @@ output: html_document ...@@ -10,24 +12,134 @@ output: html_document
knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(echo = TRUE)
``` ```
## Quelques explications
Ceci est un document R markdown que vous pouvez aisément exporter au format HTML, PDF, et MS Word. Pour plus de détails sur R Markdown consultez <http://rmarkdown.rstudio.com>. # Objectif
Il est de présenter un document reproductible réalisé selon les modes énoncés dans le MOOC. Le sujet choisi est le n°1 _Concentration de CO~2~ dans l'atmosphère depuis 1958_ avec des mesures de la concentration de CO~2~ dans l'atmosphère effectuées à l'observatoire de Mauna Loa, Hawaii, États-Unis.
Les différents points à aborder pour l'exercice sont :
1. Réaliser un graphique qui montre une oscillation périodique superposée à une évolution systématique plus lente.
2. Séparer ces deux phénomènes, en caractérisant d'une part l'oscillation périodique et d'autre part en proposant un modèle simple de la contribution lente pour pouvoir ensuite proposer une extrapolation de cette dernière évolution jusqu'à 2025 (dans le but de pouvoir valider le modèle par des observations futures).
Lorsque vous cliquerez sur le bouton **Knit** ce document sera compilé afin de ré-exécuter le code R et d'inclure les résultats dans un document final. Comme nous vous l'avons montré dans la vidéo, on inclue du code R de la façon suivante:
```{r cars}
summary(cars) # Données utilisées et importation
Les données sont disponibles sur le [site internet de l'institut Scripps]( https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.html), on télécharge le fichier
"weekly_in_situ_co2_mlo.csv" avec les observations hebdomadaires [ici](https://scrippsco2.ucsd.edu/data/atmospheric_co2/mlo.html), jusqu'à la date du 2023-09-12 (date du téléchargement). Ce fichier est en copie locale et sert maintenant de base aux analyses.
```{r download}
data_url <-
"https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv"
if(file.exists("weekly_in_situ_co2_mlo.csv"))
{
data <- read.csv("weekly_in_situ_co2_mlo.csv")
} else {
data <- read.csv(data_url, skip=44, header=FALSE, # on n'importe pas les 44emes premières lignes
col.names=c("WeekDate", "CO2")) # colonnes "WeekDate" et "CO2".
write.csv(data, "weekly_in_situ_co2_mlo.csv", row.names=FALSE)
}
str(data)
range(data$WeekDate)
```
Le fichier contient deux colonnes :
* `WeekDate`: la date de la semaine, ajustée à 12:00 pour le jour du milieu de chaque semaine
* `CO2`: la concentration en CO~2~ en $\mu$mole CO~2~ par mole (ppm).
# Mise en forme des données et vérification
La colonne `WeekDate` est en format `character`, nous allons la transformer en format `Date`, puis classer le `data.frame` selon cette variable par ordre croissant.
```{r}
library(parsedate)
data$WeekDate <- as.Date(data$WeekDate)
data <- data[order(data$WeekDate),]
class(data$WeekDate)
```
La colonne `WeekDate` est bien en format `Date`.
```{r}
table(is.na(data$CO2))
```
Il n'y a pas de donnée manquante dans le tableau.
# Représentation graphique
On représente sur un graphique la série en entier et un extrait sur 3 ans.
```{r}
library(ggplot2)
library(patchwork)
p1 <- ggplot(data) +
geom_line(aes(WeekDate, CO2), col="blue") +
labs(title=bquote("Atmospheric " ~CO[2]~" concentration"),
subtitle = "Mauna Loa Observatory, Hawaii ; data from 1958 to 2023",
x="Date (weekly measures)",
y = bquote(~CO[2]~" concentration (ppm)"))
p1 <- p1 + annotate("rect", xmin=as.Date("2015-01-01"), xmax=as.Date("2017-12-31"),
ymin=396, ymax=412, fill="red", alpha=0.2)
subd <- data[data$WeekDate >= "2015-01-01" & data$WeekDate <= "2017-12-31",] # extrait sur 3 ans
p2 <- ggplot(subd) +
geom_line(aes(WeekDate, CO2), col="blue") +
geom_point(aes(WeekDate, CO2), size=0.8, fill="blue", shape=21, col="red") +
labs(x=NULL, y = NULL) +
theme(axis.text=element_text(size=5), axis.title=element_text(size=3))
p1 + inset_element(p2, 0.01, 0.5, 0.5, 0.99)
``` ```
Et on peut aussi aisément inclure des figures. Par exemple: # Décomposition de la série
La série montre une oscillation périodique superposée à une évolution systématique plus lente. Nous allons caractériser ces deux phénomènes.
```{r pressure, echo=FALSE} ## Tendance à long terme
plot(pressure) On teste d'abord un ajustement linéaire (modèle le plus simple).
```{r}
lm1 <- lm(CO2 ~ as.numeric(WeekDate), data)
plot(lm1, which=1)
``` ```
Vous remarquerez le paramètre `echo = FALSE` qui indique que le code ne doit pas apparaître dans la version finale du document. Nous vous recommandons dans le cadre de ce MOOC de ne pas utiliser ce paramètre car l'objectif est que vos analyses de données soient parfaitement transparentes pour être reproductibles. La dispersion des résidus indique une tendance non-linéaire. On ajuste alors un modèle quadratique $CO2 = a \times WeekDate^2 + b \times WeekDate + c$ qui devrait mieux coller aux données.
```{r}
lm2 <- lm(CO2 ~I(as.numeric(WeekDate)^2) + as.numeric(WeekDate), data)
plot(lm2, which=1)
knitr::kable(AIC(lm1, lm2))
```
Le modèle avec une tendance quadratique s'ajuste bien mieux aux données.
```{r}
knitr::kable(summary(lm2)$coefficients)
```
On représente alors ce modèle avec une prédiction jusqu'en milieu d'année 2025.
```{r}
newdata <- data.frame(WeekDate=seq.Date(from = as.Date("1958-03-01"),
to = as.Date("2025-07-01"), by = "weeks"))
pred <- predict(lm2, newdata = newdata, se=TRUE)
pred <- cbind(newdata, as.data.frame(pred)[,1:2])
ggplot() +
geom_line(data=data, aes(WeekDate, CO2), col="blue") +
geom_line(data=pred, aes(WeekDate, fit), col="black", linewidth=1.3) +
labs(title=bquote("Atmospheric " ~CO[2]~" concentration"),
subtitle = "Mauna Loa Observatory, Hawaii ; data from 1958 to 2023",
x="Date (weekly measures)",
y = bquote(~CO[2]~" concentration (ppm)"))
```
La concentration de CO~2~ dans l'atmosphère prédite par ce modèle pour la semaine du `r max(pred$WeekDate)` est de `r round(pred$fit[pred$WeekDate==max(pred$WeekDate)], 2)` ppm $\pm$ `r round(pred$se.fit[pred$WeekDate==max(pred$WeekDate)]*qt(0.975, df.residual(lm2)), 2)` (CI 95%).
## Variations périodiques
Nous allons travailler sur les résidus du modèle quadratique pour analyser la périodicité des variations à plus fine échelle temporelle à l'aide d'une fonction d'autocorrelation.
```{r}
data$resid.lm2 <- residuals(lm2, arg = "pearson")
acf(data$resid.lm2, lag.max=60, main = "Auto-correlation plot for residuals")
```
On observe une autocorrelation positive maximale autour du lag 52, soit 1 année (52 semaines). La périodicité de concentrations de CO~2~ dans l'atmosphère est donc d'environs 1 an.
# Conclusion
Les données de concentrations en CO~2~ dans l'atmosphère mesurées à l'observatoire de Mauna Loa (Hawaii) de 1958 à 2023 montrent des variations importantes, avec une superposition de deux phénomènes :
* une augmentation continue positive à long terme, avec une accélaration constante. En 2025, la concentration devrait atteindre `r round(pred$fit[pred$WeekDate==max(pred$WeekDate)], 2)` ppm ;
* une périodicité de l'ordre d'un an, probablement expliquée par des variations saisonnières dans les émissions de CO~2~.
Comme les résultats ne sont pas stockés dans les fichiers Rmd, pour faciliter la relecture de vos analyses par d'autres personnes, vous aurez donc intérêt à générer un HTML ou un PDF et à le commiter.
Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment