#! /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: 01/13/2014
#
# NOTES:    1) Does all the OSCAR Data Processing
#
#               
#
# ===========================================================
# Import Modules
# ===========================================================
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from datetime import datetime
from mpl_toolkits.basemap import Basemap, addcyclic
from mpl_toolkits.basemap import shiftgrid, cm
from scipy.ndimage.filters import minimum_filter, maximum_filter
from netCDF4 import Dataset
import glob
import os
from scipy.io import netcdf
import matplotlib.transforms as mtransforms
from matplotlib.patches import FancyBboxPatch
# from matplotlib.transforms import Bboxcx
from matplotlib.path import Path
import matplotlib.patches as patches
import math as ma
import csv
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
import scipy
from scipy import stats
#
# ===========================================================
#
# FUNCTIONS
#
# ===========================================================
#
# ===========================================================
# Function to input delX, delY, return bearing and magnitude
# ===========================================================
#
def find_bearing_magnitude(X, Y):
    import math as m
    if ((Y >= 0) and (X >=0 )):
        bearing_rads = m.atan(abs(X)/abs(Y))
        bearing = m.degrees(bearing_rads)
        magnitude = abs(X)/m.sin(bearing_rads)

    elif ((Y < 0) and (X >=0)):
        bearing_rads = m.atan(abs(Y)/abs(X))
        bearing = m.degrees(bearing_rads) + 90
        magnitude = abs(Y)/m.sin(bearing_rads)

    elif ((Y < 0) and (X < 0)):
        bearing_rads = m.atan(abs(X)/abs(Y))
        bearing = m.degrees(bearing_rads) + 180
        magnitude = abs(X)/m.sin(bearing_rads)

    elif ((Y >= 0) and (X < 0)):
        bearing_rads = m.atan(abs(Y)/abs(X))
        bearing = m.degrees(bearing_rads) + 270
        magnitude = abs(Y)/m.sin(bearing_rads)
    else:
        print "#### ERROR IN: FIND BEARING AND MAG FUNCTION ####"

    return bearing_rads, bearing, magnitude
#
# ===========================================================
# Function Determine new Lat/Lon from Old + Bearing + dist
# ===========================================================
#
def destpoint(lat1, lon1, brng, d):
    import math

    R = 6374 #Radius of the Earth
    # brng = Bearing converted to radians.
    # d = Distance in km
    # lat1 = Current lat point converted to radians
    # lon1 = Current lon point converted to radians

    lat2 = math.asin(math.sin(lat1)*math.cos(d/R) + math.cos(lat1)*math.sin(d/R)*math.cos(brng))
    lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1), math.cos(d/R)-math.sin(lat1)*math.sin(lat2))
    
    lat22 = math.degrees(lat2)
    lon22 = math.degrees(lon2)

    return lat22, lon22
#
# ===========================================================
# 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
def pos2dist(pos):
    # Pos = [lon1,lat1,lon2,lat2]
    import math as m
    R_aver = 6374
    pi = 3.14
    deg2rad = pi/180

    lat1 = pos[1] * deg2rad #lat1
    lon1 = pos[0] * deg2rad # lon1
    lat2 = pos[3] * deg2rad #lat2
    lon2 = pos[2] * deg2rad  #lon2
    dlon = lon2 - lon1
    dlat = lat2 - lat1

    e = m.cos(lat1)*m.cos(lat2)*m.sin(dlon/2)**2
    d = m.sin(dlat/2)**2
    c = d + e
    b = c**0.5
    dist = 2*R_aver*m.asin(b)

    # Returns a distance in KM
    return dist
#
# ===========================================================
# Function to new lon from zonal displacemet in km
# ===========================================================
#
def lonchange(input):
    # input = [lat,lon1,dzon]
    import math as m
    R_aver = 6374
    pi = 3.14
    deg2rad = pi/180
    rad2deg = 180/pi

    lat = input[0] * deg2rad #lat1
    lon1 = input[1] * deg2rad # lon1
    dzon = input[2] # zonal displacement in km2

    c = m.cos(lat)
    b = m.sin(dzon/2/R_aver)
    a = b/c
    lon2rad = 2*m.asin(a) + lon1
    lon2 = lon2rad*rad2deg
    return lon2
#
# ===========================================================
# Function to new lat from meridional displacemet in km
# ===========================================================
#
def latchange(input):
    # input = [lat1,dmer]
    import math as m
    R_aver = 6374
    pi = 3.14
    deg2rad = pi/180
    rad2deg = 180/pi


    lat1 = input[0]*deg2rad
    dmer = input[1] # zonal displacement in km2


    lat2rad = dmer/R_aver + lat1
    lat2 = lat2rad*rad2deg

    return lat2
#
# ===========================================================
# Function to Do Regressions
# ===========================================================
#
def linearregress(x, y):
    from scipy import stats
    #
    # Linearly Regression
    slope, intercept, rvalue, pvalue, stderr = stats.linregress(x,y)
    #
    # Regression Line
    regressline = slope*x+intercept
    #
    return regressline, rvalue, pvalue, stderr
#
#
# ===================================================================================
# Function for Bilinear Interpolation
# ===================================================================================
#
def bilinear_interpolation(x, y, points):
    #
    # Order Points vy x, then by y
    points = sorted(points)               # order points by x, then by y
    (x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points
    #
    # Are the points a rectangle?
    if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2:
        raise ValueError('points do not form a rectangle')
    if not x1 <= x <= x2 or not y1 <= y <= y2:
        raise ValueError('(x, y) not within the rectangle')
    #
    # Bilinear Interp and return value.
    return (q11 * (x2 - x) * (y2 - y) +
            q21 * (x - x1) * (y2 - y) +
            q12 * (x2 - x) * (y - y1) +
            q22 * (x - x1) * (y - y1)
           ) / ((x2 - x1) * (y2 - y1) + 0.0)




# open the netcdf file to read
ncf = netcdf.netcdf_file('/home/ardavies/satdata/OSCAR/sofex2002_north.nc','r')
#
# Getting the netcdf data into the corresponding arrays
datadates = ncf.variables['time'][:].squeeze()
v = ncf.variables['v'][:].squeeze()
u = ncf.variables['u'][:].squeeze()
rawlon = ncf.variables['lon'][:].squeeze() # degrees E
rawlat = ncf.variables['lat'][:].squeeze() # degrees N
#
# Array Initailization
numdates = len(u[:,0,0])
lonlen =  len(rawlon)
latlen =  len(rawlat)
lat = np.zeros([latlen,lonlen])
lon = np.zeros([latlen,lonlen])
uavg = np.zeros([numdates-1,latlen,lonlen])
vavg = np.zeros([numdates-1,latlen,lonlen])
#
# Getting lat/lon's in i,j arrays
for i in range(0,lonlen):
    lon[:,i] = rawlon[i]
for j in range(0,latlen):
    lat[j,:] = rawlat[j]
#
# DOY of 5 day avg Oscar data
centerdates = np.zeros(numdates)
for k in range(0,numdates):
    centerdates[k] = datadates[k]

import math as ma
mke5day = np.zeros([numdates,latlen,lonlen])
#
# MKE from Raw OSCAR Data
for k in range(0,numdates):
    for i in range(0,lonlen):
        for j in range(0,latlen):
            mke5day[k,j,i] = 0.5*((u[k,j,i]*100)**2 + (v[k,j,i]*100)**2)
print rawlon
print rawlat
print mke5day[:,1,4]

#
# ===================================================================================
# Find Nearest Lat's and Lon's to float location for U and V components
# ===================================================================================
#

print u[:,0,0]
sitelat = -56.23
sitelon = 188
arraylen = numdates
#
for jj in range(0,arraylen):
    #
    # Get the closest OSCAR data latitude with respect to the float latitude
    latvalue1, latindex1 = find_nearest(lat[:,0], sitelat)
    #
    # Find the next closest OSCAR data latitude with respect to the float latitude
    upperlat = lat[latindex1+1,0]
    lowerlat = lat[latindex1-1,0]
    #
    # Closest lat + 0.33 degrees?
    if (abs(sitelat-upperlat)<=abs(sitelat-lowerlat)):
        latindex2= latindex1+ 1
    #
    # Closest lat - 0.33 degrees?
    else:
        latindex2 = latindex1 - 1
    latvalue2 = lat[latindex2,0]
    #
    # Get the closest OSCAR data longitude with respect to the float longitude
    lonvalue1, lonindex1 = find_nearest(lon[0,:], sitelon)
    #
    # Find the next closest OSCAR data longitude with respect to the float longitude
    upperlon = lon[0,lonindex1+1]
    lowerlon = lon[0,lonindex1-1]
    #
    # Closest lon + 0.33 degrees?
    if (abs(sitelon-upperlon)<=abs(sitelon-lowerlon)):
        lonindex2 = lonindex1 + 1
    #
    # Closest lon - 0.33 degrees?
    else:
        lonindex2 = lonindex1 - 1
    lonvalue2 = lon[0,lonindex2]

print latvalue1
print latvalue2
print lonvalue1
print lonvalue2
#
# ===================================================================================
# Original 5 Day Bi-Linear Interpolation
# ===================================================================================
#
# Initializing Arrays for Interpolation
siteMKE5day = np.zeros(arraylen)
#
# Loop through the float data
for k in range(0,arraylen):
    #
    # Interpolate Linearly Along Closest Latitude (latvalue1)
    x1_closest5day  = abs(pos2dist([lonvalue1,latvalue1,sitelon,latvalue1]))
    x2_closest5day  = abs(pos2dist([lonvalue2,latvalue1,sitelon,latvalue1]))

    xtot_closest5day  = abs(pos2dist([lonvalue1,latvalue1,lonvalue2,latvalue1]))
    MKEinterp_closest5day  = mke5day[k,latindex1, lonindex1]*(x2_closest5day /xtot_closest5day ) + mke5day[k, latindex1, lonindex2]*(x1_closest5day /xtot_closest5day)
    #
    # Interpolate Linearly Along Furthest Latitude (latvalue2)
    x1_furthest5day  = abs(pos2dist([lonvalue1,latvalue2,sitelon,latvalue2]))
    x2_furthest5day  = abs(pos2dist([lonvalue2,latvalue2,sitelon,latvalue2]))
    
    xtot_furthest5day  = abs(pos2dist([lonvalue1,latvalue2,lonvalue2,latvalue2]))
    MKEinterp_furthest5day  = mke5day[k,latindex2, lonindex1]*(x2_furthest5day /xtot_furthest5day ) + mke5day[k, latindex2, lonindex2]*(x1_furthest5day /xtot_furthest5day)
    #
    # Tnterpolate Along Londitude
    y_closest5day  = abs(pos2dist([sitelon,latvalue1,sitelon,sitelat]))
    y_furthest5day  = abs(pos2dist([sitelon,latvalue2,sitelon,sitelat]))
    ytot5day  = abs(pos2dist([sitelon,latvalue1,sitelon,latvalue2]))
    siteMKE5day[k] = MKEinterp_furthest5day *(y_closest5day /ytot5day ) + MKEinterp_closest5day *(y_furthest5day /ytot5day )

print siteMKE5day
print datadates