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

* Notes 
  - Lien de téléchargement : [[https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/weekly/weekly_in_situ_co2_mlo.csv][link]]
  - Date du téléchargement : 2020-11-16
  - The CO2 concentration is in [micro-mol CO2 per more] (ppm).
  - The weekly values have been adjusted to 12:00 hours at middle day
    of each weekly period.
* Tasks
** DONE Faire une FFT pour essayer de trouver un modèle pour l'évolution lente (filter contribution périodique) 
** TODO Rajouter des titres et légendes aux figures. 
* Lecture des données brutes
** Lecture 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 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. 
** 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 supperposition 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édément.

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

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

** FFT zoom

On effectue un zoom sur la partie positive du graphique pour tenter de
caractériser les différé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 amlitude 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 notre analyse fréquentielle.
** Filtrage de la contribution lente

On élimine toutes les valeurs au dessus des seuis établis
précédemment.

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

#+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.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

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


#+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.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

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

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

On élimine toutes les valeurs au dessus du 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.plot(freq[startIndx:fftzoomIndx],big[startIndx:fftzoomIndx])
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

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

** Reconstruction du signal : grandes oscillations 

#+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.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:bigOsillationsTime.png]]
** Comparaison des signaux bruts et filtrés

#+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.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

#+RESULTS:
[[file:bigOsillationsComparisons.png]]
* 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.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.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