import netCDF4
import numpy as np
import matplotlib
matplotlib.use('Agg')
from matplotlib.path import Path
import math
import utm
from track_exception import OnLandException
from track_exception import NegativeDepthException
from track_exception import NoNeighborException
from track_exception import NoDatasetException

def distance(lat1, lon1, lat2,lon2):
    """ 
    Calculates both distance in km and bearing
    note: "origin" and "destintation" are tuples (lat,lon) 
    note: if user inputs lat & lon as degrees-minutes (ddmm.m), it will convert to decimal degrees (dd.dddd)
    """
    #lat1, lon1 = origin
    #lat2, lon2 = destination
    if lat1>1000:
        (lat1,lon1)=dm2dd(lat1,lon1) # this is the conversion from degrees-minutes to decimal degrees
        (lat2,lon2)=dm2dd(lat2,lon2)
        print 'converted to from ddmm to dd.ddd'
    radius = 6371 # km
    

    dlat = math.radians(lat2-lat1)
    dlon = math.radians(lon2-lon1)
    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \
        * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = radius * c
    
    '''def calcBearing(lat1, lon1, lat2, lon2):
       dLon = lon2 - lon1
       y = math.sin(dLon) * math.cos(lat2)
       x = math.cos(lat1) * math.sin(lat2) \
           - math.sin(lat1) * math.cos(lat2) * math.cos(dLon)
       return math.atan2(y, x)
       
    bear= math.degrees(calcBearing(lat1, lon1, lat2, lon2))'''  
    return d

def dm2dd(lat,lon):
    """
    convert lat, lon from decimal degrees,minutes to decimal degrees
    """
    (a,b)=divmod(float(lat),100.)   
    aa=int(a)
    bb=float(b)
    lat_value=aa+bb/60.

    if float(lon)<0:
        (c,d)=divmod(abs(float(lon)),100.)
        cc=int(c)
        dd=float(d)
        lon_value=cc+(dd/60.)
        lon_value=-lon_value
    else:
        (c,d)=divmod(float(lon),100.)
        cc=int(c)
        dd=float(d)
        lon_value=cc+(dd/60.)
    return lat_value, -lon_value

def get_nc_data(url, *args):
    '''
    get specific dataset from url

    *args: dataset name, composed by strings
    ----------------------------------------
    example:
        url = 'http://www.nefsc.noaa.gov/drifter/drift_tcs_2013_1.dat'
        data = get_url_data(url, 'u', 'v')
    '''
    try:
        nc = netCDF4.Dataset(url)
    except:
        raise NoDatasetException({"message":"No such dataset."})
    # nc = netCDF4.Dataset('test.nc')
    data = {}
    for arg in args:
        try:
            data[arg] = nc.variables[arg]
        except (IndexError, NameError, KeyError):
            print 'Dataset {0} is not found'.format(arg)
        #print 'arg %s' % arg
        #print data[arg]
    return data

def nearest_point(lon, lat, lons, lats, length=0.06, k=1):  #0.3/5==0.06
    '''Find the nearest point to (lon,lat) from (lons,lats),
       return the nearest-point (lon,lat)
       author: Bingwei'''
    p = Path.circle((lon,lat),radius=length)
    #numpy.vstack(tup):Stack arrays in sequence vertically
    points = np.vstack((lons.flatten(),lats.flatten())).T  
    
    insidep = []
    #collect the points included in Path.
    for i in xrange(len(points)):
        if p.contains_point(points[i]):# .contains_point return 0 or 1
            insidep.append(points[i])  
    # if insidep is null, there is no point in the path.
    if not insidep:
        print 'There is no model-point near the given-point.'
        raise NoNeighborException({"message": "There is no model-point near the given-point."})
    #calculate the distance of every points in insidep to (lon,lat)
    distancelist = []
    for i in insidep:
        ss=math.sqrt((lon-i[0])**2+(lat-i[1])**2)
        distancelist.append(ss)
    # find index of the min-distance
    npdistancelist = np.array(distancelist)
    try:
        mindex = np.argpartition(npdistancelist, k)
    except Exception as err:
        print 'There is no model-point near the given-point.'
        raise NoNeighborException({"message": "There is no model-point near the given-point."})
    # location the point
    # lonp = insidep[mindex][0]; latp = insidep[mindex][1]
    
    # return lonp,latp
    npinsidep = np.array(insidep)
    return npinsidep[mindex[:k]]

def idw(values, ds):
    denominator = 0.
    numerator = 0.
    for i in range(len(values)):
        d = float(ds[i])
        v = float(values[i])
        denominator = denominator + (1. / d**2)
        numerator = numerator + (v / d**2)
    return numerator/denominator

def shrink_data(lon,lat,lons,lats):
    lont = []; latt = []
    p = Path.circle((lon,lat),radius=0.6)
    pints = np.vstack((lons.flatten(),lats.flatten())).T
    for i in range(len(pints)):
        if p.contains_point(pints[i]):
            lont.append(pints[i][0])
            latt.append(pints[i][1])
    lonl=np.array(lont); latl=np.array(latt)#'''
    if not lont:
        print 'point position error! shrink_data'
        raise Exception('point position error! shrink_data.')
    return lonl,latl

class get_model():
    def __init__(self, url):
        self.get_data(url)
        
    def get_data(self, url):
        '''
        return the data needed.
        url is from get_roms.get_url(starttime, endtime)
        '''
        data = get_nc_data(url, 'lon_rho','lat_rho','mask_rho','u','v','mask_u','mask_v','h','s_rho','zeta')
        self.lon_rho = data['lon_rho'][:]; self.lat_rho = data['lat_rho'][:] 
        # self.lon_u,self.lat_u = data['lon_u'][:], data['lat_u'][:]
        # self.lon_v,self.lat_v = data['lon_v'][:], data['lat_v'][:]
        self.h = data['h'][:]; self.s_rho = data['s_rho'][:]
        # self.mask_u = data['mask_u'][:]; self.mask_v = data['mask_v'][:]
        self.mask_rho = data['mask_rho'][:]
        self.u = data['u']; self.v = data['v']; self.zeta = data['zeta']
        #return self.fmodtime, self.emodtime
    
    def velocities_on_point(self, th_hour, lon, lat, depth):
        '''
        Find velocities u, v on (lon, lat)
        '''
        indexr = np.where(self.lon_rho[0][:]==lon)
        indeyr = np.where(self.lat_rho[:,0]==lat)
        indexu = indexv = indexr
        indeyu = indeyv = indeyr

        if not self.mask_rho[(indeyu,indexu)]:
            print 'No u velocity.'
            raise OnLandException({"message":"The given-point is on the land."})

        # if not self.mask_v[(indeyv,indexv)]:
        #     print 'No v velocity'
        #     raise OnLandException({"message":"The given-point is on the land."})
                
        waterdepth = self.h[(indeyr,indexr)]+self.zeta[th_hour][(indeyr,indexr)][0]

        if waterdepth<(abs(depth)): 
            print 'This point is too shallow. %f less than %d meter.'% (waterdepth, abs(depth))
            raise NegativeDepthException({"message":"Depth of given-point is less than %d meter." % abs(depth)})
        depth_total = self.s_rho*waterdepth  
        layer = np.argmin(abs(depth_total+depth))

        u_t = self.u[th_hour,layer][(indeyu,indexu)][0][0] 
        v_t = self.v[th_hour,layer][(indeyv,indexv)][0][0]

        return u_t, v_t, layer, waterdepth

    def idw_interpolation(self, th_hour, lon, lat, lonrho, latrho, depth, radius):
        points = nearest_point(lon, lat, lonrho, latrho, radius, 4)

        u_ts = []
        v_ts = []
        layers = []
        waterdepths = []
        distances = []
        for p in points:
            p_u, p_v, p_layer, p_waterdepth = self.velocities_on_point(th_hour, p[0], p[1], depth) #p[0]: lon, p[1]: lat
            u_ts.append(p_u)
            v_ts.append(p_v)
            layers.append(p_layer)
            waterdepths.append(p_waterdepth)
            p_d = distance(lat, lon, p[1], p[0])
            distances.append(p_d)

        # print u_ts, v_ts, layers, waterdepths, distances

        u_t = idw(u_ts, distances)
        v_t = idw(v_ts, distances)
        layer = idw(layers, distances)
        waterdepth = idw(waterdepths, distances)

        return u_t, v_t, layer, waterdepth

    def rk4_interpolation(self, th_hour, lon, lat, lonrho, latrho, track_way, depth, radius):
        h = 60.*60. #3600 secs
        p_angle = 0.
        xpar, ypar, zNumber, zLetter = utm.from_latlon(lat, lon)

        # step 1
        u_ad, v_ad, layer_ad, waterdepth_ad = self.idw_interpolation(th_hour, lon, lat, lonrho, latrho, depth, radius)
        kn1_u = u_ad
        kn1_v = v_ad
        x1 = xpar + (u_ad * math.cos(p_angle) - v_ad * math.sin(p_angle)) * h / 2.
        y1 = ypar + (u_ad * math.sin(p_angle) + v_ad * math.cos(p_angle)) * h / 2.

        # step 2
        lat, lon = utm.to_latlon(x1, y1, zNumber, zLetter)
        u_ad, v_ad, layer_ad, waterdepth_ad = self.idw_interpolation(th_hour, lon, lat, lonrho, latrho, depth, radius)
        kn2_u = u_ad
        kn2_v = v_ad
        x2 = xpar + (u_ad * math.cos(p_angle) - v_ad * math.sin(p_angle)) * h / 2.
        y2 = ypar + (u_ad * math.sin(p_angle) + v_ad * math.cos(p_angle)) * h / 2.

        # step 3
        lat, lon = utm.to_latlon(x2, y2, zNumber, zLetter)
        u_ad, v_ad, layer_ad, waterdepth_ad = self.idw_interpolation(th_hour, lon, lat, lonrho, latrho, depth, radius)
        kn3_u = u_ad
        kn3_v = v_ad
        x3 = xpar + (u_ad * math.cos(p_angle) - v_ad * math.sin(p_angle)) * h / 2.
        y3 = ypar + (u_ad * math.sin(p_angle) + v_ad * math.cos(p_angle)) * h / 2.

        # step 4
        lat, lon = utm.to_latlon(x3, y3, zNumber, zLetter)
        u_ad, v_ad, layer_ad, waterdepth_ad = self.idw_interpolation(th_hour, lon, lat, lonrho, latrho, depth, radius)
        kn4_u = u_ad
        kn4_v = v_ad

        # rk formulat
        P_U = (kn1_u + 2. * kn2_u + 2. * kn3_u + kn4_u) / 6.
        P_V = (kn1_v + 2. * kn2_v + 2. * kn3_v + kn4_v) / 6.

        dx = h * (P_U * math.cos(p_angle) - P_V * math.sin(p_angle))
        dy = h * (P_U * math.sin(p_angle) + P_V * math.cos(p_angle))

        return dx, dy, 'NA', 'NA'

    def euler_interpolation(self, th_hour, lon, lat, lonrho, latrho, track_way, depth, radius):
        u_t, v_t, layer, waterdepth = self.idw_interpolation(th_hour, lon, lat, lonrho, latrho, depth, radius)
        # print u_t, v_t, layer, waterdepth

        if track_way=='backward': # backwards case
            u_t = -1*u_t
            v_t = -1*v_t
                    
        # print 'u_t,v_t',u_t,v_t
        if np.isnan(u_t) or np.isnan(v_t): #There is no water
            print 'Sorry, the point on the land or hits the land. Info: u or v is NAN'
            return nodes
        
        dx = 60*60*u_t#float(u_p)
        dy = 60*60*v_t#float(v_p)

        return dx, dy, layer, waterdepth

    def get_track_nbi_gps(self,lon,lat,depth,radius,track_way, start, end):#, depth
        '''
        get the nodes of specific time period
        lon, lat: start point
        depth: 0~35, the 0th is the bottom.
        '''
        # self.hours = 72
        nodes = dict(lon=[lon], lat=[lat], time=[start])
        if track_way == 'backward':
            step = -1
        else:
            step = 1
        lonrho,latrho = shrink_data(lon, lat, self.lon_rho, self.lat_rho)
        # t = abs(self.hours)
        for i in range(start, end, step):  #Roms points update every 2 hour
            # print 'hour: ', i, lat, lon
            # if i%72==0 :
                #print 'layer,lon,lat,i',layer,lon,lat,i
                # lonrho,latrho = self.shrink_data(lon,lat,self.lon_rho,self.lat_rho)
                # lonu,latu = self.shrink_data(lon,lat,self.lon_u,self.lat_u)
                # lonv,latv = self.shrink_data(lon,lat,self.lon_v,self.lat_v)

            try:
                p = nearest_point(lon, lat, lonrho, latrho, radius)
                u_t, v_t, layer, waterdepth = self.velocities_on_point(i, p[0][0], p[0][1], depth)

                if track_way=='backward': # backwards case
                    u_t = -1*u_t
                    v_t = -1*v_t
                    
                # print 'u_t,v_t',u_t,v_t
                if np.isnan(u_t) or np.isnan(v_t): #There is no water
                    print 'Sorry, the point on the land or hits the land. Info: u or v is NAN'
                    return nodes
                dx = 60*60*u_t#float(u_p)
                dy = 60*60*v_t#float(v_p)
                #mapx = Basemap(projection='ortho',lat_0=lat,lon_0=lon,resolution='l')                        
                #x,y = mapx(lon,lat)
                #lon,lat = mapx(x+dx,y+dy,inverse=True)            
                lon = lon + dx/(111111*np.cos(lat*np.pi/180))
                lat = lat + dy/111111
                if track_way == 'backward':
                    nodes['time'].append(i)
                    print '%d,lat,lon,layer'%(i),lat,lon,layer,waterdepth
                else:
                    nodes['time'].append(i+1)
                    print '%d,lat,lon,layer'%(i+1),lat,lon,layer,waterdepth
                nodes['lon'].append(lon);nodes['lat'].append(lat)

            except OnLandException as err:
                raise err
            except NoNeighborException as err:
                if len(nodes['lon'])==1:
                    raise err
                else:
                    return nodes
            except NegativeDepthException as err:
                if len(nodes['lon'])==1:
                    raise err
                else:
                    return nodes
            except Exception as err:
                print err
                return nodes
        return nodes

    def get_track_idw_utm(self,lon,lat,depth,radius,track_method,track_way, start, end):
        nodes = dict(lon=[lon], lat=[lat], time=[start])
        if track_way == 'backward':
            step = -1
        else:
            step = 1
        lonrho,latrho = shrink_data(lon, lat, self.lon_rho, self.lat_rho)
        for i in range(start, end, step):
            try:
                if track_method == 'euler':
                    dx, dy, layer, waterdepth = self.euler_interpolation(i, lon, lat, lonrho, latrho, track_way, depth, radius)
                elif track_method == 'rk4':
                    dx, dy, layer, waterdepth = self.rk4_interpolation(i, lon, lat, lonrho, latrho, track_way, depth, radius)
                x, y, zNumber, zLetter = utm.from_latlon(lat, lon)
                lat, lon = utm.to_latlon(x+dx, y+dy, zNumber, zLetter)
                if track_way == 'backward':
                    nodes['time'].append(i)
                    print '%d,lat,lon,layer'%(i),lat,lon,layer,waterdepth
                else:
                    nodes['time'].append(i+1)
                    print '%d,lat,lon,layer'%(i+1),lat,lon,layer,waterdepth
                nodes['lon'].append(lon);nodes['lat'].append(lat)

            except OnLandException as err:
                raise err
            except NoNeighborException as err:
                if len(nodes['lon'])==1:
                    raise err
                else:
                    return nodes
            except NegativeDepthException as err:
                if len(nodes['lon'])==1:
                    raise err
                else:
                    return nodes
            except Exception as err:
                print err
                return nodes
        return nodes