#!/home/mesohorse/anaconda2/bin/python
'''
This code is meant to be used to process ceilo data from ceilometera
other entities

Designed by Luke Leclair-Marzolf
updated extensively  by jdh to avoid having 3 versions of very similar code
and to allow for one code to be used for real time and prior periods
and to switch to gradient as opposed to local max's
'''
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
#import all the modules
import numpy as np
import time, datetime
import calendar
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
from matplotlib import rc
import modified_methods_120917 as m
from itertools import chain
import sys, signal
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

#get which ceilometer from the command line
unit_id = sys.argv[1]
## Define output directories for realtime
dayfiledir = '/home/mesohorse/oper/mesowest/mwingest/datafiles/uunet_ceil/'
appendtxtdir = '/uufs/chpc.utah.edu/common/home/horel-group9/uunet/ceilometer/PBL_txt/'
#realtime flag only applies to copying a plot to "current"
realtime = True
#get end time
end=int(time.time())

if len(sys.argv) >= 3:
	# for past time should be entered as the " the end of the current day 2359" YYMMDDHHMM utc
	# will then do for 0000-2355 local time on that calendar day
	et = sys.argv[2]
	end = time.mktime(time.strptime(et,"%y%m%d%H%M"))
	## Define output directories for pasttime
	dayfiledir = '/home/mesohorse/oper/mesowest/mwingest/datafiles/uunet_ceil/past/'
	appendtxtdir = '/uufs/chpc.utah.edu/common/home/horel-group9/uunet/ceilometer/PBL_txt_past/'
	realtime = False

#time in epoch seconds. 
e=datetime.datetime.fromtimestamp(int(end))

#round to previous 5 min interval  
tdelt=e-datetime.timedelta(minutes=e.minute%5,seconds = e.second)
new_end = str(int(time.mktime(tdelt.timetuple())))
#start is 1 day - 5 min (86400 - 300)
start = str(int(new_end) - 86100)
print 'start', start , 'end', new_end

class mysigtimeout:
	def __init__(self, seconds=1, error_message='Timeout'):
		self.seconds = seconds
		self.error_message = error_message
	def handle_timeout(self, signum, frame):
		raise TimeoutError(self.error_message)
	def __enter__(self):
		signal.signal(signal.SIGALRM, self.handle_timeout)
		signal.alarm(self.seconds)
	def __exit__(self, type, value, traceback):
		signal.alarm(0)

# tunable parameters
depth = 2000 #Adjustable variable for how high do you want to look and plot
d2 = 1250 # for lower plots
#binsize here means 600 sec or 10 min filtering inverval when inTime=True
binsize=600
#vert_val is the vertical interval in m
# fht is the first height used 
#7 means 70 m for vaisala and 14 means  70 m for cs. had to go higher due to noise near ground
#vertbin_val is the vertical smoothing. =2 implies +-1 or 30 m. =4 for logan means +-2 or also 30 m since 5 m vertical interval
#plow is the lowest backscatter 
#phigh is the highest backscatter 
#pint is the backscatter bin interval
#grad_th is the threshold gradient (%/(2*pint))
#no_val_att: check for attenuation/saturation defined as having several values greater than phigh in the column of interest
#offset is the minimum separation between two layers in m
#very arbitrary 
# will average elevation between two layers in proximity
#or if offset = 0 then letting the elevation of each layer be independently defined

print 'unit_id',unit_id

if unit_id == 'LOGCS':
	vert_val = 5 
	fht = 24
	vertbin_val = 4
	stnm = 'LGCUT'
	stnm_lc = 'lgcut'
	plow = -8.1
	phigh = -6.2
	pint = 0.05
	grad_th = 0.75
	no_val_att = 10
	offset = 0
elif unit_id == 'UUCL0':
	vert_val = 10 
	fht = 10
	vertbin_val=2
	stnm = 'URHSC'
	stnm_lc = 'urhsc'
	plow = -7.7
	phigh = -5.3
	pint = 0.05
	grad_th = 0.75
	no_val_att = 5;
	offset =  0
elif unit_id == 'UUCL1':
	vert_val = 10 
	fht = 10
	vertbin_val=2
	stnm = 'MTMET'
	stnm_lc = 'mtmet'
	plow = -7.7
	phigh = -5.3
	pint = 0.05
	grad_th = 0.75
	no_val_att = 5
	offset = 0
elif unit_id == 'UUCLA':
	vert_val = 10 
	fht = 10
	vertbin_val=2
	stnm = 'USDR1'
	stnm_lc = 'usdr1'
	plow = -7.7
	phigh = -5.3
	pint = 0.05
	grad_th = 0.75
	no_val_att = 5;
	offset =  0
elif unit_id == 'UUCLB':
	vert_val = 10 
	fht = 10
	vertbin_val=2
	stnm = 'UUSYR'
	stnm_lc = 'uusyr'
	plow = -7.7
	phigh = -5.3
	pint = 0.05
	grad_th = 0.75
	no_val_att = 5
	offset = 0
else:
	sys.exit(0)

#first elevation in m
fht_val = fht * vert_val
#last index value
lht = depth/vert_val
lht2 = d2/vert_val

#get the data
if unit_id != 'LOGCS':
	timetrack1 = int(time.mktime(time.localtime()))
	with mysigtimeout(seconds=80):
		try:
			data, URL = m.git_api(start, new_end, unit_id, 'BS,STATUS', ceilo=True)
			status=data['STATUS']
			timetrackdiff = int(time.mktime(time.localtime()))-timetrack1
			print('API_RETURN_SUCCESS',unit_id,timetrackdiff)
			if(realtime == False):
				clb1 = []
				clb2 = []
				clb3 = []
				for t in data['STATUS']:
					clb1.append(int(float(t[2])*0.3048))
					clb2.append(int(float(t[3])*0.3048))
					clb3.append(int(float(t[4])*0.3048))
			else:
				clb1=int(float(data['STATUS'][-1][2])*0.3048 )
				clb2=int(float(data['STATUS'][-1][3])*0.3048 )
				clb3=int(float(data['STATUS'][-1][4])*0.3048 )
		except:
			timetrackdiff = int(time.mktime(time.localtime()))-timetrack1
			print('API Call or Data Parsing Error',unit_id,timetrackdiff)
			new_end = -9999
			dt='-9999,-9999'
			sys.exit(0)
else:
	try:
		end1 = str(int(start) + 86400/4)
		data, URL = m.git_api(start, end1, unit_id, 'BS,FIRST_CLOUD,SECOND_CLOUD,THIRD_CLOUD', ceilo=True)
		end2 = str(int(end1) + 86400/4)
		data2, URL = m.git_api(end1, end2, unit_id, 'BS,FIRST_CLOUD,SECOND_CLOUD,THIRD_CLOUD', ceilo=True)
		end3 = str(int(end2) + 86400/4)
		data3, URL = m.git_api(end2, end3, unit_id, 'BS,FIRST_CLOUD,SECOND_CLOUD,THIRD_CLOUD', ceilo=True)
		end4 = str(int(end3) + 86400/4)
		data4, URL = m.git_api(end3, end4, unit_id, 'BS,FIRST_CLOUD,SECOND_CLOUD,THIRD_CLOUD', ceilo=True)
		for k,v in data2.iteritems():
			data[k].extend(v)
		for k,v in data3.iteritems():
			data[k].extend(v)
		for k,v in data4.iteritems():
			data[k].extend(v)
		status={}
		status['cl1']=data['FIRST_CLOUD']
		status['cl2']=data['SECOND_CLOUD']
		status['cl3']=data['THIRD_CLOUD']
		if(realtime == False):
			clb1=[int(float(x)*0.3048) for x in status['cl1']]
			clb2=[int(float(x)*0.3048) for x in status['cl2']]
			clb3=[int(float(x)*0.3048) for x in status['cl3']]
		else:
			clb1=int(float(status['cl1'][-1])*0.3048)
			clb2=int(float(status['cl2'][-1])*0.3048)
			clb3=int(float(status['cl3'][-1])*0.3048)
	except:
		print 'API Call or Data Parsing Error'
		new_end = -9999
		dt='-9999,-9999'
		sys.exit(0)

d=np.array(data['DATTIM']) # array of epoch dates
if d.any():
	d.sort() # you gotta do this cause the API sometimes sends funky time arrays
	tdiff = d[-1]-d[0]
	if tdiff< 3600: #must have an hour of data
		print "less than an hour of data"
		sys.exit(0)
else:
	print "no dates returned"
	sys.exit(0)
	
h=np.array(data['HEIGHT']) # array of Heights
bs=np.array(data['BS']) #  array of backscatter values
#print 'len h, len bs:',len(h),len(bs)
dt = datetime.datetime.strftime(datetime.datetime.utcfromtimestamp(int(new_end)),'%m-%d-%Y,%H:%M:%S')
dtplot = datetime.datetime.strftime(datetime.datetime.fromtimestamp(int(new_end)),'%H%M Local %d %b %Y')

#smooth the data in time and vertical
s_bs,s_d,s_h=m._ComputeFieldMeans(data,binsize,unit_id,inTime=True,continuous=False,vertbin=vertbin_val)
s_bs=np.array(s_bs)
s_d=np.array(s_d)
s_h=np.array(s_h)
s_d.sort()

levels=[]
percents = []
multi = []
gradient = []

#count the gradients in  the entire image
#only going to use the highest big gradient
percent_first,first=m.count_the_highest_gradient(s_bs,unit_id,dt,fht,lht,phigh,plow,pint,grad_th)
no_times = len(s_bs)

#really want to have at least one layer
if first:
	multi.append(first)
	#it's bottom up to find the first layer less than multi[0] at each time 
	#we are limiting to within the vertical levels between fht and lht
	temp1 = (m.tLow(s_bs[:,fht:lht],s_h[fht:lht],multi[0],phigh,pint,no_val_att,vert_val,fht_val))
	#print temp1
	#get the index values for each of the  elevations of the highest layer
	ind_first = temp1/vert_val
	#print 'ind_first',ind_first
	#now get the percentages below the first layer and any other layers
	percents,gradient,multi[1:]=m.count_the_percent_gradient(s_bs,unit_id,dt,fht,ind_first,phigh,plow,pint,grad_th)
	for i in range(len(multi)):
		if i != 0:
			#any layer below the top one is found from the
			#top down from prior layer and will be the first layer that's greater than
			#that threshold until fht
			temp1,levels[i-1] = (m.tHigh(s_bs[:,fht:lht],s_h[fht:lht],multi[i],levels[i - 1],phigh,pint,no_val_att,vert_val,fht_val,offset))

		# for all layers, see if layer is above cloud layers
		if unit_id == 'LOGCS':
			temp2 = [-9999 if temp1[j] > (status['cl1'][j] * 0.3048) and status['cl1'][j] > 0 and temp1[j] > 0 else temp1[j] for j in range(len(temp1))]
			temp2 = [-9999 if temp1[j] > (status['cl2'][j] * 0.3048) and status['cl2'][j] > 0 and temp1[j] > 0 else temp1[j] for j in range(len(temp1))]
			temp2 = [-9999 if temp1[j] > (status['cl3'][j] * 0.3048) and status['cl3'][j] > 0 and temp1[j] > 0 else temp1[j] for j in range(len(temp1))]
		else:
			temp2 = [-9999 if temp1[j] > int(status[j][2]*0.3048 ) and status[j][2] > 0 and temp1[j] > 0 else temp1[j] for j in range(len(temp1))]
			temp2 = [-9999 if temp1[j] > int(status[j][3]*0.3048 ) and status[j][3] > 0 and temp1[j] > 0 else temp1[j] for j in range(len(temp1))]
			temp2 = [-9999 if temp1[j] > int(status[j][4]*0.3048 ) and status[j][4] > 0 and temp1[j] > 0 else temp1[j] for j in range(len(temp1))]
		temp3 = m.eradicate(temp2)
		levels.append(temp3)
	levels = np.array(levels)
else:
	# no levels at any time so set all to 0 elevation
	while len(levels) < len(bs[:,0]):
		levels.append(None)

## Real Time Output
		
if(realtime == True):		
	multstring = ''
	levstring = ''
	try:
		outmultis = multi[::-1]
		outlevels = levels[:,-1]
		outlevels = outlevels[::-1]
		for i in range(0,len(outmultis)):
			multstring = multstring+("%.3f" % outmultis[i])+';'
			levstring = levstring+("%.1f" % outlevels[i])+';'
		multstring = multstring[0:-1]
		levstring = levstring[0:-1]  
		print 'CLOUD LAYERS',clb1,clb2,clb3
	except:
		print 'No print'
	try:
		if(clb1 <= 0): clb1 = -9999
	except:
		clb1 = -9999
	try:
		if(clb2 <= 0): clb2 = -9999
	except:
		clb2 = -9999
	try:
		if(clb3 <= 0): clb3 = -9999
	except:
		clb3 = -9999
	with open(dayfiledir+'PBL_' + stnm + '.txt','a') as fout:
		fout.write(str(new_end)+'|'+str(dt)+'|'+str(multstring)+'|'+str(levstring)+'|'+str(clb1)+'|'+str(clb2)+'|'+str(clb3)+'\n')
	fout.close()
	with open(appendtxtdir+'PBL_' + stnm_lc + '.txt','a') as fout:
		fout.write(str(new_end)+'|'+str(dt)+'|'+str(multstring)+'|'+str(levstring)+'|'+str(clb1)+'|'+str(clb2)+'|'+str(clb3)+'\n')
	fout.close()
	m.all_the_plots(stnm,d, h[:lht],bs[:,:lht], levels, plow,phigh, pint, dtplot,unit_id,lht2,gradient,multi=multi, percent_first=percent_first,percents=percents,d_2=s_d,bs_smoothed=s_bs[:,:lht],realtime=realtime)

## Historical Output
	
else:
	print appendtxtdir+'PBL_' + stnm_lc + '_refill.txt'
	for tt in range(0,len(s_d)):
		multstring = ''
		levstring = ''		
		try:
			outmultis = multi[::-1]
			outlevels = levels[:,tt]
			outlevels = outlevels[::-1]
			for i in range(0,len(outmultis)):
				multstring = multstring+("%.3f" % outmultis[i])+';'
				levstring = levstring+("%.1f" % outlevels[i])+';'
			multstring = multstring[0:-1]
			levstring = levstring[0:-1]  
		except:
			print 'No print'
		outepoch = int(s_d[tt])-(int(s_d[tt]) % 600)
		outdate = datetime.datetime.strftime(datetime.datetime.utcfromtimestamp(outepoch),'%m-%d-%Y,%H:%M:%S')
		clb1out = -9999
		clb2out = -9999
		clb3out = -9999
		try:
			thisindex = np.where(d <= outepoch)[0][-1]
			if(clb1[thisindex] > 0): clb1out = clb1[thisindex]
			if(clb2[thisindex] > 0): clb2out = clb2[thisindex]
			if(clb3[thisindex] > 0): clb3out = clb3[thisindex]
		except:
			skip = 1
		with open(appendtxtdir+'PBL_' + stnm_lc + '_refill.txt','a') as fout:
			fout.write(str(outepoch)+'|'+str(outdate)+'|'+str(multstring)+'|'+str(levstring)+'|'+str(clb1out)+'|'+str(clb2out)+'|'+str(clb3out)+'\n')
		fout.close()
	m.all_the_plots(stnm,d, h[:lht],bs[:,:lht], levels, plow,phigh, pint, dtplot,unit_id,lht2,gradient,multi=multi, percent_first=percent_first,percents=percents,d_2=s_d,bs_smoothed=s_bs[:,:lht],realtime=realtime)
	
sys.exit(0)