#! /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
#
#			2) Calculates vertical velocities for density
#				and chlorphyll. 
#
#			3) Plots the vertical velocities over contoured 
#				gridded, interpolated w data
#				
#
# ===========================================================
# 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]

#
# ===========================================================
# Read the Data and Interpolate for Gridding
# ===========================================================
#
Ygrid = np.linspace(-1800,-10,1791*2)
interplength = len(Ygrid)
GridGrad = np.zeros([interplength,1,arraylen])
for i in range(0,arraylen):
	gradout, grows, gcols = csvread(gradnames[i])
	#
	drhodzInterp = interpolate.interp1d(gradout[0,:], gradout[6,:],kind='linear')
	#
	for j in range(0,interplength):
		GridGrad[j,0,i] = drhodzInterp(Ygrid[j])
#
#Verification Ploting 
fig = plt.figure()
ax1 = fig.add_subplot(1,1,1)
#p11 = ax1.plot(Ygrid,GridData[:,5,i],'k')
#p12 = ax1.plot(Ygrid,ChlyDumby[:,i])
#p13 = ax1.scatter(dataout[0,:], dataout[11,:])
#ax1.set_xlabel("Depth(m)")
#ax1.set_ylabel("Chlorophyll Concentration")

#pp11 = ax1.plot(Ygrid,GridData[:,5,i],'k')
pp1 = ax1.scatter(np.abs(GridGrad[:,0,:]),GridData[:,4,:])
ax1.set_xlabel("d(rho)/dz (kg/m^4)")
ax1.set_ylabel("Chlorophyll Flouresence (micro grams/liter") 
plt.savefig("C:/Documents and Settings/Alex/Documents/Research/Delaware/MUSTACHE/Jan01_Jun04_2013/plots/drhodz_v_chly.png")
#
os.chdir('C:/Documents and Settings/Alex/Documents/Research/Delaware/MUSTACHE/Jan01_Jun04_2013/code/interpolate_profiles')