#+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 = 150; 
  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. 
On possède des données toutes les semaines. 

#+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
  co2 = co2 - np.mean(co2)
  # Concentration FFT
  concentrationFFT = np.fft.fft(co2)
  # 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 que les "petites" oscillations ont une amplitude
d'environ 1400. On tentera de filtrer donc le signal au dessus de ce
seuil.
** Filtrage des grandes oscillations 

On élimine toutes les valeurs au dessus du seuil établit précédemment.

#+begin_src python :results file :session :var matplot_lib_filename="smallOsillations.png" :exports both
  # Extraction of the small oscillations
  smallAmplitude = 1400
  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]]

** Filtrage des petites oscillations 

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 = 1400
  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]]