[[!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? =--  ### 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>