#! /usr/bin/env python
#
# ===========================================================
#
# Alexander R Davies
#
# Ocean Exploration, Remote Sensing, and Biogeography Lab
# School of Marine Science and Policy
# College of Earth, Ocean, and Environment
# University of Delaware
# ardavies@udel.edu
#
# Latest Update: 10/09/13
#
# NOTES: 1) 
#
#			- Reads Processed Data, Linearly Interpolates each
#			for each profile to a known grid, writes the
#			linearly interpolated profiles (by variable)
#			to csv files
#
#			- First column in the Ygrid, across are the
#				individual profiles
#
#
# ===========================================================
# Import Modules
# ===========================================================
#
import numpy as np
import csv
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import glob
import os
from matplotlib.ticker import MaxNLocator
from matplotlib import rc, rcParams
from numpy.random import uniform, seed
from matplotlib.mlab import griddata
from numpy import arange,array,ones#,random,linalg
from pylab import plot,show
from scipy import interpolate
from matplotlib.colors import LogNorm
from matplotlib.backends.backend_pdf import PdfPages
#
# ===========================================================
# Function to read csv files
# ===========================================================
#
def csvread(filename):
	with open(filename, 'rb') as f:
	    reader = csv.reader(f, delimiter=',')
	    nrow = 0;
	    for row in reader:
	    	nrow = nrow +1
	    	ncol = len(row)
	array = np.zeros([ncol,nrow]) 
	with open(filename, 'rb') as f:
	    reader = csv.reader(f, delimiter=',')
	    counter = 0
	    for row in reader:
	    	for jj in range(0,ncol):
	    		array[jj,counter] = row[jj]
	    	counter = counter + 1
	return array, nrow, ncol   
#
# ===========================================================
# Function to Find the nearest value
# ===========================================================
#
def find_nearest(array,value):
    idx = (np.abs(array-value)).argmin()
    return array[idx], idx
#
# ===========================================================
# Sort Data and Descriptive Arrays
# ===========================================================
#
os.chdir('C:/Documents and Settings/Alex/Documents/Research/Delaware/MUSTACHE/Jan01_Jun04_2013/data/type')
fullnames = glob.glob('Full*')
gradnames = glob.glob('Grad*')
infonames = glob.glob('Text*')
arraylen =  len(fullnames)
floatdate = np.zeros(arraylen)
lat = [None]*arraylen
lon = [None]*arraylen
day = np.zeros(arraylen)
month=[None]*arraylen
year=[None]*arraylen
daystr=[None]*arraylen
#
# ===========================================================
# Build the Arrays as days past 1/1/13
# ===========================================================
#
for d in range(0,int(arraylen)):
	f = open(infonames[d], 'r')
	while 1:
		line = f.readline()
		if line[:4] == "Date":
			month[d] = str(line[5:7])
			day[d] = int(line[8:10])
			daystr[d] = str(line[8:10])
			year[d] = str(line[11:15])
			break
	f = open(infonames[d], 'r')
	while 1:
		line = f.readline()
		if line[:3] == "Lat":
			lat[d] = str(line[5:12])
			break
	f = open(infonames[d], 'r')
	while 1:
		line = f.readline()
		if line[:3] == "Lon":
			lon[d] = str(line[5:12])
			break
	
	if (year[d] == '2013'):
		yr = 0
	elif (year[d] == '2014'):
		yr = 365
	elif (year[d] == '2015'):
		yr = 365*2
	if(month[d] == '01'):
		mon = 0
	elif(month[d] == '02'):
		mon = 31
	elif(month[d] == '03'):
		mon = 31 + 28
	elif(month[d] == '04'):
		mon = (31 + 28 + 31)
	elif(month[d] == '05'):
		mon = (31 + 28 + 31 + 30)
	elif(month[d] == '06'):
		mon = (31 + 28 + 31 + 30 + 31)
	elif(month[d] == '07'):
		mon = (31 + 28 + 31 + 30 + 31 + 30)
	elif(month[d] == '08'):
		mon = (31 + 28 + 31 + 30 + 31 + 30 + 31)
	elif(month[d] == '09'):
		mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31)
	elif(month[d] == '10'):
		mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30)
	elif(month[d] == '11'):
		mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31)
	elif(month[d] == '12'):
		mon = (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30)
	else:
		print 'invalid month'

	floatdate[d] = yr + mon + day[d]
#
# ===========================================================
# Read the Data and Interpolate for Gridding
# ===========================================================
#
Ygrid = np.linspace(-1800,-10,1791*2)
interplength = len(Ygrid)
GridData = np.zeros([interplength,7,arraylen])
ChlyDumby = np.zeros([interplength,arraylen])
for i in range(0,arraylen):
	dataout, rows, cols = csvread(fullnames[i])
	#
	TemInterp = interpolate.interp1d(dataout[0,:], dataout[1,:],kind='linear')
	DenInterp = interpolate.interp1d(dataout[0,:], dataout[2,:],kind='linear')
	SalInterp = interpolate.interp1d(dataout[0,:], dataout[3,:],kind='linear')
	PreInterp = interpolate.interp1d(dataout[0,:], dataout[4,:],kind='linear')
	OxyInterp = interpolate.interp1d(dataout[0,:], dataout[7,:],kind='linear')
	ChlInterp = interpolate.interp1d(dataout[0,:], dataout[11,:],kind='linear')
	BskInterp = interpolate.interp1d(dataout[0,:], dataout[12,:],kind='linear')
	CDMInterp = interpolate.interp1d(dataout[0,:], dataout[13,:],kind='linear')
	#
	for j in range(0,interplength):
		GridData[j,0,i] = TemInterp(Ygrid[j])
		GridData[j,1,i] = DenInterp(Ygrid[j])
		GridData[j,2,i] = SalInterp(Ygrid[j])
		GridData[j,3,i] = OxyInterp(Ygrid[j])			 
		ChlyDumby[j,i] = ChlInterp(Ygrid[j])	
		GridData[j,5,i] = BskInterp(Ygrid[j])	
		GridData[j,6,i] = CDMInterp(Ygrid[j])	
		if (ChlInterp(Ygrid[j]) < 0):
			GridData[j,4,i] = GridData[j-1,4,i]
		else: 
			GridData[j,4,i] = ChlyDumby[j,i]	
#
# ===========================================================
# Initial Data Contour Plotting
# ===========================================================
contoursdirpath ='C:/Documents and Settings/Alex/Documents/Research/Delaware/MUSTACHE/Jan01_Jun04_2013/plots/contours/'
#
#pp1 = PdfPages('LinearInterpolationGriddedContourDataPlots.pdf')
#
# Do you want to plot contours?  Yes: contplt = 1, No: else
contplt = 1
if contplt == 1:
	for k in range (0,7):
		#
		# Plotting correct data
		contourdat = np.zeros([interplength,arraylen])
		for ci in range(0,interplength):
			for cj in range(0,arraylen):
				#cjj = cj + 1
				contourdat[ci,cj] = GridData[ci,k,cj]
		#
		# Plot Set-up
		fig = plt.figure()
		ax = fig.add_subplot(1, 1, 1)
		#
		# Which Plotting Information?
		if k == 0: #temperature
			contlevels= np.linspace(-1.5,3,100)
			conticks = [-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]
			figtit = 'Temperature (C)'
			savetit = "Temperature_Contour"
			cs = plt.contourf(floatdate,Ygrid,contourdat, levels = contlevels)
			plt.colorbar(cs, spacing='proportional', ticks = conticks)
		elif k == 1: # density
			contlevels= np.linspace(1026.75,1027.8,100)
			conticks = [1026.50, 1026.75, 1027.00, 1027.25, 1027.50, 1027.75, 1028.00]
			figtit = 'Density (kg/m^3)'
			savetit = "Density_Contour"
			cs = plt.contourf(floatdate,Ygrid,contourdat, levels = contlevels)
			plt.colorbar(cs, spacing='proportional', ticks = conticks)
		elif k == 2: # sal
			contlevels= np.linspace(33.5,34.75,100)
			conticks = [33.25, 33.50, 33.75, 34.00, 34.25, 34.50, 34.75, 35.00]
			figtit = 'Salinity (psu)'
			savetit = "Salinity_Contour"
			cs = plt.contourf(floatdate,Ygrid,contourdat, levels = contlevels)
			plt.colorbar(cs, spacing='proportional', ticks = conticks)
		elif k == 3: #oxygen
			contlevels= np.linspace(4,10,100)
			conticks = [4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
			figtit = 'Oxygen Concentration (mg/l)'
			savetit = "Oxygen_Contour"
			cs = plt.contourf(floatdate,Ygrid,contourdat, levels = contlevels)
			plt.colorbar(cs, spacing='proportional', ticks = conticks)
		elif k == 4: # chl
			contlevels= np.logspace(-2.5,.65,100)
			conticks = [0.01,0.25,0.50, 1.00, 1.50, 2.00, 3.00]
			figtit = 'Chlorophyll (micro gram/liter)'
			savetit = "Chlorophyll_Contour"
			cs = plt.contourf(floatdate,Ygrid,contourdat, levels= contlevels, norm = LogNorm())
			cbar = plt.colorbar(cs, spacing='proportional', norm = LogNorm())
			cbar.set_ticks(conticks)
			cbar.set_ticklabels(conticks)
		elif k == 5: #backscatter
			contlevels= np.logspace(-4,-3.3,100)
			conticks = [0.0001, 0.0002, 0.0002, 0.0003, 0.0004, 0.0004, 0.0005, 0.00075, 0.001, 0.0015]
			figtit = 'Backscatter Cross Section (m^-1 sr^-1)'
			savetit = "Backscatter_Contour"
			cs = plt.contourf(floatdate,Ygrid,contourdat, levels= contlevels, norm = LogNorm())
			cbar = plt.colorbar(cs, spacing='proportional', norm = LogNorm())
			cbar.set_ticks(conticks)
			cbar.set_ticklabels(conticks)
		elif k == 6: #CDOM
			contlevels= np.linspace(1,3,100)
			conticks = [1.5, 2.0, 2.5, 3.0]
			figtit = 'CDOM (ppb)'
			savetit = "CDOM_Contour"
			cs = plt.contourf(floatdate,Ygrid,contourdat, levels= contlevels)
			plt.colorbar(cs, spacing='proportional', ticks = conticks)
		#
		# Plotting
		ax = plt.gca()
		ax.set_xlabel("Days since 01/01/2013")
		ax.set_ylabel("Depth (m)")
		ax.set_ylim([-900,0])
		# plt.title(figtit)
		fig.suptitle(figtit) 
		plt.savefig(contoursdirpath + savetit+".png")

os.chdir('C:/Documents and Settings/Alex/Documents/Research/Delaware/MUSTACHE/Jan01_Jun04_2013/code/interpolate_profiles')