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

** TODO COMMENT Vérification des dates
Nous faisons encore une vérification: nos dates doivent être séparées d'exactement une semaine. 
#+BEGIN_SRC python :results output :session
#dates = [date for date, _ in converted_data]
#for date1, date2 in zip(dates[:-1], dates[1:]):
#    if date2-date1 != datetime.timedelta(weeks=1):
#        print(f"Il y a {date2-date1} entre {date1} et {date2}")
#+END_SRC

#+RESULTS:
* 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.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 oscillation périodique "rapide" 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.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

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

Sur cette figure on met plus en évidence les oscillations rapides 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 ans. 
* 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.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.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.plot(freq[startIndx:fftzoomIndx],small[startIndx:fftzoomIndx])
  plt.tight_layout()

  plt.savefig(matplot_lib_filename)
  matplot_lib_filename
#+end_src

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

** Reconstruction du signal :  petites oscillations 

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