--- title: 'Traitement sujet 1 : Concentration de CO~2~ dans l''atmosphère depuis 1958' author: "David Pinaud" date: "2023-09-12" output: pdf_document: df_print: paged --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # 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). # 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) ``` # 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. ## Tendance à long terme 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) ``` 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 : 1. 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 ; 2. une périodicité de l'ordre d'un an, probablement expliquée par des variations saisonnières dans les émissions de CO~2~.