Xperiments in CO2

%auto # based on http://wiki.sagemath.org/interact/web made by Marshall Hampton CC 3.0 # and the time series stock example by William Stein. This code is available # under the same conditions. The Time Series in Sage are faster # than numpy arrays though. It is a bit overkill to make it like # this - classes for MeasureDate and GHGData- but hopefully we can continue adding features and data sources. # Staffan Liljegren removed regression code - which you can see at the above URL - and # added the Time series features. # 2011-03-10. Corrected the GHGData.plot_diff() and by using month as arg to # Timeseries.diffs(). # 2011-03-30. Fixed the x-axis for plot_diff() and merged plot_ma() into that method. # 2011-06-20. Aligned the time series plots correctly # TBD Handle 1958 as a special case as it only has 10 items # TBD. Add more data sources. I am currently adding Vostok ice core co2 # TBD Make it persistent over sessions and just load data-series, when needed. # TBD correlate with SOI or MEI import urllib2 as U import time class MeasureDate: def __init__(self, year, month, decimal_date, interpolated, seasonally_corrected): self.year=float(year); self.month=float(month); self.decimal_date=float(decimal_date); self.interpolated=float(interpolated); self.seasonally_corrected= float(seasonally_corrected) def __repr__(self): return '%i %i %4.2f %4.2f %4.2f %10d'%(self.year, self.month, self.decimal_date, self.interpolated, self.seasonally_corrected) # sublcass SageObject so we can save over sessions class GHGData(SageObject): def __init__(self, chem_formula): self.chem_formula = chem_formula.upper() def __repr__(self): return "%s from NOOA %s %s"%(self.chem_formula, self.nooa_col, self.cdate) def load(self,data_url='ftp://ftp.cmdl.noaa.gov/ccg/co2/trends/co2_mm_mlo.txt'): co2data = U.urlopen(data_url).readlines() self.datalines = [] # get all monthly measurements for a_line in co2data: if a_line.find('Creation:') != -1: self.cdate = a_line if a_line[0] != '#': temp = a_line.replace('\n','').split(' ') temp = [float(q) for q in temp if q != ''] self.datalines.append(temp) #if len(self.datalines) !=0: # self.save(DATA,'co2') return self.datalines def plot_data(self): d = self.co2_data if len(d) == 0: return text('end < start or no data about %s'%self.chem_formula , (0,3)) lg = '%s' % self lg = lg.split() lg = '%s %s %s'%(lg[2],lg[7],lg[-1]) lp = list_plot(d, plotjoined=True, rgbcolor='darkred', legend_label=lg,gridlines=True) lp.axes_labels(['$Year$','$CO_2 ppm$']) return lp def plot_diff(self,months,start,end): if end-start < months/12.0: return text('insufficient data',(0,3)) tsl = '%s increases over %s months' % (self.chem_formula,months) ts = self.make_timeseries().diffs(months) tma = ts.simple_moving_average(months) t = [i+(start) + (months+1)/12.0 for i in srange(0,end,months/12.0)] tsp = list_plot(zip(t,ts.list()),plotjoined=False,rgbcolor="red",legend_label=tsl, gridlines=True) tmap = list_plot(zip(t,tma.list()),ymin=0.0, ymax = 4.0,plotjoined=True,rgbcolor="orange",legend_label="Moving average") tsp.axes_labels(['$Year$','$CO_2 ppm$']) return tsp + tmap def make_timeseries(self): d = self.co2_data if len(d) == 0: return text('end < start or no data about %s'%self.chem_formula , (0,3)) tl = [d[i][1] for i in range(0,len(d)-1)] ts = finance.TimeSeries(tl) return ts chem_formulas= ['co2', 'ch4', 'n2o','o2','o3','h2o'] ghgs = dict([(s,GHGData(s)) for s in chem_formulas]) current_year = time.localtime().tm_year @interact def co2_viewer(chem_formula = chem_formulas, months = [1,3,6,12,24], seasonally_corrected = False, start = (1958,current_year,1), end = (1959,current_year,1)): ghg = GHGData(chem_formula) if seasonally_corrected: ghg.co2_data = [[q[2],q[5]] for q in ghg.load() if start <= q[2] <= end] ghg.nooa_col = " seasonally corrected" else: ghg.co2_data = [[q[2],q[4]] for q in ghg.load() if start <= q[2] <= end] ghg.nooa_col = " " pd = ghg.plot_data() pdiff = ghg.plot_diff(months,start,end) ga = graphics_array([pd,pdiff]) pd.show(figsize = [9,3]) pdiff.show(figsize = [9,3])