Added all the QC files I've been working on

testing/correlation_test.py

+import pandas as pd
+import numpy as np
+import scipy.signal
+from aeri_tools.io.dmv import housekeeping
+from aeri_tools.io.dmv import radiance
+from sys import argv
+import os
+import time
+import matplotlib.pyplot as plt
+
+def main(path):
+
+    #get data from B1.CXS file
+    dmv = radiance.get_radiance_from_rnc(path)
+    hk = housekeeping.get_all_housekeeping(path)
+
+    data = pd.DataFrame(index=hk.index, columns=('dmv', 'hk'))
+    data['dmv'] = dmv.xs(str.encode('H'), level='scene').mean(axis=1)
+    data['hk'] = hk['SCEtemp']
+
+    data.dropna(inplace=True)
+
+    try:
+        for x in data.columns:
+            tmp = scipy.signal.savgol_filter(data[x].values[:], 201, 3)
+            data[x] = data[x] - tmp
+
+        correlation = np.correlate(data['dmv'].values[:], data['hk'].values[:], mode='same')
+
+        print('max = ', np.amax(correlation), ' : min = ', np.amin(correlation))
+
+        fig, ax = plt.subplots(1, figsize=(15,10), sharex=True)
+
+        '''
+        for x, val in ('dmv', 'hk'):
+            plt.sca(ax[x])
+            plt.plot(list(data.index), data[val].values[:])
+        '''
+
+        plt.plot(data.index, correlation)
+
+        #plt.show()
+        name = '/Users/adiebold/Documents/corr_pngs/' + path[-12:-6] + '.png'
+        plt.savefig(name)
+        plt.clf()
+
+    except:
+        print('FAIL')
+
+
+if __name__ == '__main__':
+    start_time = time.time()
+
+    pathname = argv[1]
+
+    skip_num = 0
+    curr_num = 0
+    if os.path.isdir(pathname):
+        for filename in os.listdir(pathname):
+            full_path = pathname + '/' + filename
+            if os.path.isdir(full_path):
+                for fname in os.listdir(full_path):
+                    fuller_path = full_path + '/' + fname
+                    if os.path.isfile(fuller_path) and 'B1.CXS' in fname:
+                        print(curr_num, ': ', fuller_path)
+                        if curr_num >= skip_num:
+                            main(fuller_path)
+                        curr_num += 1
+            elif os.path.isfile(full_path) and 'B1.CXS' in filename:
+                if curr_num >= skip_num:
+                    main(full_path)
+                curr_num += 1
+                print(curr_num, ': ', filename)
+    else:
+        main(pathname)
+
+    print('execution time: %d minute(s), %.2f second(s)' % ((time.time()-start_time)//60, (time.time()-start_time)%60))

testing/detector_temp_flag.py

+from __future__ import print_function
+from scipy.signal import periodogram
+import pandas as pd
+import numpy as np
+from aeri_tools.io.dmv import housekeeping
+from matplotlib import pyplot as plt
+import os
+WIN_SIZE = 40
+DETECTOR_TEMP_KEY = 'detectorTemp'
+TIME_WINDOW = (60,600) #seconds
+MIN_AMP = .02
+WIN_TYPE = 'blackmanharris'
+def sig_finder_function(rolling_window,**kwargs):
+    """Function applied to detectorTemp dataset with pd.rolling_apply
+    Calculates the dominant period and
+    """
+
+    #calculate the time difference and average frequency
+    nsamples = rolling_window.size
+    times = pd.Series(rolling_window.index)
+    win_type = kwargs.get('win_type',WIN_TYPE)
+    try:
+        #get the time span covered by the window
+        dt = (times.iloc[-1] - times.iloc[0]).seconds
+    except AttributeError:
+        #catch NaT indicies
+        return np.nan
+    else:
+        #compute the dominant frequency/period
+        fs = nsamples/dt
+        f, Pxx_den = periodogram(rolling_window,fs,window=win_type)
+        max_i = Pxx_den.argmax()
+        if f[max_i] != 0:
+            primary_T = 1./f[max_i]
+        else:
+            primary_T = float("inf")
+        return primary_T
+def series_rolling_apply(series,window,func,**kwargs):
+    """rolling applies function to series by passing subsection of series,
+    rather than np array
+    """
+    #lambda to map the passed function to a rolling window
+    map_f=lambda i:np.nan if i<window else func(series.iloc[i-window:i],**kwargs)
+    return pd.Series(map(map_f,range(series.size)),index = series.index)
+
+def flag_detector_temp_oscillations(path_or_dmv,window_size = WIN_SIZE,std=10,**kwargs):
+    """Performs power spectral density analysis on detector temperature,
+    and flags the data if a consistent period is found
+    """
+    #read in the housekeeping dataframe
+    housekeeping_df = housekeeping.get_all_housekeeping(path_or_dmv)
+    #get the detector temp variable
+    detector_temp = housekeeping_df[DETECTOR_TEMP_KEY]
+    #detrend by subtracting the day's mean
+    detrend_dtemp = detector_temp - detector_temp.mean()
+    #approximate periodicity in data using a rolling periodogram
+    periods = series_rolling_apply(detrend_dtemp,window_size,sig_finder_function,**kwargs)
+    #check for big spikes in the data (these throw off the periodogram)
+    #currently disabled
+    dtt_roller = detrend_dtemp.rolling(window_size)
+    amps = False & ((dtt_roller.max()-dtt_roller.min()) > .2)
+    #look for low standard deviation in the period (= a consistent periodicity)
+    stds = series_rolling_apply(periods,window_size,lambda x:x.std() < std)
+    #flag as bad if there's a consistent periodicity or large jump in the data
+    flags = stds | amps
+    return flags, periods
+
+def flatten_out(original_flags, pds):
+    '''incrementing windows to help flatten out the 0/1 value as either
+    0 or 1; removes small blips
+    '''
+
+    flags = original_flags.copy()
+    values_dict = flags.value_counts()
+
+    num_half_hour = int(len(flags)/48)
+    curr_flags = []
+    std_data = {}
+    std_data['curr_std'] = []
+    std_data['mean_std'] = []
+    std_data['mean'] = []
+    std_data['std_curr_std'] = []
+    std_data['std_mean_std'] = []
+
+    lower = 0
+    upper = 1
+
+    while upper <= len(flags):
+        values_dict = original_flags[lower:upper].value_counts()
+        if True in values_dict.keys():
+            mean = values_dict[True] / len(flags[lower:upper])
+        else:
+            mean = 1 - (values_dict[False] / len(flags[lower:upper]))
+        curr_std = np.std(pds[lower:upper])
+        mean_std = np.mean(pds[lower:upper])
+        std_data['mean'].append(mean)
+        std_data['curr_std'].append(curr_std)
+        std_data['mean_std'].append(mean_std)
+        std_curr_std = np.std(std_data['curr_std'][lower:upper])
+        std_data['std_curr_std'].append(std_curr_std)
+        std_mean_std = np.std(std_data['mean_std'][lower:upper])
+        std_data['std_mean_std'].append(std_mean_std)
+        if std_mean_std < 20 or np.isnan(std_mean_std):
+            flags[lower:upper] = flags[lower:upper].replace(False, True)
+        else:
+            flags[lower:upper] = flags[lower:upper].replace(True, False)
+
+        upper += 1
+        if upper - lower > num_half_hour:
+            lower += 1
+
+    new_flags = flags & (np.asarray(std_data['mean']) > 0.35)
+    # new_flags = np.asarray(std_data['mean']) > 0.8
+    # fig, ax = plt.subplots(2, figsize=(15,10),sharex=True)
+    # plt.plot(flags)
+    # plt.plot(new_flags)
+    # plt.show()
+    # plt.clf()
+    return new_flags, std_data
+
+def main(path):
+    """Test function that displays detectorTemp alongside the computed flag
+    using different windowing functions
+    """
+    dmv = housekeeping.get_all_housekeeping(path)
+    detector_temp = dmv[DETECTOR_TEMP_KEY]
+    fig,ax = plt.subplots(5,figsize=(15,10),sharex=True)
+    plt.sca(ax[0])
+    plt.plot(dmv['Time'],detector_temp,'r')
+    # plt.ylim(76.9,77.6)
+    plt.ylabel("Detector Temperature (K)")
+    plt.grid()
+
+    flags,pds = flag_detector_temp_oscillations(path,win_type=WIN_TYPE)
+    flags, std_data = flatten_out(flags.copy(), pds)
+
+    plt.sca(ax[1])
+    plt.plot(dmv['Time'],pds,lw=2,label=WIN_TYPE)
+    plt.grid(True)
+    plt.ylabel("Dominant Period (s)")
+    #plt.legend()
+
+    plt.sca(ax[2])
+    plt.plot(dmv['Time'],flags,lw=2,label=WIN_TYPE)
+    plt.grid(True)
+    plt.ylabel("Standard Deviation Flag")
+    plt.xlabel("Time (Hours UTC)")
+    plt.ylim(-0.1,1.1)
+    plt.sca(ax[3])
+    plt.plot(dmv['Time'], np.asarray(std_data['mean']))
+
+    title = os.path.split(path)[-1][:6]
+    plt.suptitle(title,fontsize=14)
+
+    plt.sca(ax[4])
+    plt.plot(dmv['Time'], std_data['curr_std'])     #blue
+    plt.plot(dmv['Time'], std_data['mean_std'])     #green
+    plt.plot(dmv['Time'], std_data['std_curr_std']) #red
+    plt.plot(dmv['Time'], std_data['std_mean_std'])
+
+    print('title: ', title)
+    plt.savefig("sgp_detect_pngs/detector_temp_flag_test_%s.png"%title,bbox_inches='tight')
+    # plt.savefig("detector_temp_flag_test_%s_2.png"%title,bbox_inches='tight')
+    # plt.show()
+    plt.clf()
+    plt.close()
+
+if __name__ == '__main__':
+    from sys import argv
+    import os
+    import time
+
+    start_time = time.time()
+
+    pathname = argv[1]
+
+    skip_num = 0
+    curr_num = 0
+    print('skip_num is ', skip_num)
+
+    if os.path.isdir(pathname):
+        for filename in os.listdir(pathname):
+            full_path = pathname + '/' + filename
+            if os.path.isdir(full_path):
+                for fname in os.listdir(full_path):
+                    fuller_path = full_path + '/' + fname
+                    if os.path.isfile(fuller_path) and 'B1.CXS' in fname:
+                        if curr_num >= skip_num:
+                            main(fuller_path)
+                        curr_num += 1
+                        print(curr_num, ': ', fuller_path)
+            elif os.path.isfile(full_path) and 'B1.CXS' in filename:
+                if curr_num >= skip_num:
+                    main(full_path)
+                curr_num += 1
+                print(curr_num, ': ', filename)
+    else:
+        main(pathname)
+
+    print('execution time: %d minute(s), %.2f second(s)' % ((time.time()-start_time)//60, (time.time()-start_time)%60))

testing/quick_vis.py

+import netCDF4
+import matplotlib.pyplot as plt
+import matplotlib.patches as patches
+import time
+from scipy import arange
+import os
+import time
+import numpy as np
+
+def run(filename):
+
+    subplot_start_time = time.time()
+
+    print('program running...')
+    plt.rcParams['figure.figsize'] = [25, 15]
+
+    old_data = netCDF4.Dataset(filename).variables
+    data = {}
+    greater_than_zero = {}
+    for x, val in reversed(list(enumerate(old_data['time']))):
+        if val > 0:
+            high_ind = x
+            break
+    for x, val in enumerate(old_data['time']):
+        if val > 0:
+            low_ind = x
+            break
+    time_increment = (old_data['time'][high_ind] - old_data['time'][low_ind]) / len(old_data['time'][low_ind:high_ind])
+    print('time_increment: ', time_increment)
+    print(time.localtime(old_data['time'][0]/1000000000))
+    variable_order = ['qc_percent', 'hatch_check', 'missing_data_flag_check', 'safing_check', 'encoder_check', 'detector_check',
+                        'hbb_thermistor_check', 'abb_thermistor_check', 'spike_check', 'hbb_temp_outlier_check', 'abb_temp_outlier_check',
+                        'bst_temp_outlier_check',
+                        'sce_temp_deviation_check', 'hbb_stable_check', 'hbb_covariance_check', 'imaginary_radiance_check',
+                        'detector_temp_check',
+                        'sky_brightness_temp_spectral_averages_ch1_check', 'sky_brightness_temp_spectral_averages_ch2_check',
+                        'hbb_std_dev_check',
+                        'hbb_lw_nen_check', 'hbb_sw_nen_check', 'lw_responsivity_check', 'sw_responsivity_check', ]
+    print('*********')
+    for var_name in old_data:
+        if var_name not in variable_order and var_name != 'time':
+            print(var_name)
+    print('&&&&&&&&&')
+    data['time'] = []
+    qc_percent_num = 0.0
+    for var_name in ['time',] + variable_order:
+        print(var_name)
+        greater_than_zero[var_name] = False
+        if var_name == 'time':
+            for x, curr_time in enumerate(old_data['time']):
+                if curr_time < 0:
+                    if x > 0:
+                        curr_time = data['time'][x-1] + time_increment
+                    else:
+                        curr_time = old_data['time'][low_ind] - time_increment*(low_ind-x)
+                data['time'].append(curr_time)
+        else:
+            data[var_name] = []
+            for x, value in enumerate(old_data[var_name]):
+                if var_name != 'qc_percent' and var_name != 'hatch_check' and var_name != 'missing_data_flag_check':
+                    if data['missing_data_flag_check'][x] == 0.0:
+                        data[var_name].append(value)
+                    else:
+                        data[var_name].append(0.0)
+                else:
+                    data[var_name].append(value)
+                    if var_name == 'qc_percent' and value == 0.0:
+                        qc_percent_num += 1.0
+                if value > 0.0 and not greater_than_zero[var_name]:
+                    greater_than_zero[var_name] = True
+
+    time_array = np.asarray(data['time'])/1000000000.0/60.0/60.0 - data['time'][0]/1000000000.0/60.0/60.0
+    print(type(time_array))
+    for x in range(len(data['missing_data_flag_check'])):
+        if data['missing_data_flag_check'][x] == data['hbb_stable_check'][x] and data['hbb_stable_check'][x] != 0:
+            print(data['missing_data_flag_check'][x], ' ----- ', data['hbb_stable_check'][x])
+    plt.figure(1)
+    curr_plot_num = 0
+    print('creating subplots...')
+    for value in variable_order:
+        if '_check' in value:
+            var_name = value[:-6] + ' '*30
+        else:
+            var_name = value + ' '*30
+        if 'spike' in var_name:
+            var_name = 'igm_spike' + ' '*30
+        if 'sky_brightness_temp_spectral_averages' in var_name:
+            if 'ch1' in var_name:
+                var_name = 'surface_bt_ch1' + ' '*30
+            elif 'ch2' in var_name:
+                var_name = 'surface_bt_ch2' + ' '*30
+        if 'qc_percent' in var_name:
+            var_name = 'qc_percent ({:3.2f}%)'.format(100.0 * qc_percent_num / len(data['qc_percent'])) + ' '*30
+        curr_plot_num += 1
+        if value not in old_data.keys():
+            print('subplot ', curr_plot_num, ' of ', len(variable_order), ' -- ', value, ' --- missing')
+        else:
+            print('subplot ', curr_plot_num, ' of ', len(variable_order), ' -- ', var_name)
+            ax = plt.subplot(len(variable_order), 1, curr_plot_num)
+            # if not greater_than_zero[value]:
+            #     ax.set_axis_bgcolor((0.87,0.87,0.87))
+            ax.set_axis_bgcolor((0.8,1.0,1.0))
+            plt.ylabel(var_name, rotation=0)
+            plt.axis([0, 24, 0, 1])
+            plt.xticks( arange(0,25) )
+            plt.yticks((0,1))
+            plt.setp(ax.xaxis.get_ticklabels()[1::3], visible=False)
+            plt.setp(ax.xaxis.get_ticklabels()[2::3], visible=False)
+            # plt.bar(data['time'], old_data[value], width=time_increment
+            # plt.plot(data['time'], data[value], color='black')
+            # print(value, ' ---- ', data[value])
+
+            if greater_than_zero[value]:
+                plt.fill_between(time_array, data[value], y2=1, color='white')
+                plt.fill_between(time_array, data[value], color='black')
+                # x_vals = arange(data['time'][len(data['time'])-1], 25)
+                # plt.fill_between(x_vals, [1]*len(x_vals), color='yellow')
+            else:
+                ax.set_axis_bgcolor((0.87,0.87,0.87))
+            plt.tight_layout(h_pad=0.1)
+            plt.subplots_adjust(wspace=1.0)
+    print('saving...')
+    # plt.savefig('/Users/adiebold/rooftop_pngs/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/qc/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/archive_test/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/archive_pngs/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/ena_data_pics/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/sgp-c1_pics/' + filename.split('/')[-1].split('.')[0] + '.png')
+    plt.savefig('/Users/adiebold/test_sgp_pngs/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/awr_pngs/' + filename.split('/')[-1].split('.')[0] + '.png')
+
+    print('finished')
+    print('subplot execution time: %d minute(s), %.2f second(s)' % ((time.time() - subplot_start_time) // 60,
+            (time.time() - subplot_start_time) % 60))
+    # plt.show()
+    plt.clf()
+
+if __name__ == '__main__':
+    start_time = time.time()
+    import argparse
+    parser = argparse.ArgumentParser()
+    parser.add_argument('filepath')
+
+    args = parser.parse_args()
+
+    print(args.filepath)
+#skip 23 for 1607*
+#skip 54 for 1608*
+#skip 62 for 1609*
+#skip 91 for 1610*
+    skip_num = 222         #172
+    curr_num = 0
+    if os.path.isdir(args.filepath):
+        for filename in os.listdir(args.filepath):
+            filename = args.filepath + '/' + filename
+            if os.path.isdir(filename):
+                for file_name in os.listdir(filename):
+                    file_name = filename + '/' + file_name
+                    if 'QC.nc' in file_name:
+                        # if not os.path.isfile('/Users/adiebold/sgp_new_pngs/' + file_name.split('/')[-1].split('.')[0] + '.png'):
+                        curr_num += 1
+                        print(curr_num, ' -- ', file_name)
+                        if curr_num >= skip_num:
+                            run(file_name)
+                        else:
+                            print('SKIP')
+            elif 'QC.nc' in filename:
+                run(filename)
+    elif os.path.isfile(args.filepath):
+        if 'QC.nc' in args.filepath:
+            run(args.filepath)
+
+    print('total execution time: %d minute(s), %.2f second(s)' % ((time.time() - start_time) // 60, (time.time() - start_time) % 60))

testing/testing_quick_vis.py

+import netCDF4
+import matplotlib.pyplot as plt
+import matplotlib.patches as patches
+import time
+from scipy import arange
+import os
+import numpy as np
+import pandas as pd
+import argparse
+
+def run(filename):
+
+    subplot_start_time = time.time()
+
+    print('program running...')
+    plt.rcParams['figure.figsize'] = [25, 15]
+
+    old_data = netCDF4.Dataset(filename).variables
+    #need to account for times when the first or last values are not a time (NaT)
+    for x, val in reversed(list(enumerate(old_data['time']))):
+        if val > 0:
+            high_ind = x
+            break
+    for x, val in enumerate(old_data['time']):
+        if val > 0:
+            low_ind = x
+            break
+
+    #calculate average time increment between times for inferring NaT values
+    time_increment = (old_data['time'][high_ind] - old_data['time'][low_ind]) / len(old_data['time'][low_ind:high_ind])
+
+    data = pd.DataFrame(index=range(len(old_data['time'])), columns=old_data.keys())
+
+    #put data into DataFrame
+    for key in old_data:
+        data[key] = old_data[key]
+
+    #makes sure the graphs are put in the correct order
+    variable_order = ['qc_percent', 'hatch_check', 'missing_data_flag_check', 'safing_check', 'encoder_check', 'detector_check',
+                        'hbb_thermistor_check', 'abb_thermistor_check', 'spike_check', 'hbb_temp_outlier_check', 'abb_temp_outlier_check',
+                        'bst_temp_outlier_check',
+                        'sce_temp_deviation_check', 'hbb_stable_check', 'hbb_covariance_check', 'imaginary_radiance_check',
+                        'detector_temp_check',
+                        'sky_brightness_temp_spectral_averages_ch1_check', 'sky_brightness_temp_spectral_averages_ch2_check',
+                        'hbb_std_dev_check',
+                        'hbb_lw_nen_check', 'hbb_sw_nen_check', 'lw_responsivity_check', 'sw_responsivity_check', ]
+
+    #checks for any variables in the DataFrame not accounted for in variable_order
+    for var_name in old_data:
+        if var_name not in variable_order and var_name != 'time':
+            print('*********\n', var_name, '\n&&&&&&&&&')
+        else:
+            data[var_name] = old_data[var_name]
+
+    #infers any times for NaT values
+    data.loc[data['time'] < 0] = 0
+    for t in range(len(data['time'])):
+        if data['time'].iloc[t] == 0 and t != 0:
+            data['time'].iloc[t] = data['time'].iloc[t-1] + time_increment
+        elif data['time'].iloc[t] == 0 and t == 0:
+            for x, val in enumerate(data['time']):
+                if val > 0:
+                    data[t] = val - x*time_increment
+                    break
+
+    #converts times from nanoseconds to hours
+    data['time'] = data['time']/1000000000/60/60 - data['time'].iloc[0]/1000000000/60/60
+
+    data.set_index('time', inplace=True)
+
+    #calculate how much of qc_percent is valid
+    qc_percent_num = 100 * (1 - sum(data['qc_percent']) / len(data['qc_percent']))
+
+    '''
+    data['time'] = []
+    qc_percent_num = 0.0
+    for var_name in ['time',] + variable_order:
+        print(var_name)
+        greater_than_zero[var_name] = False
+        if var_name == 'time':
+            for x, curr_time in enumerate(old_data['time']):
+                if curr_time < 0:
+                    if x > 0:
+                        curr_time = old_data['time'][x-1] + time_increment
+                    else:
+                        curr_time = old_data['time'][low_ind] - time_increment*(low_ind-x)
+                data['time'].append(curr_time)
+        else:
+            data[var_name] = []
+            for x, value in enumerate(old_data[var_name]):
+                if var_name != 'qc_percent' and var_name != 'hatch_check' and var_name != 'missing_data_flag_check':
+                    if data['missing_data_flag_check'][x] == 0.0:
+                        data[var_name].append(value)
+                    else:
+                        data[var_name].append(0)
+                else:
+                    data[var_name].append(value)
+                    if var_name == 'qc_percent' and value == 0:
+                        qc_percent_num += 1.0
+                if value > 0.0 and not greater_than_zero[var_name]:
+                    greater_than_zero[var_name] = True
+    '''
+
+    plt.figure(1)
+    curr_plot_num = 0
+
+    print('creating subplots...')
+    for value in variable_order:
+
+        #alter name to make it better formatted for the graph
+        #use 30 spaces because that's what works best
+        if '_check' in value:
+            var_name = value[:-6] + ' '*30
+        else:
+            var_name = value + ' '*30
+        if 'spike' in var_name:
+            var_name = 'igm_spike' + ' '*30
+        if 'sky_brightness_temp_spectral_averages' in var_name:
+            if 'ch1' in var_name:
+                var_name = 'surface_bt_ch1' + ' '*30
+            elif 'ch2' in var_name:
+                var_name = 'surface_bt_ch2' + ' '*30
+        if 'qc_percent' in var_name:
+            var_name = 'qc_percent ({:3.2f}%)'.format(qc_percent_num) + ' '*30
+
+        curr_plot_num += 1
+
+        #doublecheck that key is actually a key
+        if value not in data.keys():
+            print('subplot ', curr_plot_num, ' of ', len(variable_order), ' -- ', value, ' --- missing')
+        else:
+            print('subplot ', curr_plot_num, ' of ', len(variable_order), ' -- ', var_name)
+            ax = plt.subplot(len(variable_order), 1, curr_plot_num)
+
+            #turn background light blue so anywhere not graphed is distinguishable
+            ax.set_axis_bgcolor((0.8,1.0,1.0))
+            plt.ylabel(var_name, rotation=0)
+            # 0-24 for 24 hours, 0-1 either valid or invalid
+            plt.axis([0, 24, 0, 1])
+
+            # x ticks only on every third to reduce clutter
+            plt.xticks( arange(0,25) )
+            plt.yticks((0,1))
+            plt.setp(ax.xaxis.get_ticklabels()[1::3], visible=False)
+            plt.setp(ax.xaxis.get_ticklabels()[2::3], visible=False)
+
+            '''
+            plt.bar(data['time'], old_data[value], width=time_increment
+            plt.plot(data['time'], data[value], color='black')
+            print(value, ' ---- ', data[value])
+            '''
+
+            #fill in the area below values as black and above as white to get rid
+            #of the blue background
+            if any(data[value] > 0):
+                plt.fill_between(data.index, data[value], y2=1, color='white')
+                plt.fill_between(data.index, data[value], color='black')
+                plt.tight_layout(h_pad=0.1)
+                plt.subplots_adjust(wspace=1.0)
+
+            #if none of the values are above 0 change background to gray
+            else:
+                ax.set_axis_bgcolor((0.87,0.87,0.87))
+
+    print('saving...')
+    # plt.savefig('/Users/adiebold/rooftop_pngs/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/qc/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/archive_test/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/archive_pngs/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/ena_data_pics/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/sgp-c1_pics/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/sgp_pngs/' + filename.split('/')[-1].split('.')[0] + '.png')
+    # plt.savefig('/Users/adiebold/' + filename.split('/')[-1].split('.')[0] + '.png')
+    plt.savefig('/Users/adiebold/awr_pngs/' + filename.split('/')[-1].split('.')[0] + '.png')
+
+    print('finished')
+    print('subplot execution time: %d minute(s), %.2f second(s)' % ((time.time() - subplot_start_time) // 60,
+            (time.time() - subplot_start_time) % 60))
+    # plt.show()
+    plt.clf()
+
+if __name__ == '__main__':
+    start_time = time.time()
+    parser = argparse.ArgumentParser()
+    parser.add_argument('filepath')
+
+    args = parser.parse_args()
+
+    print(args.filepath)
+
+    #amount of files to skip
+    skip_num = 393
+    curr_num = 0
+    if os.path.isdir(args.filepath):
+        for filename in os.listdir(args.filepath):
+            filename = args.filepath + '/' + filename
+            if os.path.isdir(filename):
+                for file_name in os.listdir(filename):
+                    file_name = filename + '/' + file_name
+                    if 'QC.nc' in file_name:
+                        curr_num += 1
+                        print(curr_num, ' -- ', file_name)
+                        if curr_num >= skip_num:
+                            run(file_name)
+                        else:
+                            print('SKIP')
+            elif 'QC.nc' in filename:
+                run(filename)
+    elif os.path.isfile(args.filepath):
+        if 'QC.nc' in args.filepath:
+            run(args.filepath)
+
+    print('total execution time: %d minute(s), %.2f second(s)' % ((time.time() - start_time) // 60, (time.time() - start_time) % 60))