Redirection de la sortie =stderr= vers la sortie =stdout= afin que la première soit imprimée correctement (dans le rapport, j'ai trouvé comment faire ça sur le [[http://kitchingroup.cheme.cmu.edu/blog/2015/01/04/Redirecting-stderr-in-org-mode-shell-blocks/][Kitchin's blog]]) :
Je peux afficher les lignes allant de la 39e à la 44e -- afin d'avoir des informations sur l'organisation des données dans le fichier -- avec le programme =sed= (j'ai trouvé les « bonnes lignes » du fichier en explorant celui-ci semi-interactivement avec le [[https://fr.wikipedia.org/wiki/Visionneur_(Terminal)][visionneur]], ou /pager/, =less=) :
** Statistiques de bases et recherche de données aberrantes
Je vais utiliser [[http://www.gnuplot.info/][gnuplot]] pour faire les figures et même un peu plus. Ici, il y a quelques « spécificités » du fichier de données dont je dois tenir compte :
- les données proprement dites commencent en ligne 45, je vais donc utiliser =skip 44= lorsque la commande =plot= de =gnuplot= sera appelée ;
- la virgule est utilisée pour séparer les colonnes du fichier, d'où le =set datafile separator ","= ci-dessous ;
- l'abscisse correspond à une date ce que je précise par =set xdata time= ;
- la date sera obtenue à partir de la première (année) et de la deuxième (mois) colonne, ce que je spécifie avec =set timefmt "%Y-%m-%d"= ;
- je souhaite avoir les labels de l'abscisse sous forme "mois-année" ce que j'obtiens avec =set xtics format "%m-%Y"= ;
- le reste devrait être assez clair.
Avant d'utiliser concrètement =gnuplot=, j'inscris dans ce document la version utilisée afin d'améliorer la reproductibilité de cette analyse :
#+NAME: version-de-gnuplot
#+BEGIN_SRC sh :exports both :results output
gnuplot --version
#+END_SRC
#+RESULTS: version-de-gnuplot
: gnuplot 5.2 patchlevel 4
Ceci étant fait, je peux utiliser la fonction =stats= de =gnuplot= pour accéder à un résumé des données :
#+NAME: co2_mlo-résumé
#+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export :results output
set datafile separator ","
stats "weekly_in_situ_co2_mlo.csv" skip 44 using 1:2
#+END_SRC
#+RESULTS: co2_mlo-résumé
#+BEGIN_EXAMPLE
FILE:
Records: 3082
Out of range: 0
Invalid: 0
Column headers: 0
Blank: 0
Data Blocks: 1
COLUMNS:
Mean: 1988.3303 354.0582
Std Dev: 17.2668 27.2089
Sample StdDev: 17.2696 27.2134
Skewness: -0.0109 0.3316
Kurtosis: 1.8177 1.9420
Avg Dev: 14.9273 23.3342
Sum: 6.12803e+06 1.09121e+06
Sum Sq.: 1.21855e+10 3.88633e+08
Mean Err.: 0.3110 0.4901
Std Dev Err.: 0.2199 0.3466
Skewness Err.: 0.0441 0.0441
Kurtosis Err.: 0.0882 0.0882
Minimum: 1958.0000 [ 0] 313.0400 [ 60]
Maximum: 2018.0000 [3081] 411.8100 [3068]
Quartile: 1974.0000 329.4100
Median: 1988.0000 351.6000
Quartile: 2003.0000 375.7800
Linear Model: y = 1.559 x - 2745
Slope: 1.559 +- 0.004175
Intercept: -2745 +- 8.301
Correlation: r = 0.9891
Sum xy: 2.171e+09
#+END_EXAMPLE
À ce stade il n'y a pas de signe de présence de données aberrantes.
* Graphe global :export:
La figure montrant la suite chronologique de la concentration de $CO_2$ en fonction de la date est obtenue avec :
#+NAME: co2_mlo-fig
#+HEADERS: :file co2_mlo.png
#+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export
set datafile separator ","
set xdata time
set timefmt "%Y-%m-%d"
set xtics format "%m-%Y" rotate
set grid
set xrange [*:*]
set yrange [300:425]
set title '[CO_{2}] atmosphérique (observatoire de Mauna Loa, Hawaï)'
set xlabel 'Date'
set ylabel '[CO_{2}] (ppm)'
unset key
plot "weekly_in_situ_co2_mlo.csv" skip 44 \
using 1:2 with lines linecolor 'red'
#+END_SRC
#+RESULTS: co2_mlo-fig
[[file:co2_mlo.png]]
* Estimer la tendance :export:
** Ajustement d'un modèle paramétrique
Pour l'estimation de la tendance, je dois d'une manière ou d'une autre « supprimer » les oscillations périodiques annuelles dues à la pousse des feuilles des feuillus (phase décroissante) puis à leur chute et décomposition (phase croissante), comme expliqué dans la rubrique [[http://scrippsco2.ucsd.edu/history_legacy/early_keeling_curve][History & Legacy]] du site sur lequel se trouvent les données (c'est ce qui se passe dans l'hémisphère nord qui domine du fait de la plus grande surface de forêts à feuilles caducs qu'on y trouve). Une façon de « supprimer » des oscillations périodiques à fréquence relativement élevée (nous avons des oscillations dont la période est l'année alors que nous disposons de 60 ans d'observations) est d'intégrer le signal (l'intégrale d'une fonction sinus ou cosinus sur une période est nulle). C'est ce que je fais maintenant en prenant comme unité de temps l'année (52 semaines) et en soustrayant la « ligne de base » égale à 313 (d'après les statistiques obtenues ci-dessus) :
#+NAME: co2_mlo-cumsum-fig-1
#+HEADERS: :file co2_mlo_cumsum_1.png
#+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export
set datafile separator ","
set grid
set title '[CO_{2}] sommée (observatoire de Mauna Loa, Hawaï)'
Sur ces données intégrées, les oscillations ont disparu et la courbe a un aspect bien plus lisse que sur les données de départ. Je vais à présent utiliser la commande =fit= de =gnuplot= pour ajuster par la méthode des moindres carrés un polynôme de degré 4 aux données intégrées. J'ai choisi le degré du polynôme est les valeurs initiales des paramètres ci-dessous en répétant la procédure suivante en mode interactif:
#+NAME: co2_mlo-cumsum-fit
#+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export
set datafile separator ","
s = 0
h(x) = a*x + b*x**2 +c*x**3+d*x**4
a=1
b=12.5
c=0.32
d=0.0005
fit h(x) "weekly_in_situ_co2_mlo.csv" skip 44 \
using ($0/52):(s=s+($2-valeur_initiale)) via a,b,c,d
with lines linecolor 'red' linewidth 2 title '[CO_{2}] intégrée', \
h(x) linecolor 'black' linewidth 1 dashtype '-' title 'Polynôme de degré 4'
#+END_SRC
#+RESULTS: co2_mlo-cumsum-fig
[[file:co2_mlo_cumsum.png]]
** Estimation d'un modèle non-paramétrique
Une autre façon d'estimer la tendance consiste à travailler année par année et, en utilisant le fait que la croissance / dégénérescence des feuilles caducs -- responsable des oscillations annuelles -- est de somme nulle sur l'année, à intégrer sur chaque année les concentrations. On peut alors approcher la vraie tendance par une fonction constante par morceaux -- la constante de chaque année étant la moyenne calculée sur cette même année -- ou par une fonction linéaire par morceaux -- la pente étant constante sur l'année. Pour faire ces calculs, j'ai besoin d'outils un petit peu plus sophistiqués que ceux mis à ma disposition par =gnuplot= et je vais employer =Python=.
*** Lecture des données dans =Python=
Je lis les données ligne pas ligne en faisant attention au fait que les 44 premières lignes du fichier contiennent des méta-données. Je crée trois [[https://en.wikipedia.org/wiki/Python_syntax_and_semantics#Lists,_tuples,_sets,_dictionaries][listes]] :
- la_date :: contient la date sous forme de chaîne caractères, telle qu'elle se trouve dans le fichier ;
- annee :: contient l'année (un entier) où chaque mesure a été effectuée (nécessaire pour calculer la moyenne annuelle) ;
- CO2 :: contient la concentration (un réel) de CO$_{2}$.
Pour « m'assurer » que l'importation s'est bien passée j'imprime les dix premières et les dix dernières valeurs de [CO$_{2}$] pour les comparer avec les valeurs précédentes (affichées avec =sed= et =tail=) :
#+NAME: lecture-des-données-dans-Python
#+BEGIN_SRC python :session *Mauna-Loa* :results output :exports both :eval no-export
f = open("weekly_in_situ_co2_mlo.csv","r")
titi = [f.readline() for i in range(44)]
toto = f.readline()
CO2 =[]
annee = []
la_date = []
while toto:
la_date.append(toto.split(',')[0])
CO2.append(float(toto.split(',')[1]))
annee.append(int(toto.split('-')[0]))
toto = f.readline()
f.close()
print("\nLes dix premières valeurs de [CO2] sont:\n")
print(CO2[:10])
print("Les dix dernières valeurs de [CO2] sont:\n")
J'utilise =gnuplot= pour faire le graphe de la tendance basée sur l'approximation constante par morceau (moyenne annuelle) :
#+NAME: co2_mlo-tendance-moyenne-fig
#+HEADERS: :file co2_mlo_tendance-moyenne.png
#+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export
set datafile separator ","
set xdata time
set timefmt "%Y-%m-%d"
set xtics format "%m-%Y" rotate
set grid
set xlabel 'Date'
set ylabel '[CO_{2}] moyenne sur une année (ppm)'
unset key
plot "co2_mlo_hebdo_non_para.csv" using 1:2 \
with lines lc 'red' lw 3
#+END_SRC
#+RESULTS: co2_mlo-tendance-moyenne-fig
[[file:co2_mlo_tendance-moyenne.png]]
De même, je construis le graphe de la tendance basée sur l'approximation linéaire par morceaux (les estimations des deux années extrêmes ne sont pas fiables puisque les premiers mois de l'année 1958 manquent, de même que les derniers mois de l'année 2018) :
#+NAME: co2_mlo-tendance-lineaire-fig
#+HEADERS: :file co2_mlo_tendance-lineaire.png
#+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export
set datafile separator ","
set xdata time
set timefmt "%Y-%m-%d"
set xtics format "%m-%Y" rotate
set grid
set xlabel 'Date'
set ylabel '[CO_{2}] approx. linéaire de tendance (ppm)'
unset key
plot "co2_mlo_hebdo_non_para.csv" using 1:3 \
with lines lc 'red' lw 3
#+END_SRC
#+RESULTS: co2_mlo-tendance-lineaire-fig
[[file:co2_mlo_tendance-lineaire.png]]
Si nous comparons ces données avec la première figure de la page Wikipédia [[https://fr.wikipedia.org/wiki/Ressources_et_consommation_%C3%A9nerg%C3%A9tiques_mondiales][Ressources et consommation énergétiques mondiales]], tirée du [[https://www.bp.com/en/global/corporate/energy-economics/statistical-review-of-world-energy.html][BP Statistical Review of World Energy]], nous constatons que la relation entre [CO$_2$] (mesurée à Mauna Loa) et consommation d'énergie fossile au niveau mondiale est loin d'être simple (en principe il faudrait intégrer la consommation, si nous supposons que tout le carbone rejeté s'accumule dans l'atmosphère). Les doubles chocs pétroliers des années 70 (à la suite desquels la pente de la consommation de pétrole diminue nettement) ne sont suivis d'aucun effet sur la [CO$_2$] -- on devine à peine un léger coude après 1973. La forte accélération de la consommation de charbon après 2000 ne se voit pas non plus sur la [CO$_2$]. Il faut clairement prendre en compte le rôle de tampon / réserve de CO$_2$ joué par les océans (qui deviennent plus [[https://fr.wikipedia.org/wiki/Acidification_des_oc%C3%A9ans][acides]]), mais là, ça se complique...
Enfin, « l'inconvénient » de ces méthodes non-paramétriques locales est qu'elles ne permettent pas de faire d'extrapolation (sensée).
* Estimer les oscillations :export:
** Avec l'approximation constante par morceaux de la tendance
Je construis le graphe des résidus de [CO$_2$] obtenu en soustrayant l'approximation constante par morceaux de la tendance :
#+NAME: co2_mlo-oscillation-moyenne-fig
#+HEADERS: :file co2_mlo_oscillation-moyenne.png
#+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export
set datafile separator ","
set xdata time
set timefmt "%Y-%m-%d"
set xtics format "%m-%Y" rotate
set grid
set title 'Variation annuelle de [CO_{2}] (tendance approx. par fct cst par morceaux)'
set xlabel 'Date'
set ylabel '[CO_{2}] (ppm)'
unset key
plot "co2_mlo_hebdo_non_para.csv" using 1:4 \
with lines lc 'red' lw 2
#+END_SRC
#+RESULTS: co2_mlo-oscillation-moyenne-fig
[[file:co2_mlo_oscillation-moyenne.png]]
** Avec l'approximation linéaire par morceaux de la tendance
Je construis le graphe des résidus de [CO$_2$] obtenu en soustrayant l'approximation linéaire par morceaux de la tendance :
#+NAME: co2_mlo-oscillation-lineaire-fig
#+HEADERS: :file co2_mlo_oscillation-lineaire.png
#+BEGIN_SRC gnuplot :exports both :session *gnuplot* :eval no-export
set datafile separator ","
set xdata time
set timefmt "%Y-%m-%d"
set xtics format "%m-%Y" rotate
set grid
set title 'Variation annuelle de [CO_{2}] (tendance approx. par fct linéaire par morceaux)'
set xlabel 'Date'
set ylabel '[CO_{2}] (ppm)'
unset key
plot "co2_mlo_hebdo_non_para.csv" using 1:5 \
with lines lc 'red' lw 2
#+END_SRC
#+RESULTS: co2_mlo-oscillation-lineaire-fig
[[file:co2_mlo_oscillation-lineaire.png]]
Les deux graphes (surtout le second) suggèrent une augmentation de l'amplitude des oscillations.
