Etude de la concentration de CO2 dans l’atmosphère depuis 1958

Fichier télécharger le 18/02/2021

Atmospheric CO2 concentrations (ppm) derived from in situ air measurements at Mauna Loa, Observatory, Hawaii: Latitude 19.5°N Longitude 155.6°W Elevation 3397m.

Source: R. F. Keeling, S. J. Walker, S. C. Piper and A. F. Bollenbacher. Scripps CO2 Program ( http://scrippsco2.ucsd.edu ). Scripps Institution of Oceanography (SIO). University of California. University of California. La Jolla, California USA 92093-0244.

Données hebdomadaires.

Les analyses effectuées pour décomposer la série temporelle ont été réalisées à partir du livre Numerical Ecology écrit par Pierre et Louis Lengendre et édité par Elsevier en 1998.

url_data_co2_web = lien pour télécharger les données. url_data_co2_local = lien des données en local.

url_data_co2_web = "https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv"

url_data_co2_locale = "D:/PhD formations/formation recherche reproductible/mooc-rr/module3/exo3/weekly_in_situ_co2_mlo.csv"

1. Réalisez un graphique qui vous montrera une oscillation périodique superposée à une évolution systématique plus lente.

Après un examen visiuel du fichier csv, l’on constate que les données commencent à la ligne 45. L’on peut donc “skipper” les 44 premières lignes. La première colonne correspond aux dates et la seconde aux concentration de CO2 en ppm.

data_co2 = read.csv(url_data_co2_locale, skip = 44, col.names = c("dates","CO2"), header = FALSE)

On vérifie le type de données dans chaque colonne.

class(data_co2$dates)
## [1] "factor"
class(data_co2$CO2)
## [1] "numeric"

Il faut convertir le vecteur dates au format date avec la fonction as.Date. Les dates renseignées étant déjà au bon format, cela ne devrait pas poser de problème. Les résultats sont placés dans un vecteur dates2:

data_co2$dates2 = as.Date(data_co2$dates, format = "%Y-%m-%d")

On vérifie que cela a fonctionné:

class(data_co2$dates2)
## [1] "Date"

Y a-t-il des points manquants dans nos données ?

na_records = apply(data_co2, 1, function (x) any(is.na(x)))
data_co2[na_records,]
## [1] dates  CO2    dates2
## <0 rows> (or 0-length row.names)

On peut à présent représenter graphiquement la relation entre la concentration en CO2 et le temps:

plot(data_co2$CO2 ~ data_co2$dates2, type = "l", col = "blue", xlab = "Time", ylab = "CO2 in ppm")

On zoom sur la fin des données pour mieux voir la périodicité.

plot(data_co2$CO2[c(3000:3208)] ~ data_co2$dates2[c(3000:3208)], type = "l", col = "blue", xlab = "Time", ylab = "CO2 in ppm")

On peut donc constater la présence d’une périodicité, qui semble à première vue être annuelle, et d’une tendance générale à l’augmentation.

2. Séparez ces deux phénomènes. Caractérisez l’oscillation périodique. Proposez un modèle simple de la contribution lente, estimez ses paramètres et tentez une extrapolation jusqu’à 2025 (dans le but de pouvoir valider le modèle par des observations futures).

Pour décomposer cette série, deux modèles de décomposition s’offrent à nous: un modèle additif ou multiplicatif. La périodicité étant annuelle, nous allons estimer les moyennes et les écarts types annuels. Si aucune relation n’apparait entre ces deux paramètres, on appliquera un modèle additif pour décomposer la série, sinon on appliquera un modèle multiplicatif.

On charge le package lubridate qui va nous permettre d’utiliser la fonction year que nous allons utiliser pour isoler les concentrations en CO2 qui appartiennent toutes à la même année.

On créé deux vecteurs vides, moyenne_ann et et_ann pour stocker la concentration en CO2 moyenne de chaque année et l’écart type associé

moyenne_ann = rep(NA, length(min(year(data_co2$dates2)):max(year(data_co2$dates2))))
et_ann = rep(NA, length(min(year(data_co2$dates2)):max(year(data_co2$dates2))))

On effectue ensuite une boucle de 1958 (année minimale (min(year(data_co2$dates2)))) à 2021 (année maximale (max(year(data_co2$dates2)))). On définit un conteur a égal à 1 dont la valeur changera à chaque itération pour nous permettre de stocker les moyennes et les écarts types à la bonne position dans nos vecteurs vides moyenne_ann et et_ann. On définit également un vecteur logique test qui va nous permettre d’isoler dans la matrice data_co2$CO2 les concentration en CO2 correspondant à l’année de l’itérations.

a=1

for (i in min(year(data_co2$dates2)):max(year(data_co2$dates2))){
  test = (year(data_co2$dates2) == i) == TRUE
  
  moyenne_ann[a] = mean(data_co2$CO2[test == TRUE])
  et_ann[a] = sd(data_co2$CO2[test == TRUE])
  a=a+1
}

On peut à présent plotter les résultats:

plot(moyenne_ann ~ et_ann, xlab = "annual mean", ylab = "standard deviation")

Et calculer le coefficient de corrélation de Pearson r entre les deux paramètres avec la fonction cor.

r = cor(moyenne_ann, et_ann, method = "pearson")
r
## [1] 0.02799553

Il n’y a pas de relation entre les deux variables, on peut donc utiliser un modèle additif pour décomposer la série.

Pour extraire la tendance des données, on applique une moyenne mobile simple d’ordre n = 52 (nombre de semaines dans une année) avec la fonction movavg du package pracma

library(pracma)
trend = movavg(data_co2$CO2, 52, type = "s")

On plot le résultat

plot(trend ~ data_co2$dates2, type = "l", col = "red", xlab = "", ylab = "", axes = FALSE)
par(new = T) 
plot(data_co2$CO2 ~ data_co2$dates2, type = "l", col = "blue", xlab = "Time", ylab = "CO2 in ppm")
legend(1, 400, legend = c("tendance", "raw data"), col = c("red", "blue"), lty = 1, cex = 1)

On extrait ensuite la tendance des données avec la formule suivante \(Y_{i Detrend} = Y_{i} - \hat{Y_{i}}\), avec \(Y_{i}\) les concentrations en CO2 et \(\hat{Y_{i}}\) la tendance estimé avec la moyenne mobile.

data_co2$CO2dt = data_co2$CO2 - trend

On plot le résultat pour voir si cela a marché.

plot(data_co2$CO2dt ~ data_co2$dates2, type = "l", col = "blue", xlab = "Time", ylab = "Detrend CO2 in ppm")

Pour caractériser l’oscillation périodique, nous allons effectué une mesure d’autocorrelation avec la fonction acf sur les données détentancées:

auto_cor = acf(data_co2$CO2dt, type = "correlation", lag.max = 100)

On observe un décalage d’environs 52-53 pas de temps, ce qui correspond à la durée d’une année en semaine. Nous avons donc une périodicité annuelle.

On va ensuite ajuster un polynome de degré 3 avec la fonction polyfit du package pracma.

model = polyfit(seq(1,3208), data_co2$CO2, 3)

Voici les paramètres du modèle:

model
## [1] 4.024541e-10 2.900286e-06 1.787978e-02 3.143707e+02

Et sa représentation graphique.

output = polyval(model, seq(1,3208)) 
plot(output ~ data_co2$dates2, type = "l", col = "red", axes = FALSE, xlab = "", ylab = "", lwd = 2)
par(new = T) 
plot(data_co2$CO2 ~ data_co2$dates2, type = "l", col = "blue", xlab = "Time", ylab = "CO2 in ppm")
legend(1, 400, legend = c("model", "raw data"), col = c("red", "blue"), lty = 1, cex = 1)

Malheureusement, nous ne sommes pas parvenu à ajuster le polynome en considérant les dates en tant que tel. La dernière valeur de CO2 enregistrée correspond à la date:

last_date = data_co2$dates2[length(data_co2$dates2)]

Pour extrapoler jusqu’en 2025, il nous faut estimer le nombre de semaine entre le 30 janvier 2021 et le 30 janvier 2025.

date1 = as.Date("2021-01-30")
date2 =  as.Date("2025-01-30")
difftime(date2, date1, units="weeks")
## Time difference of 208.7143 weeks

En arrondissant à 209, on peut relancer notre code précédent pour extrapoler jusqu’en 2025:

output2 = polyval(model, seq(1,3208+209))

Avec le package lubridate et la fonction weeks l’on va pouvoir compléter notre vecteur dates2 avec les dates manquantes entre 30-01-2021 et le 30-01-2025.

library(lubridate)
date_new = c(data_co2$dates2,last_date+weeks(1:209))
last_date+weeks(1:209)
##   [1] "2021-02-06" "2021-02-13" "2021-02-20" "2021-02-27" "2021-03-06"
##   [6] "2021-03-13" "2021-03-20" "2021-03-27" "2021-04-03" "2021-04-10"
##  [11] "2021-04-17" "2021-04-24" "2021-05-01" "2021-05-08" "2021-05-15"
##  [16] "2021-05-22" "2021-05-29" "2021-06-05" "2021-06-12" "2021-06-19"
##  [21] "2021-06-26" "2021-07-03" "2021-07-10" "2021-07-17" "2021-07-24"
##  [26] "2021-07-31" "2021-08-07" "2021-08-14" "2021-08-21" "2021-08-28"
##  [31] "2021-09-04" "2021-09-11" "2021-09-18" "2021-09-25" "2021-10-02"
##  [36] "2021-10-09" "2021-10-16" "2021-10-23" "2021-10-30" "2021-11-06"
##  [41] "2021-11-13" "2021-11-20" "2021-11-27" "2021-12-04" "2021-12-11"
##  [46] "2021-12-18" "2021-12-25" "2022-01-01" "2022-01-08" "2022-01-15"
##  [51] "2022-01-22" "2022-01-29" "2022-02-05" "2022-02-12" "2022-02-19"
##  [56] "2022-02-26" "2022-03-05" "2022-03-12" "2022-03-19" "2022-03-26"
##  [61] "2022-04-02" "2022-04-09" "2022-04-16" "2022-04-23" "2022-04-30"
##  [66] "2022-05-07" "2022-05-14" "2022-05-21" "2022-05-28" "2022-06-04"
##  [71] "2022-06-11" "2022-06-18" "2022-06-25" "2022-07-02" "2022-07-09"
##  [76] "2022-07-16" "2022-07-23" "2022-07-30" "2022-08-06" "2022-08-13"
##  [81] "2022-08-20" "2022-08-27" "2022-09-03" "2022-09-10" "2022-09-17"
##  [86] "2022-09-24" "2022-10-01" "2022-10-08" "2022-10-15" "2022-10-22"
##  [91] "2022-10-29" "2022-11-05" "2022-11-12" "2022-11-19" "2022-11-26"
##  [96] "2022-12-03" "2022-12-10" "2022-12-17" "2022-12-24" "2022-12-31"
## [101] "2023-01-07" "2023-01-14" "2023-01-21" "2023-01-28" "2023-02-04"
## [106] "2023-02-11" "2023-02-18" "2023-02-25" "2023-03-04" "2023-03-11"
## [111] "2023-03-18" "2023-03-25" "2023-04-01" "2023-04-08" "2023-04-15"
## [116] "2023-04-22" "2023-04-29" "2023-05-06" "2023-05-13" "2023-05-20"
## [121] "2023-05-27" "2023-06-03" "2023-06-10" "2023-06-17" "2023-06-24"
## [126] "2023-07-01" "2023-07-08" "2023-07-15" "2023-07-22" "2023-07-29"
## [131] "2023-08-05" "2023-08-12" "2023-08-19" "2023-08-26" "2023-09-02"
## [136] "2023-09-09" "2023-09-16" "2023-09-23" "2023-09-30" "2023-10-07"
## [141] "2023-10-14" "2023-10-21" "2023-10-28" "2023-11-04" "2023-11-11"
## [146] "2023-11-18" "2023-11-25" "2023-12-02" "2023-12-09" "2023-12-16"
## [151] "2023-12-23" "2023-12-30" "2024-01-06" "2024-01-13" "2024-01-20"
## [156] "2024-01-27" "2024-02-03" "2024-02-10" "2024-02-17" "2024-02-24"
## [161] "2024-03-02" "2024-03-09" "2024-03-16" "2024-03-23" "2024-03-30"
## [166] "2024-04-06" "2024-04-13" "2024-04-20" "2024-04-27" "2024-05-04"
## [171] "2024-05-11" "2024-05-18" "2024-05-25" "2024-06-01" "2024-06-08"
## [176] "2024-06-15" "2024-06-22" "2024-06-29" "2024-07-06" "2024-07-13"
## [181] "2024-07-20" "2024-07-27" "2024-08-03" "2024-08-10" "2024-08-17"
## [186] "2024-08-24" "2024-08-31" "2024-09-07" "2024-09-14" "2024-09-21"
## [191] "2024-09-28" "2024-10-05" "2024-10-12" "2024-10-19" "2024-10-26"
## [196] "2024-11-02" "2024-11-09" "2024-11-16" "2024-11-23" "2024-11-30"
## [201] "2024-12-07" "2024-12-14" "2024-12-21" "2024-12-28" "2025-01-04"
## [206] "2025-01-11" "2025-01-18" "2025-01-25" "2025-02-01"

Et représenter graphiquement le résultat.

plot(output2 ~ date_new, type = "l", col = "red", xlab = "", ylab = "", lwd = 2)
par(new = T) 
plot(c(data_co2$CO2, rep(NA, length(last_date+weeks(1:209)))) ~ date_new, type = "l", col = "blue", axes = FALSE, xlab = "Time", ylab = "CO2 in ppm")
legend(1, 400, legend = c("model", "raw data"), col = c("red", "blue"), lty = 1, cex = 1)

Notre modèle sous estime beaucoup l’accroissement de la concentration en CO2.