#+TITLE:  Concentration de CO2 dans l'atmosphère depuis 1958
#+AUTHOR: Miguel Arpa Perozo
#+DATE:   \today
#+LANGUAGE: fr
# #+PROPERTY: header-args :eval never-export
#+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 : 
  - 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."
    
** Démarche générales de l'étude
   1. Dans un premier temps effectuer une analyse fréquentielle des
      données afin de tenter de modéliser certains phénomènes.
   2. Utiliser un outil d'optimisation pour déterminer certains
      paramètres du modèle proposé en utilisant les 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 des données brutes ou téléchargement
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 premières 45 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

Visualisation des premières colonnes du tableau : 

#+BEGIN_SRC python :results value :session 
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 traiter les données plus facilement.

** Converstion 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 dates 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 [semaines]")
  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. 
2. Une croissance lente avec une forme qui ressemble à un début de
   parabole.
** /Zoom/ sur l'affichage des données 

Effectuons un /zoom/ sur une période de temps plus courte caractériser
les oscillations rapides observés 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 [semaines]")
  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 plus en évidence le phénomène périodique sur
le niveau de concentration de CO2. D'après la figure, dans un premier
temps on peut supposer que les oscillations ont une période de un an et une
amplitude entre cinq et sept.
* 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. 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 unité de temps")
  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 quelques pics qui devraient correpondre au phénomène
périodique. On tentera de caractériser plus précisément ses pics.
** FFT zoom

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 unité de temps")
  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'amplitudes 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 le graphique temporel des
données. En effet on a supposé précédemment que les oscillations
avaient une période de un an, notre unité de temps étant une semaine,
on peut déduire la fréquence des oscillation en supposant que une
année a 52 semaine : \( f = \frac{1}{52} \approx 0.019 \) ce qui reste
cohérent avec les résultats dans le graphique.
** Filtrage de la contribution lente

On élimine toutes les valeurs supérieures aux seuils établis précédemment
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 unité de temps")
  plt.ylabel("Amplitude")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

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

** Reconstruction du 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 sur 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 unité de temps")
  plt.ylabel("Amplitude")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

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

** Reconstruction du signal : contribution lente

Finalement on reconstruit le signal temporel puis l'on 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 [semaines]")
  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 nous permet pas de déduire grand chose quant à cette
contribution lente. Tentons de le comparer avec les mesures. 
** Comparaison des signaux filtré 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 nous permet pas d'avoir plus
d'informations pertinentes quant à cette contribution lente. Donc par
la suite on décidera de modéliser 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. 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 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 peut supposer dans un premier temps 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 [semaines]")
  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 [semaines]")
  plt.ylabel("Concentration de CO2 [ppm]")
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

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

Calcul de l'augmentation de concentration en pourcentage : 

#+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