Code that creates quicklooks

creates aoss tower quicklooks

aosstower/tower_quicklooks/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)
+                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)
+            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()    

aosstower/tower_quicklooks/test_cases.txt

+time python create_quicklook.py /Users/kgao/Desktop/ncFiles/05/aoss_tower.2013-05-01.day.nc /Users/kgao/Desktop/ncFiles/05/aoss_tower.2013-05-02.day.nc -o test --varnames air_temp dewpoint rh pressure wind_speed wind_dir accum_precip solar_flux --daily
+
+time python create_quicklook.py -s 2013-05-01 /Users/kgao/Desktop/ncFiles/05/aoss_tower.2013-05-01.day.nc -o test --varnames air_temp dewpoint rh pressure wind_speed wind_dir accum_precip solar_flux
+
+time python create_quicklook.py -s 2013-05-01 /Users/kgao/Desktop/ncFiles/05/aoss_tower.2013-05-01.day.nc -o test --varnames air_temp dewpoint
+
+time python create_quicklook.py -s 2013-05-01 /Users/kgao/Desktop/ncFiles/05/aoss_tower.2013-05-01.day.nc -o test --varnames pressure