rendered paste body#!/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 reimport mathfrom operator import itemgetterimport pickleimport random# Set wind event threshold valuesthreshold = 12 # Minimum wind magnitude (knots) required to start an eventduration = 2 # Minimum number of hours (6 hour intervals) required for an eventweakduration = 1 # Maximum duration of weak winds during an event (6 hour intervals) synoptictime = 8 # Synoptic timescale required to separate eventsspringnweventfile = 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 containernwevents = list() # Create event dictionary list containerseevents = list() # Create event dictionary list containerspringnw = list()springse = list()summernw = list()summerse = list()fallnw = list()fallse = list()winternw = list()winterse = list() # Create revised event containersspringnwmaxwinds = list()springsemaxwinds = list()summernwmaxwinds = list()summersemaxwinds = list()fallnwmaxwinds = list()fallsemaxwinds = list()winternwmaxwinds = list()wintersemaxwinds = list() # Create ordered maxwinds containersspringnwcomp = 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 = 0seincrementor = 0nwmonthcount = 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 = 0nweventno = 0seeventno = 0springnweventno = 0springseeventno = 0summernweventno = 0summerseeventno = 0fallnweventno = 0fallseeventno = 0winternweventno = 0winterseeventno = 0sepincrementor = 0durations = {'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%2def 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=dirnpatt=re.compile("^(\d\d)\s+(\d\d)\s+(\d\d)\s+(\d\d)\s+(\d+\.\d+)\s+(\d+\.\d+)$")# Open all datamagfile = 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 eventsfor 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 = 0for 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 filesdurationsfile.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 objectmspringnw = 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.'