'''
This script reads tabular data of Earth and planetary
soundings, and makes plots of temperature and potential
temperature.  It illustrates how to read in and plot
tabular data, and set the axis options for a vertical
profile plot.  This script plots cepex soundings, but
other soundings are plotted similarly, just by changing
the dataset (and sometimes the variable names).
'''

#Data on section of text which this script is associated with
Chapter = '2.**'
Figure = '**'

#Set the data path here for your installation
datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/'
dataset = 'Chapter2Data/cepex_sondes/'

from ClimateUtilities import *
from glob import glob #Imported to illustrate use of glob to find files


#Read in the data
filename = '93032700.txt'
c = readTable(datapath+dataset+'93032700.txt')
#
#If you don't know the variable names, you would
#execute the following to see if you need to extract data
#print c.listVariables()

#CEPEX sounding files contain relative humidity (rh)
#as well as p and T. If we just want to plot T,
#we want to subset
c1 = c.extract(['p','T'])
#
#Now set up the axes the way we want for a sounding
c1.switchXY = True #Interchange axes so pressure axis is verticle
c1.reverseY = True #Reverse pressure axis so pressure decreases upward
c1.YlogAxis = True #Use logarithmic pressure axis
c1.Xlabel = 'Pressure' #Put on a y axis label
c1.Ylabel = 'Temperature' #Put on an x axis label
c1.PlotTitle = 'Sounding '+filename #Plot title 
plot(c1)

#Now we'll illustrate how to build a filename from the
#date data, and use that to get the data for the plot
year = 93
month = 3
day = 20
hour = 0 #0 or 12 GMT
#Build the filename
#The format '%02d' says to format as a two-digit integer
#padded at the front by zero. E.g. '%02d'%7 is the string '07'
filename = '%2d%02d%02d%02d.txt'%(year,month,day,hour)
c = readTable(datapath+dataset+filename)
c1 = c.extract(['p','T'])
c1.switchXY=c1.reverseY = c1.YlogAxis = True
c1.PlotTitle = filename
plot(c1)


#Now we'll illustrate how to use glob to select a collection
#of files and plot them.
#This selects all March data for March 20-29
files = glob(datapath+dataset+'93032*.txt')
#
#Now get all the data and put it in a single curve
#for plotting. This example assumes all the data is on
#the same set of pressure levels
c1 = Curve()
start = True
for f in files:
    c = readTable(f)
    if start: #Put in the pressure grid
        c1.addCurve(c['p'],'p')
        start = False
    fname = f.split('/')[-1] #Extract file name without directory
    c1.addCurve(c['T'],fname)
c1.switchXY=c1.reverseY = c1.YlogAxis = True
plot(c1)