#!/usr/bin/python3 import sys sys.path.append("/usr/dd/common/python_lib") import pytz, time_utils from datetime import datetime, timedelta ############################################################################### # Compute daily high 12-hour averages, 12-hour running averages, and maximum # %-TDG for a given time range given atmospheric and TDG pressure time series def compute_avgs( atm, tdg, stime, etime, tz, thresh=None ): ( high12s, roll12s, maxes ) = ( {}, {}, {} ) ( pcts, rolls, rollcts ) = _hourly_tdgs( atm, tdg, stime, etime, tz ) one_day = timedelta( days=1 ) one_hour = timedelta( hours=1 ) # Ugliness here: Step from the start of the day so that we never ask for # an 0100 value for the day DST falls back, because that time occurs twice # in the local time zone, and astimezone() will return the second one. # Add the subtracted hour back to each datetime in the range. for dt in time_utils.local_range( stime - one_hour, etime - one_hour, one_day, tz ): dt += one_hour date = dt.date() ( exceed, flag24, flag36, hi_roll, max ) = ( 0, False, False, None, None ) pct_list = [] end_of_day = time_utils.local_increment( dt, one_day - one_hour, tz ) for time in time_utils.datetime_range( dt, end_of_day, one_hour ): if time in rolls: ( val, ct, flag ) = rolls[time] if flag: flag36 = True avg = val / ct if hi_roll is None or hi_roll < avg: hi_roll = avg if time in pcts: ( pct, valid ) = pcts[time] if pct > 136: flag24 = True if valid: pct_list.append( pct ) if max is None or pct > max: max = pct if thresh is not None and pct > thresh: exceed += 1 if date in rollcts: roll12s[date] = [ hi_roll, rollcts[date], flag36 ] ( count, sum, total ) = ( 0, 0, 0 ) pct_list.sort( reverse=True ) for x in pct_list: total += 1 if total > 12: break sum += x count += 1 if count > 0: high12s[date] = [ sum / count, len( pct_list ), flag24 ] if max is not None: maxes[date] = [ max, exceed, len( pct_list ), flag24 ] return( high12s, roll12s, maxes ) ############################################################################### # Determine if an Atmospheric Pressure or Total Dissolved Gas value should be # considered valid. def is_valid_atm( atm ): return( atm > 609 and atm < 899 ) def is_valid_tdg( tdg ): return( tdg > 601 and tdg < 1250 ) ############################################################################### # Calculate the % Total Dissolved Gas, determine if the value should be # considered valid or not. def percent( tdg, atm ): if atm is None or tdg is None: return( None, None ) pct = None if atm > 0 and tdg > 0: pct = tdg / atm * 100 valid = is_valid_atm( atm ) and is_valid_tdg( tdg ) return( pct, valid ) ############################################################################### # Compute hourly TDG percentages (pcts), 12-hour rolling averages (rolls), and # the number of rolling averages available in a given calendar day (rollcts) def _hourly_tdgs( atm, tdg, stime, etime, tz ): ( pcts, rolls, rollcts ) = ( {}, {}, {} ) one_hour = timedelta( hours=1 ) eleven_hours = timedelta( hours=11 ) for dt in time_utils.datetime_range( stime - eleven_hours, etime, one_hour ): if dt in tdg and dt in atm: ( pct, valid ) = percent( tdg[dt], atm[dt] ) else: ( pct, valid ) = ( None, None ) pcts[dt] = [ pct, valid ] if valid: # Since hours run 1-24, count hour 0 as affecting yesterday. # Count it for today as well, since it will roll forward 12 hours and we # won't make a date change later. if dt.astimezone( tz ).hour == 0: yesterday = ( dt - one_hour ).astimezone( tz ).date() if yesterday not in rollcts: rollcts[yesterday] = 0 rollcts[yesterday] += 1 date = dt.astimezone( tz ).date() if date in rollcts: rollcts[date] += 1 else: rollcts[date] = 1 # Add this hour's percentage to the 12 rolling averages it affects, # keep count for averaging puposes, and flag if over 136% last_affected = dt + eleven_hours for affect in time_utils.datetime_range( dt, last_affected, one_hour ): # Add to a new day's rolling count if we're *going* to affect its # high rolling average next time around. But if 00 tomorrow is our # last affected hour, that's technically just hour 24 of today. # This increment occurs here because we want to do it a maximum of # once per rolling average, right *after* a date transition. if( affect.astimezone( tz ).hour == 0 and affect != dt and affect != last_affected ): date = affect.astimezone( tz ).date() if date in rollcts: rollcts[date] += 1 else: rollcts[date] = 1 if affect not in rolls: rolls[affect] = [ 0, 0, False ] rolls[affect][0] += pct rolls[affect][1] += 1 if pct > 136: rolls[affect][2] = True return( pcts, rolls, rollcts ) ############################################################################### ############################################################################### if __name__ == '__main__': pacific = pytz.timezone( 'US/Pacific' ) atm_val = [ 738.5, 737.6, 736.9, 736.4, 736.2, 736.2, 736.6, 737.2, 737.7, 738.0, 738.4, 739.0, 739.1, 739.4, 739.5, 739.7, 740.2, 740.9, 740.9, 741.2, 742.5, 742.1, 740.9, 740.1, 739.8, 739.2, 738.7, 738.3, 738.2, 738.2, 738.7, 739.3, 739.8, 740.3, 740.6 ] tdg_val = [ 853, 850, 837, 848, 848, 845, 805, 793, 786, 781, 780, 779, 779, 779, 779, 778, 778, 794, 804, 822, 831, 833, 838, 841, 844, 845, 845, 842, 839, 834, 811, 786, 778, 779, 780 ] stime = datetime.strptime( '2015-08-12 14:00', '%Y-%m-%d %H:%M' ) etime = stime + timedelta( hours=34 ) ( atm, tdg ) = ( {}, {} ) for dt in time_utils.datetime_range( stime, etime, timedelta( hours=1 ) ): atm[dt] = atm_val.pop( 0 ) tdg[dt] = tdg_val.pop( 0 ) stime += timedelta( hours=11 ) ( pcts, rolls, rollcts ) = _hourly_tdgs( atm, tdg, stime, etime ) for tm in sorted( pcts.keys() ): print("%s: [%6.3f, %5s] [%8.3f, %2d, %s]" % ( tm, pcts[tm][0], pcts[tm][1], rolls[tm][0], rolls[tm][1], rolls[tm][2] )) for tm in rollcts: print("%s: %s" %( tm, rollcts[tm] )) ( hi, roll, max ) = compute_avgs( atm, tdg, stime, etime, pacific, 120 ) for tm in sorted( hi.keys() ): print(( "%s: [%6.3f, %d, %s], [%6.3f, %d, %s], [%6.3f, %d, %d, %s]" % ( tm, hi[tm][0], hi[tm][1], hi[tm][2], roll[tm][0], roll[tm][1], roll[tm][2], max[tm][0], max[tm][1], max[tm][2], max[tm][3] ) ))