#!/usr/bin/python

# Sorts wind data into events
# This version was started on Jan 22 2009 and will fix mistakes
# related to categorization of events
# It will record data into new files.

import re
import math
from operator import itemgetter
import pickle
import random

# Set wind event threshold values
threshold = 12		# Minimum wind magnitude (knots) required to start an event
duration = 2		# Minimum number of hours (6 hour intervals) required for an event
weakduration = 1	# Maximum duration of weak winds during an event (6 hour intervals) 
synoptictime = 8	# Synoptic timescale required to separate events

springnweventfile = open('springnwevents_04.txt', 'w')
springseeventfile = open('springseevents_04.txt', 'w')
summernweventfile = open('summernwevents_04.txt', 'w')
summerseeventfile = open('summerseevents_04.txt', 'w')
fallnweventfile = open('fallnwevents_04.txt', 'w')
fallseeventfile = open('fallseevents_04.txt', 'w')
winternweventfile = open('winternwevents_04.txt', 'w')
winterseeventfile = open('winterseevents_04.txt', 'w')
springnwpickle = open('springnw_04.pkl', 'wb')
springsepickle = open('springse_04.pkl', 'wb')
summernwpickle = open('summernw_04.pkl', 'wb')
summersepickle = open('summerse_04.pkl', 'wb')
fallnwpickle = open('fallnw_04.pkl', 'wb')
fallsepickle = open('fallse_04.pkl', 'wb')
winternwpickle = open('winternw_04.pkl', 'wb')
wintersepickle = open('winterse_04.pkl', 'wb')
durationsfile = open('durations.txt', 'w')
maxwindsfile = open('maxwinds.txt', 'w')
data = list()		# Create data list container
nwevents = list()		# Create event dictionary list container
seevents = list()		# Create event dictionary list container
springnw = list()
springse = list()
summernw = list()
summerse = list()
fallnw = list()
fallse = list()
winternw = list()
winterse = list()	# Create revised event containers
springnwmaxwinds = list()
springsemaxwinds = list()
summernwmaxwinds = list()
summersemaxwinds = list()
fallnwmaxwinds = list()
fallsemaxwinds = list()
winternwmaxwinds = list()
wintersemaxwinds = list()	# Create ordered maxwinds containers
springnwcomp = list()
springsecomp = list()
summernwcomp = list()
summersecomp = list()
fallnwcomp = list()
fallsecomp = list()
winternwcomp = list()
wintersecomp = list()
springnwcomp2 = list()
springsecomp2 = list()
summernwcomp2 = list()
summersecomp2 = list()
fallnwcomp2 = list()
fallsecomp2 = list()
winternwcomp2 = list()
wintersecomp2 = list()
nwincrementor = 0
seincrementor = 0
nwmonthcount = list((0,0,0,0,0,0,0,0,0,0,0,0))
semonthcount = list((0,0,0,0,0,0,0,0,0,0,0,0))
mcfinal = list((0,0,0,0,0,0,0,0,0,0,0,0))
linenumber = 0
nweventno = 0
seeventno = 0
springnweventno = 0
springseeventno = 0
summernweventno = 0
summerseeventno = 0
fallnweventno = 0
fallseeventno = 0
winternweventno = 0
winterseeventno = 0
sepincrementor = 0

durations = {'winternw': 'winternwdurations=[', 'winterse': 'wintersedurations=[', 'fallnw': 'fallnwdurations=[', 'fallse': 'fallsedurations=[', 'springnw': 'springnwdurations=[', 'springse': 'springsedurations=[', 'summernw': 'summernwdurations=[', 'summerse': 'summersedurations=['}
maxwinds = {'winternw': 'winternwmaxwinds=[', 'winterse': 'wintersemaxwinds=[', 'fallnw': 'fallnwmaxwinds=[', 'fallse': 'fallsemaxwinds=[', 'springnw': 'springnwmaxwinds=[', 'springse': 'springsemaxwinds=[', 'summernw': 'summernwmaxwinds=[', 'summerse': 'summersemaxwinds=['}

def eventwriter(eventlist):
	eventlist.append(prev)	

def eventwriter2(X, eventfile):
	line1 = "Event number " + str(linenumber2) +":\n"
	line2 = X['startDate'] + " to " + X['endDate'] + " : approx " + str(X['duration']*6) + " hours.\n"
	line3 = str(X['maxWind']) + " knots max wind on " + X['maxWindDate'] + ". " + str(X['minDirection']) + " to " + str(X['maxDirection']) + " degrees.\n"
	line4 = "\n"		
	eventfile.write(line1)
	eventfile.write(line2)
	eventfile.write(line3)
	eventfile.write(line4)

def isodd(n):
	return n%2

def mwrite(maxwindlist, eventfile):
	l = len(maxwindlist)
	if isodd(l):
		a = int(math.ceil(l/2))
		m = maxwindlist[a]
	else:
		a = maxwindlist[l//2] + maxwindlist[l//2 + 1]
		m = a/2

	mline = "Median: " + str(m) + "\n"
	eventfile.write(mline)
	print mline
	return m

# R.E. for combined mag and dirn data. 1=yr 2=mn 3=dd 4=hr 5=mag 6=dirn
patt=re.compile("^(\d\d)\s+(\d\d)\s+(\d\d)\s+(\d\d)\s+(\d+\.\d+)\s+(\d+\.\d+)$")


# Open all data
magfile = open('mag.txt', 'r')
dirnfile = open('dirn.txt', 'r')
while 1:	# Add a tab and direction data to the end of each line of magnitude data
	magline = magfile.readline()
	if not magline: break
	dirnline = dirnfile.readline()
	if not dirnline: break
	line = magline + '\t' + dirnline
	data.append(line)


# Test for wind events
for i in range(len(data)):
	line = data[i]
	if not line: break
	linenumber+=1
	matchobj = patt.match(line)
	if matchobj:
		# If the pattern matches then set matched values
		yr = matchobj.group(1)	
		mn = matchobj.group(2)
		dd = matchobj.group(3)
		hr = matchobj.group(4)
		mag = eval(matchobj.group(5))	# Mag and Direction set as numbers
		dirn = eval(matchobj.group(6))
		if mag >= threshold and 270 <= dirn <= 360:	
			# If the wind magnitude is event-strength (as previously defined) and northwesterly
			nwincrementor+=1
			seincrementor = 0
			nwlastdate = yr + " " + mn + " " + dd + " " + hr
			nwlastline = linenumber
			mnn = int(mn)
			nwmonthcount[mnn-1]+=1
			if nwincrementor == 1:	
				# If this is the first consecutive event-strength wind
				nwfirstdate = yr + " " + mn + " " + dd + " " + hr
				nwfirstline = linenumber
				nwdirnmin = dirn
				nwdirnmax = dirn

			if mag >= nwmaxwind:	
				# Set the maximum wind if it exists 
				nwmaxwind = mag
				nwmaxwinddate = yr + " " + mn + " " + dd + " " + hr

			if dirn < nwdirnmin:	
				nwdirnmin = dirn
			elif dirn > nwdirnmax:
				nwdirnmax = dirn

		elif mag >= threshold and 90 <= dirn <= 180:
			# If the wind magnitude is event-strength (as previously defined) and southeasterly
			seincrementor+=1
			nwincrementor = 0
			selastdate = yr + " " + mn + " " + dd + " " + hr
			selastline = linenumber
			mnn = int(mn)
			semonthcount[mnn-1]+=1
			if seincrementor == 1:	
				# If this is the first consecutive event-strength wind
				sefirstdate = yr + " " + mn + " " + dd + " " + hr
				sefirstline = linenumber
				sedirnmin = dirn
				sedirnmax = dirn

			if mag >= semaxwind:	
				# Set the maximum wind if it exists 
				semaxwind = mag
				semaxwinddate = yr + " " + mn + " " + dd + " " + hr

			if dirn < sedirnmin:	
				sedirnmin = dirn
			elif dirn > sedirnmax:
				sedirnmax = dirn

		else:
			# If the wind magnitude is not event-strength
			if nwincrementor >= duration:
				nwevent = {'startDate': nwfirstdate, 'endDate': nwlastdate, 'monthcount': nwmonthcount, 'maxWind': nwmaxwind, 'maxWindDate': nwmaxwinddate, 'minDirection': nwdirnmin, 'maxDirection': nwdirnmax, 'firstLine': nwfirstline, 'lastLine': nwlastline, 'duration': nwincrementor}
				nwevents.append(nwevent)
				#print nwmonthcount #checker
				#print nwmonthcount.index(max(nwmonthcount))+1
			elif seincrementor >= duration:
				seevent = {'startDate': sefirstdate, 'endDate': selastdate, 'monthcount': semonthcount, 'maxWind': semaxwind, 'maxWindDate': semaxwinddate, 'minDirection': sedirnmin, 'maxDirection': sedirnmax, 'firstLine': sefirstline, 'lastLine': selastline, 'duration': seincrementor}
				seevents.append(seevent)
				if semonthcount.index(max(semonthcount))+1 == 12 or semonthcount.index(max(semonthcount))+1 == 1 or semonthcount.index(max(semonthcount))+1 == 2:
					print seevent['startDate']
					print semonthcount #checker
					print semonthcount.index(max(semonthcount))+1

			nwincrementor = 0
			seincrementor = 0
			nwmonthcount = list((0,0,0,0,0,0,0,0,0,0,0,0))
			semonthcount = list((0,0,0,0,0,0,0,0,0,0,0,0))
			nwmaxwind = 0
			semaxwind = 0

# Join events that are separated by less than weakduration
# If two events are separated by less than synoptictime, throw out the event with the smallest maxwind
# Events are printed to a new event file, divided up depending on season and wind direction 
for i in range(1,len(nwevents)):
	curr = nwevents[i]
	prev = nwevents[i-1]
	linesep = curr['firstLine']-prev['lastLine']-1
	if linesep > synoptictime:	
		# If the time between 2 events is greater than synoptic time scale, sort and print previous event
		majmonth = prev['monthcount'].index(max(prev['monthcount']))+1
		if majmonth == 3 or majmonth == 4 or majmonth == 5:
			season = 'spring'
		elif majmonth == 6 or majmonth == 7 or majmonth == 8:
			season = 'summer'
		elif majmonth == 9 or majmonth == 10 or majmonth == 11:
			season = 'fall'
		elif majmonth == 12 or majmonth == 1 or majmonth == 2:
			season = 'winter'

		if season == 'spring':
			springnw.append(prev)
			springnweventno+=1

		elif season == 'summer':
			summernw.append(prev)
			summernweventno+=1

		elif season == 'fall':
			fallnw.append(prev)
			fallnweventno+=1
		
		elif season == 'winter':
			winternw.append(prev)
			winternweventno+=1

		sepincrementor = 0
	elif linesep <= weakduration:
		# If the time between 2 events is less than or equal to weakduration, join the previous and current event
		curr['startDate'] = prev['startDate']
		curr['firstLine'] = prev['firstLine']
		for mci in range(0,12):
			mcfinal[mci] = prev['monthcount'][mci] + curr['monthcount'][mci]

		#print 'prev: ' + str(prev['monthcount'])
		#print 'curr: ' + str(curr['monthcount'])
		#print 'final: ' + str(mcfinal) #checker
		curr['monthcount'] = mcfinal
		curr['duration'] = curr['duration']+prev['duration']+weakduration
		sepincrementor+=1
		if curr['maxWind'] < prev['maxWind']:
			# If the current event's max wind is less than the previous, set max wind to previous value
			curr['maxWind'] = prev['maxWind']
			curr['maxWindDate'] = prev['maxWindDate']

		if curr['minDirection'] > prev['minDirection']:
			curr['minDirection'] = prev['minDirection']
			
		if curr['maxDirection'] < prev['maxDirection']:
			curr['maxDirection'] = prev['maxDirection']

		nwevents[i] = curr

	else:
		# If the two events are during the same synoptic timescale, throw out the weakest event
		if curr['maxWind'] >= prev['maxWind']:
			nwevents[i] = curr
		elif curr['maxWind'] < prev['maxWind']:
			nwevents[i] = prev

	if i == len(nwevents):
		majmonth = curr['monthcount'].index(max(curr['monthcount']))+1
		if majmonth == 3 or majmonth == 4 or majmonth == 5:
			season = 'spring'
		elif majmonth == 6 or majmonth == 7 or majmonth == 8:
			season = 'summer'
		elif majmonth == 9 or majmonth == 10 or majmonth == 11:
			season = 'fall'
		elif majmonth == 12 or majmonth == 1 or majmonth == 2:
			season = 'winter'

		if season == 'spring':
			springnw.append(curr)
			springnweventno+=1

		elif season == 'summer':
			summernw.append(curr)
			summernweventno+=1

		elif season == 'fall':
			fallnw.append(curr)
			fallnweventno+=1
		
		elif season == 'winter':
			winternw.append(curr)
			winternweventno+=1
		
for i in range(1,len(seevents)):
	curr = seevents[i]
	prev = seevents[i-1]
	linesep = curr['firstLine']-prev['lastLine']-1
	if linesep > synoptictime:	
		# If the time between 2 events is greater than synoptic time scale, sort and print previous event
		majmonth = prev['monthcount'].index(max(prev['monthcount']))+1
		if majmonth == 3 or majmonth == 4 or majmonth == 5:
			season = 'spring'
		elif majmonth == 6 or majmonth == 7 or majmonth == 8:
			season = 'summer'
		elif majmonth == 9 or majmonth == 10 or majmonth == 11:
			season = 'fall'
		elif majmonth == 12 or majmonth == 1 or majmonth == 2:
			season = 'winter'

		if season == 'spring':
			springse.append(prev)
			springseeventno+=1

		elif season == 'summer':
			summerse.append(prev)
			summerseeventno+=1

		elif season == 'fall':
			fallse.append(prev)
			fallseeventno+=1
		
		elif season == 'winter':
			winterse.append(prev)
			winterseeventno+=1

		sepincrementor = 0
	elif linesep <= weakduration:
		# If the time between 2 events is less than or equal to weakduration, join the previous and current event
		curr['startDate'] = prev['startDate']
		curr['firstLine'] = prev['firstLine']
		for mci in range(0,12):
			mcfinal[mci] = prev['monthcount'] + curr['monthcount']

		curr['monthcount'] = mcfinal
		curr['duration'] = curr['duration']+prev['duration']+weakduration
		sepincrementor+=1
		if curr['maxWind'] < prev['maxWind']:
			# If the current event's max wind is less than the previous, set max wind to previous value
			curr['maxWind'] = prev['maxWind']
			curr['maxWindDate'] = prev['maxWindDate']

		if curr['minDirection'] > prev['minDirection']:
			curr['minDirection'] = prev['minDirection']
			
		if curr['maxDirection'] < prev['maxDirection']:
			curr['maxDirection'] = prev['maxDirection']

		seevents[i] = curr

	else:
		# If the two events are during the same synoptic timescale, throw out the weakest event
		if curr['maxWind'] >= prev['maxWind']:
			seevents[i] = curr
		elif curr['maxWind'] < prev['maxWind']:
			seevents[i] = prev

	if i == len(seevents):
		majmonth = curr['monthcount'].index(max(curr['monthcount']))
		if majmonth == 3 or majmonth == 4 or majmonth == 5:
			season = 'spring'
		elif majmonth == 6 or majmonth == 7 or majmonth == 8:
			season = 'summer'
		elif majmonth == 9 or majmonth == 10 or majmonth == 11:
			season = 'fall'
		elif majmonth == 12 or majmonth == 1 or majmonth == 2:
			season = 'winter'

		if season == 'spring':
			springse.append(curr)
			springseeventno+=1

		elif season == 'summer':
			summerse.append(curr)
			summerseeventno+=1

		elif season == 'fall':
			fallse.append(curr)
			fallseeventno+=1
		
		elif season == 'winter':
			winterse.append(curr)
			winterseeventno+=1
		
print "length of spring nw array:",len(springnw)," length of spring se array:",len(springse)
print "length of summer nw array:",len(summernw)," length of summer se array:",len(summerse)
print "length of fall nw array:",len(fallnw)," length of fall se array:",len(fallse)
print "length of winter nw array:",len(winternw)," length of winter se array:",len(winterse)

allseasons = list()
allseasons = (springnw, springse, summernw, summerse, fallnw, fallse, winternw, winterse)
si = 0
for s in range(len(allseasons)):
	linenumber2 = 0
	currseason = allseasons[s]
	if not currseason: break
	currseasonordered = sorted(currseason, key=itemgetter('maxWind'), reverse=True)
	for j in range(len(currseasonordered)):
		x = currseasonordered[j]
		y = x['maxWind']
		z = {'startDate': x['startDate'], 'endDate': x['endDate'], 'maxWind': x['maxWind'], 'maxWindDate': x['maxWindDate']}
		if not x: break
		if not y: break
		if not z: break
		linenumber2+=1
		if s == 0:
			eventwriter2(x, springnweventfile)
			springnwmaxwinds.append(y)
			springnwcomp.append(z)
			durations['springnw'] = durations['springnw'] + str(x['duration']) + ","
			maxwinds['springnw'] = maxwinds['springnw'] + str(x['maxWind']) + ","
		elif s == 1:
			eventwriter2(x, springseeventfile)
			springsemaxwinds.append(y)
			springsecomp.append(z)
			durations['springse'] = durations['springse'] + str(x['duration']) + ","
			maxwinds['springse'] = maxwinds['springse'] + str(x['maxWind']) + ","
		elif s == 2:
			eventwriter2(x, summernweventfile)
			summernwmaxwinds.append(y)
			summernwcomp.append(z)
			durations['summernw'] = durations['summernw'] + str(x['duration']) + ","
			maxwinds['summernw'] = maxwinds['summernw'] + str(x['maxWind']) + ","
		elif s == 3:
			eventwriter2(x, summerseeventfile)
			summersemaxwinds.append(y)
			summersecomp.append(z)
			durations['summerse'] = durations['summerse'] + str(x['duration']) + ","
			maxwinds['summerse'] = maxwinds['summerse'] + str(x['maxWind']) + ","
		elif s == 4:
			eventwriter2(x, fallnweventfile)
			fallnwmaxwinds.append(y)
			fallnwcomp.append(z)
			durations['fallnw'] = durations['fallnw'] + str(x['duration']) + ","
			maxwinds['fallnw'] = maxwinds['fallnw'] + str(x['maxWind']) + ","
		elif s == 5: 
			eventwriter2(x, fallseeventfile)
			fallsemaxwinds.append(y)
			fallsecomp.append(z)
			durations['fallse'] = durations['fallse'] + str(x['duration']) + ","
			maxwinds['fallse'] = maxwinds['fallse'] + str(x['maxWind']) + ","
		elif s == 6:
			eventwriter2(x, winternweventfile)
			winternwmaxwinds.append(y)
			winternwcomp.append(z)
			durations['winternw'] = durations['winternw'] + str(x['duration']) + ","
			maxwinds['winternw'] = maxwinds['winternw'] + str(x['maxWind']) + ","
		elif s == 7:
			eventwriter2(x, winterseeventfile)
			wintersemaxwinds.append(y)
			wintersecomp.append(z)
			durations['winterse'] = durations['winterse'] + str(x['duration']) + ","
			maxwinds['winterse'] = maxwinds['winterse'] + str(x['maxWind']) + ","

# Write durations and maxwinds files
durationsfile.write(durations['springnw'] + "\n" + durations['springse'] + "\n" + durations['summernw'] + "\n" + durations['summerse'] + "\n" + durations['fallnw'] + "\n" + durations['fallse'] + "\n" + durations['winternw'] + "\n" + durations['winterse'] + "\n")
durationsfile.close()
maxwindsfile.write(maxwinds['springnw'] + "\n" + maxwinds['springse'] + "\n" + maxwinds['summernw'] + "\n" + maxwinds['summerse'] + "\n" + maxwinds['fallnw'] + "\n" + maxwinds['fallse'] + "\n" + maxwinds['winternw'] + "\n" + maxwinds['winterse'] + "\n")
maxwindsfile.close()

# Write median maxwind to the end of each event file and object
mspringnw = mwrite(springnwmaxwinds, springnweventfile)
mspringse = mwrite(springsemaxwinds, springseeventfile)
msummernw = mwrite(summernwmaxwinds, summernweventfile)
msummerse = mwrite(summersemaxwinds, summerseeventfile)
mfallnw = mwrite(fallnwmaxwinds, fallnweventfile)
mfallse = mwrite(fallsemaxwinds, fallseeventfile)
mwinternw = mwrite(winternwmaxwinds, winternweventfile)
mwinterse = mwrite(wintersemaxwinds, winterseeventfile)

allseasons2 = list()
allseasons2 = (springnwcomp, springsecomp, summernwcomp, summersecomp, fallnwcomp, fallsecomp, winternwcomp, wintersecomp)
medians = (mspringnw, mspringse, msummernw, msummerse, mfallnw, mfallse, mwinternw, mwinterse)
for s in range(len(allseasons2)):
	allevents = allseasons2[s]
	if not allevents: break
	if len(allevents) > 50:
		for k in range(len(allevents)):
			if allevents[k]['maxWind'] > medians[s]:
				# Check next value
				if allevents[k+1]['maxWind'] <= medians[s]:
					mindex = k+1
					break
		
		a_index = sorted(random.sample(range(0,mindex+1), 25))
		b_index = sorted(random.sample(range(mindex+1,len(allevents)), 25))
		for n in a_index:
			if s == 0:
				springnwcomp2.append(allevents[n])
			elif s == 1:
				springsecomp2.append(allevents[n])
			elif s == 2:
				summernwcomp2.append(allevents[n])
			elif s == 3:
				summersecomp2.append(allevents[n])
			elif s == 4:
				fallnwcomp2.append(allevents[n])
			elif s == 5:
				fallsecomp2.append(allevents[n])
			elif s == 6:
				winternwcomp2.append(allevents[n])
			elif s == 7:
				wintersecomp2.append(allevents[n])

		for n in b_index:
			if s == 0:
				springnwcomp2.append(allevents[n])
			elif s == 1:
				springsecomp2.append(allevents[n])
			elif s == 2:
				summernwcomp2.append(allevents[n])
			elif s == 3:
				summersecomp2.append(allevents[n])
			elif s == 4:
				fallnwcomp2.append(allevents[n])
			elif s == 5:
				fallsecomp2.append(allevents[n])
			elif s == 6:
				winternwcomp2.append(allevents[n])
			elif s == 7:
				wintersecomp2.append(allevents[n])
	
	else:
		if s == 0:
			springnwcomp2 = springnwcomp
		elif s == 1:
			springsecomp2 = springsecomp
		elif s == 2:
			summernwcomp2 = summernwcomp
		elif s == 3:
			summersecomp2 = summersecomp
		elif s == 4:
			fallnwcomp2 = fallnwcomp
		elif s == 5:
			fallsecomp2 = fallsecomp
		elif s == 6:
			winternwcomp2 = winternwcomp
		elif s == 7:
			wintersecomp2 = wintersecomp

# Pickling!
pickle.dump(springnwcomp2, springnwpickle, 0)
pickle.dump(springsecomp2, springsepickle, 0)
pickle.dump(summernwcomp2, summernwpickle, 0)
pickle.dump(summersecomp2, summersepickle, 0)
pickle.dump(fallnwcomp2, fallnwpickle, 0)
pickle.dump(fallsecomp2, fallsepickle, 0)
pickle.dump(winternwcomp2, winternwpickle, 0)
pickle.dump(wintersecomp2, wintersepickle, 0)





springnweventfile.close()
springseeventfile.close()
summernweventfile.close()
summerseeventfile.close()
fallnweventfile.close()
fallseeventfile.close()
winternweventfile.close()
winterseeventfile.close()
springnwpickle.close()
springsepickle.close()
summernwpickle.close()
summersepickle.close()
fallnwpickle.close()
fallsepickle.close()
winternwpickle.close()
wintersepickle.close()

print 'End.'