root/sodar/trunk/plotSodar.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
    numAltitudes = arrayDataObject.numAltitudes
    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 (hh:mm UTC)')
        axe.set_ylabel('Altitude (m)')
        axe.set_xbound(upper=numIntervals)
        
        xticks = axe.get_xticks()
        xticklabels = [(int(x) * timeInterval) + beginStamp
                       for x in xticks]
        xticklabels = [':'.join(('%02u' % x.hour, '%02u' % x.minute))
                       for x in xticklabels]
        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 for %s' %
                (str(beginStamp)[:10],),
               '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 for %s' %
                (str(beginStamp)[:10],),
               '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 for %s' %
                (str(beginStamp)[:10],),
               'wComponents.png',
               'Speed (cm/sec)')
    
    echoStrengths = arrayDataObject.echoStrengths
    maskedEchoStrenths = n.ma.masked_where(n.isnan(echoStrengths),
                                                   echoStrengths)
    makePcolor(maskedEchoStrenths.T,
               'Echo Strength for %s' %
                (str(beginStamp)[:10],),
               'echoStrengths.png',
               'Strength (no units)')
    
    # timeComponent = n.array(range(numIntervals))
    # altComponent = n.array(range(numAltitudes))
    # fig = p.figure(1)
    # axe = fig.add_subplot(1, 1, 1)
    # qv = axe.quiver(altComponent,
    #                 timeComponent,
    #                 maskedUComponents.T,
    #                 maskedVComponents.T)
    # 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 for %s' %
    #               (str(beginStamp)[:10],))
    # 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
    """
    
    if source.endswith(os.path.sep):
        source = source.rstrip(os.path.sep)
    
    if destination.endswith(os.path.sep):
        destination = destination.rstrip(os.path.sep)
    
    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='%prog 1.0',
                                   description=__description__)
    parser.set_defaults(forceUpdate=False)
    parser.add_option('-f', '--force-update',
                      action='store_true',
                      dest='forceUpdate')
    (values, args) = parser.parse_args()
    (source, destination) = tuple(args)
    plotAll(source, destination)
    
    return 0

if __name__ == "__main__":
    _main()