From be34b9cd87b2987d99e9936758a48cd7d6cd3000 Mon Sep 17 00:00:00 2001 From: kgao <kenny.gao@ssec.wisc.edu> Date: Tue, 23 Aug 2016 15:35:06 +0000 Subject: [PATCH] Code that creates quicklooks creates aoss tower quicklooks --- aosstower/tower_quicklooks/.DS_Store | Bin 0 -> 6148 bytes .../tower_quicklooks/create_quicklook.py | 590 ++++++++++++++++++ aosstower/tower_quicklooks/test_cases.txt | 7 + 3 files changed, 597 insertions(+) create mode 100644 aosstower/tower_quicklooks/.DS_Store create mode 100644 aosstower/tower_quicklooks/create_quicklook.py create mode 100644 aosstower/tower_quicklooks/test_cases.txt diff --git a/aosstower/tower_quicklooks/.DS_Store b/aosstower/tower_quicklooks/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6 GIT binary patch literal 6148 zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**<q8>++&mCkOWA81W14cNZ<zv;LbK1Poaz?KmsK2CSc!( z0ynLxE!0092;Krf2c+FF_Fe*7ECH>lEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0<F0fCPF1$Cyrb|F7^5{eNG?83~ZUUlGt@xh*qZDeu<Z%US-OSsOPv j)R!Z4KLME7ReXlK;d!wEw5GODWMKRea10D2@KpjYNUI8I literal 0 HcmV?d00001 diff --git a/aosstower/tower_quicklooks/create_quicklook.py b/aosstower/tower_quicklooks/create_quicklook.py new file mode 100644 index 0000000..d799d64 --- /dev/null +++ b/aosstower/tower_quicklooks/create_quicklook.py @@ -0,0 +1,590 @@ +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() diff --git a/aosstower/tower_quicklooks/test_cases.txt b/aosstower/tower_quicklooks/test_cases.txt new file mode 100644 index 0000000..f39ab64 --- /dev/null +++ b/aosstower/tower_quicklooks/test_cases.txt @@ -0,0 +1,7 @@ +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 -- GitLab