#+TITLE:  Concentration de CO2 dans l'atmosphère depuis 1958
#+AUTHOR: Miguel Arpa Perozo
#+DATE:   \today
#+LANGUAGE: fr
# #+PROPERTY: header-args :eval never-export
#+latex_header: \usepackage{minted}
#+startup: indent inlineimages

#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/readtheorg.css"/>
#+HTML_HEAD: <script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.3/jquery.min.js"></script>
#+HTML_HEAD: <script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.4/js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/lib/js/jquery.stickytableheaders.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/readtheorg/js/readtheorg.js"></script>

* Introduction
** Objectifs 
   1. Proposer un modèle de la variation de CO2 dans l'atmosphère en
      fonction du temps depuis 1958 à l'aide de différentes mesures.
   2. À l'aide du modèle, prédire la concentration future de CO2
      en 2025.
** Généralités sur les données utilisées

Voici quelques points généraux sur les données utilisées dans cette
étude : 
  - Le lien de téléchargement des données est le suivant : [[https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv][link]]
  - Les données ont été téléchargées le 2020-11-16 (date au format ISO). 
  - D'après le site de téléchargement, la concentration de CO2 est
    donnée en [micro-mol CO2 per more] (ppm).
  - Finalement, d'après le site : "The weekly values have been
    adjusted to 12:00 hours at middle day of each weekly period."
    
** Description des démarches adoptées. 
Deux démarches ont été adoptées afin de déterminer un modèle de la
concentration de C02 en fonction des mesures :

   1. Dans un premier temps, effectuer une analyse fréquentielle en
      utilisant la transformée de Fourier des mesures afin de
      caractériser certains phénomènes.
   2. Utiliser un outil d'optimisation pour déterminer les paramètres
      d'un modèle proposé par simple observation des mesures.

** Outils utilisés
1. Python 3.8.6 (GCC 10.2.0)
2. numpy 1.19.2
3. scipy 1.5.3
4. matplotlib 3.3.2

* Lecture ou téléchargement des données brutes
Les données ont été téléchargées le 2020-11-16. Le lien de
téléchargement utilisé est le suivant :

#+name: data-url
https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv

On télécharge les données si celles-ci ne sont pas déjà disponibles
localement.

#+BEGIN_SRC python :results silent :var data_url=data-url :session 
  data_file ="weekly_in_situ_co2_mlo.csv"
  import os
  import urllib.request

  if not os.path.exists(data_file):
      urllib.request.urlretrieve(data_url, data_file)
#+END_SRC

Les données commencent à la ligne 45, donc on ne prend pas en compte
les 45 premières lignes du fichier.

#+BEGIN_SRC python :results silent :var data_url=data-url :session 
  data = open(data_file, 'rb').read()
  lines = data.decode('latin-1').strip().split('\n')
  data_lines = lines[45:]
  table = [line.split(',') for line in data_lines]
#+END_SRC

On visualise les valeurs du tableau : 

#+BEGIN_SRC python :results value :session :exports both 
table[:5]
#+END_SRC

#+RESULTS:
| 1958-04-05 | 317.31 |
| 1958-04-12 | 317.69 |
| 1958-04-19 | 317.58 |
| 1958-04-26 | 316.48 |
| 1958-05-03 | 316.95 |

* Traitement et affichage des données. 

Dans cette partie, on convertit les données dans des formats qui
permettront de les traiter plus facilement.

** Conversion des /string/ en valeurs numériques. 
Les données dans ~table~ sont des string. On va convertir la première
colonne en objets ~datetime~ de Python, et la deuxième colonne en
~float~.

#+begin_src python :results silent :session :exports both
  import datetime
  convertedData = [(datetime.datetime.strptime(yearWeekDay, "%Y-%m-%d"), float(co2)) for yearWeekDay, co2 in table]

  dates = [dates for dates, concentration in convertedData]
  concentration = [concentration for dates, concentration in convertedData]
#+end_src

** Affichage des données. 
Ensuite, on convertit les ~datetime~ dans un format qui peut être utilisé
pour l'affichage des données, et on affiche les données brutes. 

#+begin_src python :results file :session :var matplot_lib_filename="rawData.png" :exports both
  import matplotlib.pyplot as plt
  import matplotlib.dates as pltDates 

  plotDates = pltDates.date2num(dates)

  plt.figure(figsize=(10,5))
  plt.plot_date(plotDates,concentration)
  plt.title("Concentration de CO2 en fonction du temps")
  plt.xlabel("Temps [années]")
  plt.ylabel("Concentration de CO2 [ppm]")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:rawData.png]]

On remarque sur la figure la superposition de deux phénomènes :
1. Un phénomène périodique de faible amplitude à haute fréquence. 
2. Une croissance globale, "lente" de forme parabolique. 

** /Zoom/ sur l'affichage des données 

Effectuons un /zoom/ sur une période de temps plus courte pour
caractériser les oscillations rapides observées précédemment.

#+begin_src python :results file :session :var matplot_lib_filename="zoomRawData.png" :exports both
  # Window's size where we want to zoom. 
  windowSize = 450; 
  plotDatesZoom = pltDates.date2num(dates[:windowSize])

  plt.figure(figsize=(10,5))
  plt.plot_date(plotDatesZoom,concentration[:windowSize])
  plt.title("Concentration de CO2 en fonction du temps")
  plt.xlabel("Temps [années]")
  plt.ylabel("Concentration de CO2 [ppm]")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:zoomRawData.png]]

Sur cette figure, on met en évidence le phénomène périodique sur le
niveau de concentration de CO2. On peut faire les premières hypothèses
suivantes :

1. La période des oscillations est de l'ordre d'une année.
2. L'amplitude varie entre cinq et 7 ppm. 

* Analyse fréquentielle de la variation de la concentration de CO2. 

Dans cette partie, on souhaite faire une analyse fréquentielle des
données en effectuant une transformée de Fourier sur les mesures afin
de mieux caractériser les phénomènes observés précédemment. Pour cela,
on utilise le module ~numpy~ et l'on affiche les résultats avec
~matplotlib~.

#+begin_src python :results file :session :var matplot_lib_filename="dataFFT.png" :exports both
  import numpy as np

  # Conversion of the data to numpy array
  co2 = np.array(concentration)

  #Remove mean value
  co2NoMean = co2 - np.mean(co2)
  # Concentration FFT
  concentrationFFT = np.fft.fft(co2NoMean)
  # We take the norm of the FFT
  concentrationFFTNorm  = np.absolute(concentrationFFT)
  concentrationFFTAngle = np.angle(concentrationFFT)

  # Frequency axis
  N  = co2.shape[0] # Number of samples
  Te = 1            # Sampling interval
  t  = np.arange(N)  # time reference

  #freq = np.arange(N)/(N*Te)
  freq = np.fft.fftfreq(N)*N*(1/N*Te)


  plt.figure(figsize=(10,5))
  plt.plot(freq,concentrationFFTNorm)
  plt.title("Spectre des données mesurés")
  plt.xlabel("Fréquence par semaine")
  plt.ylabel("Amplitude")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:dataFFT.png]]

Sur ce graphique, on peut voir l'ensemble du spectre des mesures. On
remarque quatres pics qui devraient correpondre au phénomène
périodique. On tentera dans la prochaine partie de caractériser plus
précisément ces pics.

** /Zoom/ sur le spectre du signal 

On effectue un zoom sur la partie positive du graphique pour tenter de
caractériser les différentes fréquences d'oscillations du signal.

#+begin_src python :results file :session :var matplot_lib_filename="fftZoom.png" :exports both
  startIndx = 10
  fftzoomIndx = int(N/6)

  plt.figure(figsize=(10,5))
  plt.plot(freq[startIndx:fftzoomIndx],concentrationFFTNorm[startIndx:fftzoomIndx])
  plt.title("Zoom sur le spectre des données mesurés")
  plt.xlabel("Fréquence par semaine")
  plt.ylabel("Amplitude")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:fftZoom.png]]

On remarque sur le spectre du signal que le phénomène périodique
possède deux pics d'amplitude d'environ 1400 et 4000. On suppose que
le premier pic avec une plus grande amplitude correspond au phénomène
périodique que l'on a mis en évidence dans l'observation des données
brutes. En effet, on a supposé précédemment que les oscillations
avaient une période d'une année, notre unité de temps étant une semaine,
on peut déduire la fréquence des oscillations en considérant qu'une
année a 52 semaine : \( f = \frac{1}{52} \approx 0.019 \) [/semaines] ce qui reste
cohérent avec les résultats dans le graphique.

** Filtrage de la contribution lente

On élimine toutes les valeurs supérieures au seuil établi
précédemment (4000) dans le spectre pour éliminer la
contribution lente dans le signal.

#+begin_src python :results file :session :var matplot_lib_filename="smallOsillations.png" :exports both
  # Extraction of the small oscillations
  smallAmplitude = 4000
  smallNorm = np.copy(concentrationFFTNorm)
  smallNorm[smallNorm > smallAmplitude] = 0

  plt.figure(figsize=(10,5))
  plt.plot(freq,smallNorm)
  plt.title("Spectre des données sans la contribution lente")
  plt.xlabel("Fréquence par semaine")
  plt.ylabel("Amplitude")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:smallOsillations.png]]

** Reconstruction du phénomène périodique
L'objectif est de déterminer l'amplitude des oscillations en prenant
en compte toutes les mesures et non pas juste sur une fenêtre
d'observation comme on l'a fait précédemment.

#+begin_src python :results file :session :var matplot_lib_filename="smallOsillationsTime.png" :exports both
  smallFreqSignal = smallNorm*np.exp(1j*concentrationFFTAngle)
  smallSignal = np.fft.ifft(smallFreqSignal)
  smallSignal = np.real(smallSignal)

  plt.figure(figsize=(10,5))
  plt.plot(t,smallSignal)
  plt.title("Reconstruction du phénomène périodique de la concentration de CO2 en fonction du temps")
  plt.xlabel("Temps [semaines]")
  plt.ylabel("Concentration de CO2 [ppm]")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:smallOsillationsTime.png]]

Effectuons un zoom sur le signal reconstruit : 

#+begin_src python :results file :session :var matplot_lib_filename="smallOsillationsTimeZoom.png" :exports both
  plt.figure(figsize=(10,5))
  plt.plot(t,smallSignal)
  plt.ylim(-10,10)
  plt.title("Reconstruction du phénomène périodique sur la concentration de CO2 en fonction du temps")
  plt.xlabel("Temps [semaines]")
  plt.ylabel("Concentration de CO2 [ppm]")
  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:smallOsillationsTimeZoom.png]]

On peut déduire que le phénomène périodique a une amplitude de cinq d'après le graphique. 

** Filtrage du phénomène périodique

On souhaite maintenant filtrer le phénomène périodique ce qui pourrait
nous aider a mieux comprendre la contribution lente. On élimine donc
toutes les valeurs inférieures au seuil établit précédemment.

#+begin_src python :results file :session :var matplot_lib_filename="bigOsillations.png" :exports both
  # Extraction of the big oscillations
  bigAmplitude = 4000
  bigNorm = np.copy(concentrationFFTNorm)
  bigNorm[bigNorm < bigAmplitude] = 0

  plt.figure(figsize=(10,5))
  plt.plot(freq,bigNorm)
  plt.title("Spectre des données sans le phénomène périodique")
  plt.xlabel("Fréquence par semaine")
  plt.ylabel("Amplitude")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:bigOsillations.png]]

** Reconstruction de la contribution lente du signal

Finalement on reconstruit le signal temporel puis on l'affiche afin de
mieux le caractériser.

#+begin_src python :results file :session :var matplot_lib_filename="bigOsillationsTime.png" :exports both
  bigFreqSignal = bigNorm*np.exp(1j*concentrationFFTAngle)
  bigSignal = np.fft.ifft(bigFreqSignal)
  bigSignal = np.real(bigSignal)

  plt.figure(figsize=(10,5))
  plt.plot(t,bigSignal)
  plt.title("Reconstruction du signal sans phénomène périodique sur la concentration de CO2 en fonction du temps")
  plt.xlabel("Temps [années]")
  plt.ylabel("Concentration de CO2 [ppm]")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:bigOsillationsTime.png]]

Ce graphique ne nous permet pas de déduire grand chose quant à cette
contribution lente. Tentons de le comparer avec les mesures.
** Comparaison des signaux filtrés et brut

#+begin_src python :results file :session :var matplot_lib_filename="bigOsillationsComparisons.png" :exports both
  plt.figure(figsize=(10,5))
  plt.plot(t,co2,label="measures")
  plt.plot(t,bigSignal+np.mean(co2),label="filtered")
  plt.legend()
  plt.title("Comparaisons des signaux filtré et brut")
  plt.xlabel("Temps [semaines]")
  plt.ylabel("Concentration de CO2 [ppm]")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:bigOsillationsComparisons.png]]

Malheuresement, cette approche ne nous permet pas d'avoir plus
d'informations pertinentes quant à cette contribution lente. Donc par
la suite, on modélisera cette contribution par un polynôme
d'ordre 2.

* Modélisation des phénomènes et extrapolations 
** Oscillations périodiques
L'analyse fréquentielle nous a permis caractériser le phénomène
périodique par un signal de fréquence \( f = 0.02 \) par unité de
temps (par semaine), et avec une amplitude de 5 [ppm]. Ainsi, soit \(s_p(t)\) le
signal périodique avec \(t\) en semaine, on a :

\begin{equation}
s_p(t) = 5 \cos (0.02 t)
\end{equation}

** Contribution lente du signal 
L'analyse fréquentielle ne nous a pas permis d'avoir plus
d'information quant à la contribution lente du signal. Cependant, si
l'on regarde l'allure générale des données, on peut supposer que le
concentration de CO2 varie de façon quadratique par rapport au
temps. Donc, on va modéliser cette contribution lente par un polynôme
d'ordre deux \( y(t) = at^2 +bt +c \). On utilisera un outil de /curve
fitting/ pour caractériser ce polynôme.

Compte tenu de la modélisation du phénomène périodique, celui-ci
possède une valeur moyenne nulle. Or les outils de /curve fitting/
calculent la norme deux de l'erreur entre le modèle et les données,
donc il est inutile de filtrer les données brute pour tenter d'enlever
le phénomène périodique pour caractériser la contribution lente.

Dans un premier temps, on définit une fonction de notre modèle pour la
contribution lente

#+begin_src python :results silent :session  :exports both
  def slowContributionModel(t,a,b,c):
      return a*t*t + b*t + c
#+end_src

Puis l'on utilise le module ~scipy~ plus spécifiquement le module
~scipy.optimize~ pour déterminer les paramètres de notre modèle.

#+begin_src python :results output :session  :exports both
  from scipy import optimize
  params, params_covariance = optimize.curve_fit(slowContributionModel,t,co2, p0=[1, 1, 1])
  print("Parameters :\t",params)
  print("Parameters covariance :\t", params_covariance)
#+end_src

#+RESULTS:
: Parameters :	 [4.83379437e-06 1.54227807e-02 3.15061312e+02]
: Parameters covariance :	 [[ 2.85416803e-15 -9.05056674e-12  4.78171606e-09]
:  [-9.05056674e-12  3.06138429e-08 -1.81982548e-05]
:  [ 4.78171606e-09 -1.81982548e-05  1.44289419e-02]]

Finalement on compare le modèle obtenu avec les données des mesures :

#+begin_src python :results file :session :var matplot_lib_filename="modelBigOsillations.png" :exports both
  co2Model = slowContributionModel(t,params[0], params[1], params[2])

  plt.figure(figsize=(10,5))
  plt.plot(t,co2,label="Measures")
  plt.plot(t,co2Model,label="Model")
  plt.legend()
  plt.title("Modélisation de la variation lente de la concentration de CO2 en fonction du temps")
  plt.xlabel("Temps [années]")
  plt.ylabel("Concentration de CO2 [ppm]")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:modelBigOsillations.png]]

On constate que l'on obtient des résultats satisfaisants pour la
caractérisation de la contribution lente.

** Extrapolation de la concentration de CO2 jusqu'à 2025

Dans un premier temps, on crée une liste qui contient l'ensemble des
données et l'on rajoute des dates jusqu'à la première semaine de 2025.

#+begin_src python :results value :session  :exports both
  initialDate = dates[0]
  finalDate   = datetime.datetime(year=2025,month=1,day=1)
  deltaWeeks = finalDate -initialDate 
  deltaWeeks = int(deltaWeeks.days/7) + 2

  newDates = [initialDate + datetime.timedelta(weeks=i) for i in range(deltaWeeks)]
  newDates[-1]
#+end_src

#+RESULTS:
: 2025-01-04 00:00:00

Finalement, on extrapole à l'aide du modèle puis on affiche les
résultats.

#+begin_src python :results file :session :var matplot_lib_filename="co2Prediction.png" :exports both
  timePrediction = np.arange(len(newDates))
  co2Prediction  = slowContributionModel(timePrediction,params[0], params[1], params[2])
  newPlotDates   = pltDates.date2num(newDates)

  plt.figure(figsize=(10,5))
  plt.plot_date(newPlotDates,co2Prediction,label="Prediction")
  plt.plot_date(plotDates,concentration,label="Measurements")
  plt.legend()
  plt.title("Extrapolation de la variation lente de la concentration de CO2 en fonction du temps")
  plt.xlabel("Temps [années]")
  plt.ylabel("Concentration de CO2 [ppm]")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:co2Prediction.png]]

On finit par calculer l'augmentation de concentration en pourcentage
par rapport à notre dernière mesure :

#+begin_src python :results output :session  :exports both
  augmentation = (co2Prediction[-1] - concentration[-1])/ concentration[-1] *100
  print("Augmentation de la concentration en % : ", augmentation)
#+end_src

#+RESULTS:
: Augmentation de la concentration en % :  2.4518838408355386


* Conclusion  

Cette étude propose un modèle de la variation de concentration de CO2,
puis utilise ce modèle afin de prédire l'augmentation de la
concentration de C02 en 2025 par rapport à la dernière mesure.

On considère que l'on peut caractériser la variation de la
concentration de C02 en deux phénomènes dont on propose un modèle
mathématique : 

1. Une variation périodique dont la période correspond à une année et
   d'amplitude 5 ppm.
2. Une variation "globale" modélisée par un polynôme d'ordre 2.