#---------------------------------------------------------------
#Utility routine to process radiosond datasets downloaded from
#the Integrated Gloal Radiosonde Archive dataset.  The datasets
#included in this directory are subsets of the full station datasets,
#but this script can be used to process full station datasets
#you download from the IGRA web site,
#
#The script reads through a station file, extracts the profiles
#that meet the criteria specified, and writes each profile out
#to a separate file. The station datasets contain twice-daily
#data so there is a lot there to look at.
#---------------------------------------------------------------


from ClimateUtilities import *
from glob import glob #Useful for making a list of station files

#Define functions

#Function to extract header data
#Input argument is the text string of the header record
#Output is an object containing station and date data
def extractHeaderData(HeaderString):
    station = int(HeaderString[1:6])
    year = int(HeaderString[6:10])
    month = int(HeaderString[10:12])
    day = int(HeaderString[12:14])
    hour = int(HeaderString[14:16])
    nlevs = int(HeaderString[20:24])
    header = Dummy() #Empty object to return data in
    header.station = station
    header.year = year
    header.month = month
    header.day = day
    header.hour = hour
    header.nlevs = nlevs
    return header

#Make an output file name from the header information
def outname(header):
    return '%d.%d.%d.%d.%d.txt'%\
           (header.station,header.year,header.month,header.day,header.hour)

#Main script starts hear. Loop through file, look for
#header records, and extract data if it meets the criteria

#Loop over all the data files in the directory
fnames = glob('*.dat') #List of station files
for fname in fnames:
    infile = open(fname)
    lines = infile.readlines()
    p = [] #Just to make the first record work
    for line in lines:
        if line[0] == "#": #This is a header record
            #Write out the old data list if it meets the
            #criterion.
            if len(p) > 0:            
            #Edit the following line according to what data you want
                if (header.year == 2006) & (header.month in [1,7]):
                    print outname(header)
                    #Put the data in a Curve object and write it out
                    c = Curve()
                    c.addCurve(p,'p')
                    c.addCurve(T,'T')
                    c.addCurve(dTdew,'DewPointDepression')
                    c.dump(outname(header))
                    
            #Start a new data list
            p = []
            T = []
            dTdew = []
            #
            header = extractHeaderData(line)
            
        else:
            #This is not a header record. Append data to the list
            #This extracts just pressure, temperature and dew
            #point depression (a measure of moisture content). There
            #is more data in the file, including geopotential height
            #and wind, but we don't use it.
            TT = float(line[15:20])
            TTdew = float(line[21:26])
            #The following lines only enter a data record
            #if the temperature data is available. In the
            #rare case in which 
            if TT > -900.:
                p.append(float(line[2:8])) #Pressure (Pa)
                T.append(TT/10.) #data is good 
                if (TTdew > -900.):
                    dTdew.append(TTdew/10.) #data is good
                else:
                    dTdew.append(-999.) #Missing data

    infile.close()