create_quicklook.py

import os
from datetime import datetime as dt
from datetime import timedelta as delta
import logging
import pandas as pd
from netCDF4 import MFDataset, MFTime
import numpy as np
import matplotlib.pyplot as plt
import math

CHOICES = ['air_temp', 'dewpoint',
   'rh', 'pressure', 'wind_speed', 'wind_dir', 'accum_precip', 'solar_flux']

TITLES = {'air_temp': 'Temperature and Dewpoint(°C)',
         'rh': 'Relative Humidity(%)', 'pressure': 'Pressure(hpa)', 
         'wind_speed': 'Wind Speed(m/s)', 'wind_dir': 'Wind Direction(°)',
         'accum_precip': 'Accumulated Precipitation Since 00Z (mm)', 
         'solar_flux': 'Solar Flux(W/m^2)'}

IND_TITLES = {'air_temp': 'Temperature and Dewpoint',
         'rh': 'Relative Humidity', 'pressure': 'Pressure', 
         'wind_speed': 'Wind Speed', 'wind_dir': 'Wind Direction',
         'accum_precip': 'Accumulated Precipitation Since 00Z', 
         'solar_flux': 'Solar Flux'}             

# The purpose of this method is to take a string in the format
# YYYY-MM-DDTHH:MM:SS and convert that to a datetime object
# or a string in the format YY-MM-DD and convert that to a datetime object
# used in coordination with argparse -s and -e params
# @param datetime string
# @return datetime object

def _dt_convert(datetime_str):
    #parse datetime string, return datetime object
    try:
        return dt.strptime(datetime_str, '%Y-%m-%dT%H:%M:%S')
   
    except:
        return dt.strptime(datetime_str, '%Y-%m-%d')

# The purpose of this method is to take a list of level_b1 filepaths
# and convert those into a pandas frame with only good data
# @param input_files - list of level_b1 filepaths
# @return frame - pandas dataframe with only good data

def get_data(input_files):
    files = MFDataset(input_files)

    #dictionary with var_name: good data pairs
    dataDict = {}

    #get the data from the files
    for name in CHOICES:
        data = files.variables[name][:]

        qc_data = files.variables['qc_' + name][:]
        
        dataDict[name] = data
        dataDict['qc_' + name] = qc_data

    #for some reason, time has floating point precision issues
    #so I'm using base_time and offsets instead
    #might need to fix later
    base_time = files.variables['base_time'][:]

    #convert base_time epoch format into date_time object
    base_time_obj = dt(1970,1,1) + delta(seconds=int(base_time))
    
    #get offsets from files
    off_sets = files.variables['time_offset']

    #convert that into offsets from the first file's base_time
    off_sets = MFTime(off_sets)[:]  

    #for each offset, convert that into a datetime object
    stamps = [base_time_obj + delta(seconds=int(s)) for s in off_sets]

    #append all stamps into the data frame
    dataDict['stamps'] = stamps

    return pd.DataFrame(dataDict).set_index(['stamps'])

# The purpose of this method is to use the qc mask
# to only get the good data
# @param qc_data - np array with nan for good values and various set bits for others
# @param stamps - all stamps corresponding to data
# @param data - data
# @return good data and their corresponding stamps   

def get_good_data(qc_data, stamps, data):
    #all good data and stamps associated with them
    good_data = []
    good_stamps = []

    #get good indices
    good_indices = np.where(np.isnan(qc_data))[0]

    for idx in good_indices:
        #get good_data
        good_data.append(data[idx])
        good_stamps.append(stamps[idx])

    return [good_data, good_stamps]

# The purpose of this method is to determines all the 12:00:00 days
# within a start and end date
# @param cur_dt - start time
# @param dates - end time
# @return dates in lists

def get_dates_in_range(cur_dt, end):
    curr = cur_dt
    while curr <= end:
        yield curr
        curr += delta(days=1)

# The purpose of this method is to get the min and max of 
# an array of stamps and determines all the 12:00:00 days
# in that array
# @param dewpoint stamps - time stamps for valid dewpoint data. Only specified if
# air temp and dew point need to be ploted in the same plot
# @param dates - time stamps for data
# @return datesInRange - all days with 12:00:00
def find_half_days(dewpoint_stamps, dates):
    date_max = np.amax(dates)
    date_min = np.amin(dates)

    if(dewpoint_stamps):
        dew_max = np.amax(dewpoint_stamps)
        dew_min = np.amin(dewpoint_stamps)

        date_max = max(date_max, dew_max)
        date_min = min(date_min, dew_min)    

    cur_dt = date_min.replace(hour=12, minute=0, second=0)
    
    if cur_dt < date_min:
        cur_dt += delta(days=1)

    datesInRange = list(get_dates_in_range(cur_dt, date_max))    

    return datesInRange           

# The purpose of this method is to create a thumbnail plot
# @param frame - all data I need in a data frame
# @param ymin - lower limit of the y axis if specified. set only when --daily is given
# @param ymax - upper limit of yaxis if specified. Set only when --daily is given
# @param dewpoint_stamps - timestamps for valid dewpoint data. Only specified if 
# air temp and dew point need to be plotted in the same plot
# @param dewpoint_data - valid dewpoint data. Only specified if 
# air temp and dew point need to be plotted in the same plot
#@param output - output filename pattern
def thumbnail_plot(dates, data, ymin, ymax,
    dewpoint_stamps, dewpoint_data, output, wind_dir):

    if(dewpoint_stamps):
        plt.plot(dates, data, 'r', dewpoint_stamps, dewpoint_data, 'b')
        plt.axis([min(dates), max(dates), ymin, ymax])
        #plt.plot(dewpoint_stamps, dewpoint_data, 'b')

    elif(wind_dir):
        plt.plot(dates, data, 'ko', markersize=0.1)
        plt.axis([min(dates), max(dates), 0, 360])

    else:
        plt.plot(dates, data, 'k')
        plt.axis([min(dates), max(dates), ymin, ymax])

    #scale
    fig = plt.gcf()
    dpi = fig.get_dpi()
    fig.set_size_inches(80/float(dpi), 80/float(dpi))

    #create axes
    axes = plt.gca()
    axes.axes.get_xaxis().set_ticklabels([])
    axes.axes.get_xaxis().set_tick_params(length=7)
    axes.axes.get_yaxis().set_ticklabels([])
    axes.axes.get_yaxis().set_visible(False)
    axes.spines['right'].set_visible(False)
    axes.spines['bottom'].set_visible(False)

    #got the basic shape of the plot
    #still need to draw a line in the middle
    datesInRange = find_half_days(dewpoint_stamps, dates)
    
    for date in datesInRange:
        plt.axvline(x=date, color='k')       

# The purpose of this method is to create a thumbnail
# @param frame - all data I need in a data frame
# @param ymin - lower limit of the y axis if specified. set only when --daily is given
# @param ymax - upper limit of yaxis if specified. Set only when --daily is given
# @param create_air_dew_plot - boolean specifying if air temp and dewpoint
# @param output - output filename pattern 
# should be in same plot
def create_thumbnail(frame, ymin, ymax, create_air_dew_plot, output):
    #see if we're going to need subplots
    #subplots needed when len is greater than 2 for one variable
    #or greater than 4 for two variables

    need_subplots = len(list(frame.columns.values)) > 2 and not create_air_dew_plot
    need_subplots = need_subplots or len(list(frame.columns.values)) > 4 and create_air_dew_plot

    if need_subplots:
        plots_created = 1

        if create_air_dew_plot:
            numberPlots = (len(list(frame.columns.values)) - 2) / 2
            plotNumber = numberPlots * 100 + 10

        else:
            numberPlots = len(list(frame.columns.values)) / 2
            plotNumber = numberPlots * 100 + 10    

    stamps = list(frame.index)

    #get dewpoint data
    if create_air_dew_plot:
        all_data = frame['dewpoint']
        qc_data = frame['qc_dewpoint']

        good_list = get_good_data(qc_data, stamps, all_data)
        dewpoint_data = good_list[0]
        dewpoint_stamps= good_list[1]    

    for name in CHOICES:
        if name == 'dewpoint' and create_air_dew_plot:
            continue

        if name not in list(frame.columns.values):
            continue        

        elif need_subplots:
            plt.subplot(plotNumber + plots_created)
            plots_created += 1


        all_data = frame[name]
        qc_data = frame['qc_' + name]

        good_list = get_good_data(qc_data, stamps, all_data)
        good_data = good_list[0]
        good_stamps= good_list[1]

        #if we don't need two lines in same plot, plot
        if not (create_air_dew_plot and name == 'air_temp'):
            if name != 'wind_dir':    
                thumbnail_plot(good_stamps, good_data, ymin[name], ymax[name], 
                    None, None, output, False)

            if name == 'wind_dir':
               thumbnail_plot(good_stamps, good_data, ymin[name], ymax[name], 
                    None, None, output, True)      

        else:        
            thumbnail_plot(good_stamps, good_data, ymin[name], ymax[name], 
                dewpoint_stamps, dewpoint_data, output, False)

    plt.savefig(output + '.thumbnail.png')

def full_plot_stamps(stamps):
    first_stamp = np.sort(stamps)[0]

    first_stamp = first_stamp.replace(hour=0, minute=0, second=0)

    #print([((s - first_stamp).total_seconds()/3600) for s  in stamps])

    #print([(int((s - first_stamp).total_seconds()) / 3600) for s  in stamps])

    return [((s - first_stamp).total_seconds()/3600) for s  in stamps]

# The purpose of this method is to change how the y ticks are shown
# @param min - minimum of the y range right now
# @param max - maximum of the y range right now
# @return new tick labels

def get_new_labels(mini, maxi):
    delta = math.ceil((maxi - mini) / 6)

    return np.arange(mini, (mini + delta*6), delta)                 
        
# The purpose of this method is to create a full_size plot
# @param frame - all data I need in a data frame
# @param ymin - lower limit of the y axis if specified. set only when --daily is given
# @param ymax - upper limit of yaxis if specified. Set only when --daily is given
# @param dewpoint_stamps - timestamps for valid dewpoint data. Only specified if 
# air temp and dew point need to be plotted in the same plot
# @param dewpoint_data - valid dewpoint data. Only specified if 
# air temp and dew point need to be plotted in the same plot
# @param output - output filename pattern
# @param wind_direction - says different yaxis labels
# @param accum_precip - different yaxis handling
# @param see if we need subplots
def full_plot(dates, data, ymin, ymax,
    dewpoint_stamps, dewpoint_data, output, wind_dir,
    accum_precip, need_subplots):

    dates = full_plot_stamps(dates)

    if(dewpoint_stamps):
        dewpoint_stamps = full_plot_stamps(dewpoint_stamps)

        plt.plot(dates, data, 'r', dewpoint_stamps, dewpoint_data, 'b')
        plt.axis([min(dates), max(dates), ymin, ymax])
        

        #labels = [math.ceil(float(item.get_text())) for item in ticks]
        #plt.plot(dewpoint_stamps, dewpoint_data, 'b')


    elif(wind_dir):
        plt.plot(dates, data, 'ko', markersize=1)
        plt.axis([min(dates), max(dates), 0, 360])

        axes = plt.gca()
        axes.get_yaxis().set_ticks([0,90,180,270])


    else:
        plt.plot(dates, data, 'k')
        plt.axis([min(dates), max(dates), ymin, ymax])

    if not need_subplots:
        return

    #get axes
    axes = plt.gca()

    #if not wind_direction, reset yaxis ticks
    if not wind_dir:
        if not accum_precip:
            mini, maxi = axes.get_ylim()
            new_labels = get_new_labels(mini, maxi)

        else:
            mini, maxi = axes.get_ylim()
            delta = ((maxi - mini)/6)
            new_labels = np.arange(mini, (mini + delta*6), delta)
            
            if len(new_labels) == 0:
                new_labels = [0, 0.05, 0.1]   

        plt.yticks(new_labels)
        axes.get_yaxis().get_major_ticks()[-1].set_visible(False)

    #hid last tick because for subplots above the last one, it's peeping out...
    # will show for the last subplot or when subplots are not needed
    # in other words, this code should only affect plots above the last plot
    # only when sublots are needed
    axes.get_xaxis().get_major_ticks()[-1].set_visible(False)
    axes.get_xaxis().get_major_ticks()[0].set_visible(False)

# The purpose of this method is to get the subtitle's y coordinate
# based upon how many subplots there are
# @param numSubPlots - number of subplots
# @return y-coordinate of subtitle
def get_subtitle_location(numSubPlots):
    return 1 - 0.05*numSubPlots

# The purpose of this method is to create a full size quicklook
# @param frame - all data I need in a data frame
# @param ymin - lower limit of the y axis if specified. set only when --daily is given
# @param ymax - upper limit of yaxis if specified. Set only when --daily is given
# @param create_air_dew_plot - boolean specifying if air temp and dewpoint 
# should be in same plot
def create_full_plot(frame, ymin, ymax, create_air_dew_plot, output):
    #see if we're going to need subplots
    #subplots needed when len is greater than 2 for one variable
    #or greater than 4 for two variables

    need_subplots = len(list(frame.columns.values)) > 2 and not create_air_dew_plot
    need_subplots = need_subplots or len(list(frame.columns.values)) > 4 and create_air_dew_plot

    #get need supplot vars
    if need_subplots:
        plots_created = 1

        if create_air_dew_plot:
            numberPlots = (len(list(frame.columns.values)) - 2) / 2
            plotNumber = numberPlots * 100 + 10

        else:
            numberPlots = len(list(frame.columns.values)) / 2
            plotNumber = numberPlots * 100 + 10    

    stamps = list(frame.index)

    #get dewpoint data
    if create_air_dew_plot:
        all_data = frame['dewpoint']
        qc_data = frame['qc_dewpoint']

        good_list = get_good_data(qc_data, stamps, all_data)
        dewpoint_data = good_list[0]
        dewpoint_stamps= good_list[1]       

    for name in CHOICES:
        #if dewpoint with temp or a name not in frame
        # continue
        if name == 'dewpoint' and create_air_dew_plot:
            continue

        if name not in list(frame.columns.values):
            continue 

        all_data = frame[name]
        qc_data = frame['qc_' + name]

        good_list = get_good_data(qc_data, stamps, all_data)
        good_data = good_list[0]
        good_stamps= good_list[1]

        dateString = good_stamps[0].strftime('%m/%d/%Y')

        # create plots
        if need_subplots:
            if plots_created == 1:
                plt.figure().suptitle('AO&SS Building Tower Meteorogram ' + dateString, fontsize=13)
                ax1 = plt.subplot(plotNumber + plots_created)
                ax1.set_title(TITLES[name], x=0.5, y=get_subtitle_location(numberPlots), fontsize=8)


            else:
                curr_plot = plt.subplot(plotNumber + plots_created, sharex=ax1)
                curr_plot.set_title(TITLES[name], x=0.5, y=get_subtitle_location(numberPlots), fontsize=8)

            plots_created  += 1

        else:
            plt.figure().suptitle('AO&SS Building Tower ' + IND_TITLES[name] + ' ' + dateString, fontsize=13)
            plt.ylabel(TITLES[name])

        #if we don't need two lines in same plot, plot
        if not (create_air_dew_plot and name == 'air_temp'):
            full_plot(good_stamps, good_data, ymin[name], ymax[name], 
                        None, None, output, name == 'wind_dir',
                        name == 'accum_precip', need_subplots)      

        else:        
            full_plot(good_stamps, good_data, ymin[name], ymax[name], 
                dewpoint_stamps, dewpoint_data, output, False, False,
                need_subplots)                  

    if need_subplots:
        plt.subplots_adjust(hspace=0, bottom=0.125)

    plt.xlabel('Time (UTC)')
 
    #reverse code done in create_full_plot
    axes = plt.gca()
    axes.get_xaxis().get_major_ticks()[-1].set_visible(True)            
    axes.get_xaxis().get_major_ticks()[0].set_visible(True)            
    
    plt.savefig(output + '.png')         

def main():
    import argparse

    #argparse description
    parser = argparse.ArgumentParser(description="Use data from level_b1 netCDF files to create netCDF files")

    #argparse verbosity info
    parser.add_argument('-v', '--verbose', action='count', 
                       default=int(os.environ.get("VERBOSITY", 2)), 
                       dest='verbosity',
                       help=('each occurence increases verbosity 1 level through'
                             + ' ERROR-WARNING-INFO-DEBUG (default INFO)'))

    #argparse start and end times
    parser.add_argument('-s', '--start-time', type=_dt_convert,
    help="Start time of plot. If only -s is given, a plot of " +
        "only that day is created. Formats allowed: \'YYYY-MM-DDTHH:MM:SS\', \'YYYY-MM-DD\'")

    parser.add_argument('-e', '--end-time', type=_dt_convert,
    help="End time of plot. If only -e is given, a plot of only that day is " +
          "created. Formats allowed: \'YYYY-MM-DDTHH:MM:SS\', \'YYYY-MM-DD\'")

    #netcdf file paths
    parser.add_argument("input_files", nargs="+", help="aoss_tower_level_b1 paths")

    #output filename pattern
    parser.add_argument('-o', '--output', help="filename pattern")

    #thumbnail or full
    parser.add_argument('-t', '--thumb-nail', action='store_true', help="if specified, script creates a thumbnail")

    #plot names
    parser.add_argument('--varnames', nargs="+", choices=CHOICES,
    required=True,
    help="the variable names for the desired plot. Valid choices: air_temp, " +
    "dewpoint, rh, pressure, wind_speed, wind_dir, accum_precip, solar_flux")

    #--daily flag
    parser.add_argument('-d', '--daily', action='store_true', 
    help="creates a plot for every day. Usually used to create plots " +
    "that will line up for aoss tower quicklooks page")

    args = parser.parse_args()

    levels = [logging.ERROR, logging.WARN, logging.INFO, logging.DEBUG]

    level = levels[min(3, args.verbosity)]
   
    logging.basicConfig(level=level)

    var_names = args.varnames

    create_air_dew_plot = False

    #see if we need to plot air_temp and dewpoint on the same plot
    if 'air_temp' in var_names and 'dewpoint' in var_names:
        create_air_dew_plot = True
        var_names.remove('dewpoint')

    #get data
    frame = get_data(args.input_files)
    
    #only have the data we need
    for name in CHOICES:
        if name == 'dewpoint' and create_air_dew_plot:
            continue

        if name not in var_names:
           del frame[name]
           del frame['qc_' + name]

    #frame only contains data from start-end times
    if(args.start_time and args.end_time):
        start_filter = args.start_time.strftime('%Y-%m-%d %H:%M:%S')
        end_filter = args.end_time.strftime('%Y-%m-%d %H:%M:%S')
        frame = frame[start_filter: end_filter]

    #frame only contains data from start-end of that day
    elif(args.start_time):
        end_time = args.start_time.replace(hour=23, minute=59, second=59)
        start_filter = args.start_time.strftime('%Y-%m-%d %H:%M:%S')
        end_filter = end_time.strftime('%Y-%m-%d %H:%M:%S')
        frame = frame[start_filter: end_filter]

    ymin = {}
    ymax = {}    

    if not args.daily:
        for name in CHOICES:
            if name == 'dewpoint' and create_air_dew_plot:
                continue

            if name not in frame:
                continue

            ymin[name] = None
            ymax[name] = None    

        if(args.thumb_nail):
            create_thumbnail(frame, ymin, ymax, create_air_dew_plot, args.output)

        else:
            create_full_plot(frame, ymin, ymax, create_air_dew_plot, args.output)

    else:
        for name in CHOICES:
            if name == 'dewpoint' and create_air_dew_plot:
                continue

            if name not in frame:
                continue

            if name == 'accum_precip':    
                ymin[name] = frame[name].min()
                ymax[name] = frame[name].max()

            else:
                ymin[name] = math.floor(frame[name].min())
                ymax[name] = math.ceil(frame[name].max())    

        frameList = [(group[1]) for group in frame.groupby(frame.index.day)] 

        counter = 0

        for frame in frameList:
            if args.thumb_nail:
                create_thumbnail(frame, ymin, ymax, create_air_dew_plot, args.output.format(str(counter)))
            else:
                create_full_plot(frame, ymin, ymax, create_air_dew_plot, args.output.format(str(counter)))
            counter += 1
            plt.gcf().clear()          

if __name__ == "__main__":
    main()