#Script illustrating how to access, plot and do simple
#data manipulations with the zonal mean text version of
#the ERBE (Earth Radiation Budget Experiment) data set
#
#
#Data on section of text which this script is associated with
Chapter = '3.**'
Figure = '**'
#
#**ToDo:
#           Add a README file stating what the data means
#
#Set the data path here for your installation
datapath = '/Users/rtp1/Havsornen/Teaching/GeoSci232/WorkbookDatasets/'
dataset = 'Chapter3Data/ERBE/'

#Imports we'll need
from ClimateUtilities import *
import phys

#Let's use the glob utility to look at what's in the data
#directory
from glob import glob
print 'Available data sets:'
for filename in glob(datapath+dataset+'*.txt'):
    shortname = filename.split('/')[-1] #Split up into directories, and
                                        #just take the last part (the name)
    print shortname

#These are ascii tables, which can be read in a word processor and
#directly read into a spreadsheet if you want. Here, we'll read it
#into a Curve object.
#
#Get the 1988 OLR data; read it into a Curve object
cOLR = readTable(datapath+dataset+'OLR1988.txt')
#
#Print out what's in the data set
print cOLR.listVariables()
#
#Looking at the output, you see you have the latitude, and columns
#representing the data for January, Feb, etc., and also the annual
#mean (Ann)
#
#As an example, let's make a new curve and plot the January OLR vs
#latitude.
c = Curve()
c.addCurve(cOLR['lat'],'latitude')
c.addCurve(cOLR['Jan'],'OLR_Jan')
#
#The following statement will switch the axes in the plot,
#if you like it better that way.
c.switchXY = True
#
#Now plot it
w1 = plot(c)
#
#Now get the 1988 surface temperature data
cT = readTable(datapath+dataset+'TsBar1988.txt')
#
#Next, compute the January surface upward radiation. Note that we
#can do an operation on a whole column at a time!
surfRad = phys.sigma*(cT['Jan'])**4
#
#Note that the result is not just a single value. It is an
#array of values -- specifically a Numeric array, on which you
#can do arithmetic, rather than a regular Python list.  The
#first three values are surfRad[0],surfRad[1],surfRad[2],
#and so forth.  You could compute the fourth root of the list,
#for example, by typing the command: surfRad**.25 
#
#Now add this to our curve and plot again
c.addCurve(surfRad,'SurfaceRad')
w2 = plot(c)
#Voila! Now we have the plot the way it looks in the book