'''

'''

from numpy import *
from commands import getoutput

dom = [143,209]

def write_sno(file):

	ncl_content = '\
	begin							\n\
	in = addfile("%s","r")			\n\
	swe = in->SNOW					\n\
	asciiwrite("SNW_test.txt", swe)	\n\
	end' % file

	f = open("temp_ncl.ncl","w")
	f.write(ncl_content)
	f.close()

	getoutput('ncl temp_ncl.ncl')


def moving_avg(list, points):
	
	if points%2 == 0:
		padding = (points)/2
	else:
		padding = (points - 1)/2
	
	moving_avgs = []
	for i in range(padding):
		moving_avgs.append(list[i])
	
	for i in range(padding, len(list)-padding):
		
		sample = list[i-padding:i+padding]
		moving_avgs.append(average(sample))
		
	for i in range(padding):
		moving_avgs.append(float('nan'))
			
	return moving_avgs


def snoLine(file, threshold=5.):
	
	f = open(file, "r")
	raw = map(float, f)
	f.close()
	
	for i in range(len(raw)):
		if raw[i] > 5.:
			raw[i] = 5.
			
	fiveDeg_y = 1#18	# From the bottom to the top of the WRF domain is about 39 degrees
					# Divide that by 143 and you get .27 degrees per y-position
					# So that's about 18 grid spaces per 5 degrees latitude
					
	fiveDeg_x = 12	# Similar process which looks at longitudes at either end of domain
					# at 45 degrees latitude
	
	sno = array(raw).reshape(dom[0],dom[1])
	
	snow_line = []
	for x in range(len(sno[0])):
		
		y_slice = []
		for y in range(len(sno)):
			y_slice.append(sno[y][x])
			
		for y in range(len(y_slice)):
			if y > len(y_slice)-fiveDeg_y:
				if y_slice[y] >= threshold:
					snow_line.append(y)
					found = True
					break
			elif y_slice[y] >= threshold and average(y_slice[y:y+fiveDeg_y]) >= threshold:
				snow_line.append(y)
				found = True
				break
			else:
				found = False
		if found == False:
			snow_line.append(float('nan'))
	
	return moving_avg(snow_line, fiveDeg_x)


def test_plot(file, threshold=5.):
	
	import matplotlib.pyplot as plt
	
	f = open(file, "r")
	raw = map(float, f)
	f.close()
	
	for i in range(len(raw)):
		if raw[i] > 1e+10:
			raw[i] = 0.
		elif raw[i] > threshold*2:
			raw[i] = threshold*2	
	
	sno = array(raw).reshape(dom[0],dom[1])
	
	SL = snoLine(file)
	
	plt.figure()
	
	plt.imshow(sno,cmap='gray',interpolation='nearest')
	
	plt.plot(range(209),SL,'r-',linewidth=3.)
	
	plt.ylim(0,143)
	plt.xlim(0,209)
	
	plt.show()
	

def delta_plot(file, d_file, threshold=5.):
	
	import matplotlib.pyplot as plt
	
	f = open(file, "r")
	raw = map(float, f)
	f.close()
	
	f = open(d_file, "r")
	dlt = map(float, f)
	f.close()
	
	for i in range(len(raw)):
		if raw[i] > 1e+10:
			raw[i] = 0.
		elif raw[i] > threshold*2:
			raw[i] = threshold*2	
	
	sno = array(raw).reshape(dom[0],dom[1])
	
	SL = moving_avg(snoLine(file), 20)
	d_SL = array(SL)+array(dlt)
	for i in range(len(d_SL)):
		if d_SL[i] > -1:
			d_SL[i] = int(d_SL[i])
	
	plt.figure()
	
	plt.imshow(sno,cmap='gray',interpolation='nearest')
	
	plt.plot(range(209),SL,'r-',linewidth=3.)
	plt.plot(range(209),d_SL,'b-',linewidth=3.)
	
	plt.ylim(0,143)
	plt.xlim(0,209)
	
	plt.show()


def delta_writer(alter_file, d_file, threshold=5.):
	
	f = open(alter_file, "r")
	raw = map(float, f)
	f.close()
	
	f = open(d_file, "r")
	dlt = map(float, f)
	f.close()
	
	for i in range(len(raw)):
		if raw[i] > 1e+10:
			raw[i] = 0.
		elif raw[i] > threshold*2:
			raw[i] = threshold*2	
	
	sno = array(raw).reshape(dom[0],dom[1])
	
	SL = moving_avg(snoLine(alter_file), 20)
	d_SL = array(SL)+array(dlt)
	for i in range(len(d_SL)):
		if d_SL[i] > -1:
			d_SL[i] = int(d_SL[i])
			
			for j in range(d_SL[i]):
				sno[j][i] = 0.
		
		else:
			# Remove snow beyond boundaries of deltas by border snow extent value
			if i < len(d_SL)/2.:
				for j in d_SL:
					if j > -1:
						y_val = j
						break
			else:
				for j in range(len(d_SL)):
					if d_SL[-j] > -1:
						y_val = d_SL[-j]
						break
			for j in range(int(y_val)):
				sno[j][i] = 0.
				
	import matplotlib.pyplot as plt
	
	plt.figure()
	
	plt.imshow(sno,cmap='gray',interpolation='nearest')
	
	plt.plot(range(209),d_SL,'b-',linewidth=3.)
	
	plt.ylim(0,143)
	plt.xlim(0,209)
	
	plt.show()