Experiments in carbon cycle with Sage in The Azimuth Project

[[!redirects Experiments in Carbon cycle with Sage]]
# Contents
* cc
{:toc}

## Idea 

Create some plots related to the [[Carbon cycle]] and do some simple time analysis in both time and frequency domain using [[Sage]] and its built in support for [[time series]].

## Details

### Diagrams

#### Over time

The classic Keeling curve, plotted up to January 2011 in Sage with default aspect ratio for Sage List_plot:

<div align = "center">
<img src = "http://www.azimuthproject.org/azimuth/files/mloaco2.png">
</div>

Here is a plot of the whole series but showing season corrected values (also using the default aspect ratio):

<div align = "center">
<img src = "http://www.azimuthproject.org/azimuth/files/mloaco2trend.png">
</div>

I think it can be better with setting it default aspect ratio to one, so we'll try that as well and add some major grid lines:

<div align = "center">
<img src = "http://www.azimuthproject.org/azimuth/files/mloaco2aspect1.png">
</div>

Trend - season corrected also with SAGE aspect ratio one:

<div align = "center">
<img src = "http://www.azimuthproject.org/azimuth/files/mloaco2trendaspect1.png">
</div>

Here is one with the code below showing the interpolated date moving average
over six months and also axis labels:

<div align = "center">
<img src = "http://www.azimuthproject.org/azimuth/files/mloaco26mavg.png">
</div>

Let us look at some typical properties of time series, the [[Hurst exponent]] and auto-correlation to detect degree of randomness and tendencies. The Hurst exponent is 0.23 and it indicates that the $CO_2$ time series has a tendency to decrease between values will probably be followed by an increase. We can conclude that the time series from Mauna Loa is not a random walk as well, due to the fact that random walk processes always has a Hurst exponent of 0.5.

#### Frequency Domain

Here is a draft plot of the [[FFT]] of the time series - i intend to show next how you can achieve smoothing by using hamming windows which by using Sage and the inluded libraries [Scipy/Numpy](http://www.scipy.org) which have support for that. One can also use R, but I am not very proficient in R, so I leave that to others! 

+-- {: .standout}
Question. What are those narrow peaks you see after the series has gone on for a short while?
=--

![sage timeseriesfft](http://www.sagenb.org/home/pub/2683/cells/18/sage0.png?3)

### Part of the draft Code

The code below reflects the plotting with aspect ratio equal to one and also using Timeseries for months larger than two. I am just checking that it could also be used for one

<pre class="prettyprint">
# based on http://wiki.sagemath.org/interact/web made by Marshall Hampton CC 3.0
# Staffan Liljegren removed regression code - which you can see at the above URL - and
# added the Time series features. TBD enable the fft and smoothing as part of the viewer

import urllib2 as U
import time

current_year = time.localtime().tm_year
co2data = U.urlopen('ftp://ftp.cmdl.noaa.gov/ccg/co2/trends/co2_mm_mlo.txt').readlines()
datalines = []

for a_line in co2data:
    if a_line.find('Creation:') != -1:
        cdate = a_line
    if a_line[0] != '#':
        temp = a_line.replace('\n','').split(' ')
        temp = [float(q) for q in temp if q != '']
        datalines.append(temp)

@interact
def azimco2viewer(ma_month = [1,2,6,12,24],trend= True, fft= False):
    if trend:
        co2_data = [q[5] for q in datalines]
    else:
        co2_data = [q[4] for q in datalines]
    yr_data = [q[2] for q in datalines]    
    if ma_month>2:
        co2ts = finance.TimeSeries(co2_data)
        co2ts = co2ts.simple_moving_average(ma_month)
        co2_data = co2ts.list()        
    lp = list_plot(zip(yr_data,co2_data), plotjoined=True, rgbcolor=(1,0,0))
    lp.axes_labels(['$Year$','$CO_2ppm$'])
    lp.show(axes = True, figsize = [7,7], aspect_ratio=1.0,gridlines=True)</pre>