root/sodar/trunk/sodar.py

#!/usr/bin/python
"""
Module to process sodar data.
"""

import optparse
from sodar.utils import openAnything, findMissing 
from sodar import arrayData as a
import numpy as n
import pylab as p
import os

def plotSingle(filein, pathout):
    """
    Plot singe sodar raw data file.
    
    plotSingle(filein, pathout) -> rc
    
    filein - path to raw data file.
    fileout - path to plot images directory.
    rc - 0 if successful
    """
    
    if not os.path.exists(pathout):
        os.makedirs(pathout, mode=0775)
        
        try:
            rawDataHandle = openAnything.openAnything(filein)
            rawDataString = rawDataHandle.read()
            rawDataHandle.close()
        except:
            raise IOError("Failure to read raw data.")
        
        rawDataObject = a.rawData.RawData(rawDataString)
        formattedDataObject = a.formattedData.FormattedData(rawDataObject)
        arrayDataObject = a.ArrayData(formattedDataObject)
        
        beginStamp = arrayDataObject.beginStamp
        endStamp = arrayDataObject.endStamp
        numIntervals = arrayDataObject.numIntervals
        timeInterval = arrayDataObject.timeInterval
        timeInterval = (((float(timeInterval.days) * 24.0) +
                         (timeInterval.seconds / 3600.0)) * 60.0)
        minAltitude = arrayDataObject.minAltitude
        altInterval = arrayDataObject.altInterval
        
        def makePcolor(vector, title, fileout, cbLabel):
            fig = p.figure(1)
            axe = fig.add_subplot(1, 1, 1)
            pc = axe.pcolor(vector)
            
            axe.set_xlabel('Time (min)')
            axe.set_ylabel('Altitude (m)')
            axe.set_xbound(upper=numIntervals)
            
            xticks = axe.get_xticks()
            xticklabels = [str(x * timeInterval)
                           for x in xticks]
            axe.set_xticklabels(xticklabels)
            
            yticks = axe.get_yticks()
            yticklabels = [str(y * altInterval + minAltitude)
                           for y in yticks]
            axe.set_yticklabels(yticklabels)
            
            axe.set_title(title)
            cb = p.colorbar(pc)
            cb.set_label(cbLabel)
            
            fig.savefig(os.path.join(pathout, fileout))
            fig.clear()
            
            return 0
        
        uComponents = arrayDataObject.uComponents
        maskedUComponents = n.ma.masked_where(n.isnan(uComponents),
                                              uComponents)
        makePcolor(maskedUComponents.T,
                   'U Component of Wind Veloctiy\nFrom %s To %s' %
                   (beginStamp, endStamp),
                   'uComponents.png',
                   'Speed (cm/sec)')
        
        vComponents = arrayDataObject.vComponents
        maskedVComponents = n.ma.masked_where(n.isnan(vComponents),
                                              vComponents)
        makePcolor(maskedVComponents.T,
                   'V Component of Wind Veloctiy\nFrom %s To %s' %
                   (beginStamp, endStamp),
                   'vComponents.png',
                   'Speed (cm/sec)')
        
        wComponents = arrayDataObject.wComponents
        maskedWComponents = n.ma.masked_where(n.isnan(wComponents),
                                              wComponents)
        makePcolor(maskedWComponents.T,
                   'W Component of Wind Veloctiy\nFrom %s To %s' %
                   (beginStamp, endStamp),
                   'wComponents.png',
                   'Speed (cm/sec)')
        
        echoStrengths = arrayDataObject.echoStrengths
        maskedEchoStrenths = n.ma.masked_where(n.isnan(echoStrengths),
                                                       echoStrengths)
        makePcolor(maskedEchoStrenths.T,
                   'Echo Strength\nFrom %s To %s' %
                   (beginStamp, endStamp),
                   'echoStrengths.png',
                   'Strength (no units)')
        
        fig = p.figure(1)
        axe = fig.add_subplot(1, 1, 1)
        qv = axe.quiver(maskedUComponents,
                        maskedVComponents,
                        maskedWComponents)
        axe.set_xlabel('Time (min)')
        axe.set_ylabel('Altitude (m)')
        axe.set_xbound(upper=numIntervals)
        xticks = axe.get_xticks()
        xticklabels = [str(x * timeInterval)
                       for x in xticks]
        axe.set_xticklabels(xticklabels)
        yticks = axe.get_yticks()
        yticklabels = [str(y * altInterval + minAltitude)
                       for y in yticks]
        axe.set_yticklabels(yticklabels)
        axe.set_title('Wind Velocty\nFrom %s To %s' %
                      (beginStamp, endStamp))
        cb = p.colorbar(qv)
        cb.set_label('W Component Speed (cm/sec)')
        fig.savefig(os.path.join(pathout, 'quiver.png'))
        fig.clear()
        
        return 0

def plotAll(source, destination):
    """
    Plot all sodar raw data files.
    
    plotAll(source, destination) -> rc
    
    source - path to raw data directory in NCCOOS format.
    detination - path to plot images directory in NCCOOS format.
    rc - 0 if successful
    """
    
    walkList = findMissing.findMissing(source, destination)
    
    for filein, pathout in walkList:
        plotSingle(filein, pathout)
    return 0

def _main():
    """bin/python %prog [options] /path/to/raw/files/ /path/to/plots/dir/"""
    
    __description__ = 'Plot all sodar raw data files.'
    
    parser = optparse.OptionParser(usage=_main.__doc__,
                                   version='1.0',
                                   description=__description__)
    (values, args) = parser.parse_args()
    (source, destination) = tuple(args)
    plotAll(source, destination)
    return 0

if __name__ == "__main__":
    _main()