import sys
import netCDF4
#import calendar
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from datetime import datetime,timedelta
import numpy as np
import pandas as pd
from dateutil.parser import parse
#import pytz
from matplotlib.path import Path
import math
from mpl_toolkits.basemap import Basemap
import colorsys
from sympy import *
from sympy.geometry import *
from fractions import Fraction
from track_exception import OnLandException as OnLandException
from track_exception import NegativeDepthException as NegativeDepthException
from track_exception import NoNeighborException as NoNeighborException
from track_exception import NoDatasetException as NoDatasetException

def distance(lat1, lon1, lat2,lon2):
    """ 
    Calculates both distance 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 getrawdrift(did,filename):
    '''
    routine to get raw drifter data from ascii files posted on the web
    '''
    url='http://nefsc.noaa.gov/drifter/'+filename
    df=pd.read_csv(url,header=None, delimiter="\s+")
    # make a datetime
    dtime=[]
    index = np.where(df[0]==int(did))[0]
    newData = df.ix[index]
    for k in newData[0].index:
        dt1=datetime(2015, newData[2][k],newData[3][k],newData[4][k],newData[5][k])
        dtime.append(dt1)
    #print dtime
    return newData[8],newData[7],dtime,newData[9] # lat,lon,time,

def getdrift(did):
    """
    routine to get drifter data from archive based on drifter id (did)
    -assumes "import pandas as pd" has been issued above
    -get remotely-stored drifter data via ERDDAP
    -input: deployment id ("did") number where "did" is a string
    -output: time(datetime), lat (decimal degrees), lon (decimal degrees), depth (meters)
    
    note: there is another function below called "data_extracted" that does a similar thing returning a dictionary
    
    Jim Manning June 2014
    """
    url = 'http://comet.nefsc.noaa.gov:8080/erddap/tabledap/drifters.csv?time,latitude,longitude,depth&id="'+did+'"&orderBy("time")'
    df=pd.read_csv(url,skiprows=[1]) #returns a dataframe with all that requested
    #print df    
    # generate this datetime 
    for k in range(len(df)):
       df.time[k]=parse(df.time[k]) # note this "parse" routine magically converts ERDDAP time to Python datetime
    return df.latitude.values,df.longitude.values,df.time.values,df.depth.values  

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 input_with_default(data, v_default):
    '''
    data: string, could be name of value you want to get
    v_default
    '''
    l = (data, str(v_default))
    try:
        data_input = input('Please input %s(default %s)(If don\'t want to make change, press "Enter"): ' % l)
    except SyntaxError:
        data_output = v_default
    else:
        data_output = data_input
    return data_output
    
def shrink(a,b):
    """Return array shrunk to fit a specified shape by triming or averaging.
    
    a = shrink(array, shape)
    
    array is an numpy ndarray, and shape is a tuple (e.g., from
    array.shape). a is the input array shrunk such that its maximum
    dimensions are given by shape. If shape has more dimensions than
    array, the last dimensions of shape are fit.
    
    as, bs = shrink(a, b)
    
    If the second argument is also an array, both a and b are shrunk to
    the dimensions of each other. The input arrays must have the same
    number of dimensions, and the resulting arrays will have the same
    shape.
    Example
    -------
    
    >>> shrink(rand(10, 10), (5, 9, 18)).shape
    (9, 10)
    >>> map(shape, shrink(rand(10, 10, 10), rand(5, 9, 18)))        
    [(5, 9, 10), (5, 9, 10)]   
       
    """

    if isinstance(b, np.ndarray):
        if not len(a.shape) == len(b.shape):
            raise Exception, \
                  'input arrays must have the same number of dimensions'
        a = shrink(a,b.shape)
        b = shrink(b,a.shape)
        return (a, b)

    if isinstance(b, int):
        b = (b,)

    if len(a.shape) == 1:                # 1D array is a special case
        dim = b[-1]
        while a.shape[0] > dim:          # only shrink a
#            if (dim - a.shape[0]) >= 2:  # trim off edges evenly
            if (a.shape[0] - dim) >= 2:
                a = a[1:-1]
            else:                        # or average adjacent cells
                a = 0.5*(a[1:] + a[:-1])
    else:
        for dim_idx in range(-(len(a.shape)),0):
            dim = b[dim_idx]
            a = a.swapaxes(0,dim_idx)        # put working dim first
            while a.shape[0] > dim:          # only shrink a
                if (a.shape[0] - dim) >= 2:  # trim off edges evenly
                    a = a[1:-1,:]
                if (a.shape[0] - dim) == 1:  # or average adjacent cells
                    a = 0.5*(a[1:,:] + a[:-1,:])
            a = a.swapaxes(0,dim_idx)        # swap working dim back
    return a

def data_extracted(filename,drifter_id=None,starttime=None):
    '''
    get a dictionary called "data" made of time, lon, lat from local file.
    filename: local file diretory
    drifter_id: the specific data of some id you want.
    starttime: have to be input with drifter_id, or just drifter_id.
    '''
    data = {}
    did, dtime, dlon, dlat = [], [], [], []
    with open(filename, 'r') as f:
        for line in f.readlines():
            try:
                line = line.split()
                did.append(int(line[0]))
                dtime.append(datetime(year=2013,
                                      month=int(line[2]),day=int(line[3]),
                                      hour=int(line[4]),minute=int(line[5])))
                dlon.append(float(line[7]))
                dlat.append(float(line[8]))
            except IndexError:
                continue
    if drifter_id is not None:
        i = index_of_value(did, drifter_id)
        if starttime is not None:
            dtime_temp = dtime[i[0]:i[-1]+1]
            j = index_of_value(dtime_temp, starttime)
            data['time'] = dtime[i[0]:i[-1]+1][j[0]:]
            data['lon'] = dlon[i[0]:i[-1]+1][j[0]:]
            data['lat'] = dlat[i[0]:i[-1]+1][j[0]:]
        else:
            data['time'] = dtime[i[0]:i[-1]+1]
            data['lon'] = dlon[i[0]:i[-1]+1]
            data['lat'] = dlat[i[0]:i[-1]+1]
    elif drifter_id is None and starttime is None:
        data['time'] = dtime
        data['lon'] = dlon
        data['lat'] = dlat
    else:
        raise ValueError("Please input drifter_id while starttime is input")
    return data

def index_of_value(dlist,dvalue):
    '''
    return the indices of dlist that equals dvalue
    '''
    index = []
    startindex = dlist.index(dvalue)
    i = startindex
    for v in dlist[startindex:]:
        if v == dvalue:
            index.append(i)
        i+=1
    return index

class track(object):
    def __init__(self, startpoint):
        '''
        gets the start point of the water, and the location of datafile.
        '''
        self.startpoint = startpoint
        
    def get_data(self, url):
        '''
        calls get_data
        '''        
        pass                                 
        
    def nearest_point_index(self, lon, lat, lons, lats,rad):  #,num=4
        '''
        Return the nearest point(lonp,latp) and distance to origin point(lon,lat).
        lon, lat: the coordinate of start point, float
        latp, lonp: the coordinate of points to be calculated.
        '''
        def bbox2ij(lon, lat, lons, lats, rad):  
            
            '''bbox = [lon-length, lon+length, lat-length, lat+length]
            bbox = np.array(bbox)
            mypath = np.array([bbox[[0,1,1,0]],bbox[[2,2,3,3]]]).T#'''
            p = Path.circle((lon,lat),radius=rad)
            points = np.vstack((lons,lats)).T  #numpy.vstack(tup):Stack arrays in sequence vertically
            #print 'lons',lons
            #tshape = np.shape(lons)
            
            inside = []
            
            for i in range(len(points)):
                inside.append(p.contains_point(points[i]))  # .contains_point return 0 or 1
                
            sidex = np.array(inside, dtype=bool)#.reshape(tshape)
            index = np.where(sidex==True)
            
            '''check if there are no points inside the given area'''        
            
            if not index[0].tolist():          # bbox covers no area
                print 'This point is out of the model area.'
                raise Exception()
                
            else:
                return index        
        
        def min_distance(lon,lat,lons,lats):
            '''Find out the nearest distance to (lon,lat),and return lon.distance units: meters'''
            #mapx = Basemap(projection='ortho',lat_0=lat,lon_0=lon,resolution='l')
            dis_set = []
            #x,y = mapx(lon,lat)
            for i,j in zip(lons,lats):
                #x2,y2 = mapx(i,j)
                ss=math.sqrt((lon-i)**2+(lat-j)**2)
                #ss=math.sqrt((x-x2)**2+(y-y2)**2)
                dis_set.append(ss)
            dis = min(dis_set)
            p = dis_set.index(dis)
            lonp = lons[p]; latp = lats[p]
            return lonp,latp,dis       
        index = bbox2ij(lon, lat, lons, lats,rad)
        lon_covered = lons[index];  lat_covered = lats[index]       
        lonp,latp,distance = min_distance(lon,lat,lon_covered,lat_covered)
        #index1 = np.where(lons==lonp)
        #index2 = np.where(lats==latp)
        #index = np.intersect1d(index1,index2)        
        #points = np.vstack((lons.flatten(),lats.flatten())).T 
        #index = [i for i in xrange(len(points)) if ([lonp,latp]==points[i]).all()]
        '''index = []
        for i in len(points):
            if np.all([[lonp,latp],points[i]]):
                index.append(i)'''
                
        #index = np.where(points==[lonp,latp])
        #print 'index',index
        return lonp,latp,distance  #,lonp,latp
        
    def get_track(self, timeperiod, data):
        pass
    
class get_timcom(track):
    def __init__(self):
        pass
    
    def nearest_point(self, lon, lat, lons, lats, length=0.06):  #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
        mindex = np.argmin(distancelist)
        # location the point
        lonp = insidep[mindex][0]; latp = insidep[mindex][1]
        
        return lonp,latp

    def shrink_data(self,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
        
    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 get_track(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 = self.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:
                lonrp,latrp = self.nearest_point(lon,lat,lonrho,latrho,radius)
                # lonup,latup = self.nearest_point(lon,lat,lonu,latu)
                # lonvp,latvp = self.nearest_point(lon,lat,lonv,latv)

                # indexu = np.where(self.lon_u[0][:]==lonup)
                # indeyu = np.where(self.lat_u[:,0]==latup)
                # indexv = np.where(self.lon_v[0][:]==lonvp)
                # indeyv = np.where(self.lat_v[:,0]==latvp)
                indexr = np.where(self.lon_rho[0][:]==lonrp)
                indeyr = np.where(self.lat_rho[:,0]==latrp)
                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."})
                
                if track_way=='backward': # backwards case
                    waterdepth = self.h[(indeyr,indexr)]+self.zeta[i][(indeyr,indexr)][0]
                    # nodes['time'].append(self.mTime[t-i-1])
                else:
                    waterdepth = self.h[(indeyr,indexr)]+self.zeta[i][(indeyr,indexr)][0]
                    #nodes['time'].append(self.mTime[i+1])

                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))
            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

            if track_way=='backward': # backwards case
                u_t = -1*self.u[i,layer][(indeyu,indexu)][0][0] 
                v_t = -1*self.v[i,layer][(indeyv,indexv)][0][0]
            else:
                u_t = self.u[i,layer][(indeyu,indexu)][0][0] 
                v_t = self.v[i,layer][(indeyv,indexv)][0][0]
            # 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)
            
            
        return nodes

class get_roms(track):
    '''
    ####(2009.10.11, 2013.05.19):version1(old) 2009-2013
    ####(2013.05.19, present): version2(new) 2013-present
    (2006.01.01 01:00, 2014.1.1 00:00)
    '''
    
    def __init__(self):
        pass
    
    def nearest_point(self, lon, lat, lons, lats, length=0.06):  #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
        mindex = np.argmin(distancelist)
        # location the point
        lonp = insidep[mindex][0]; latp = insidep[mindex][1]
        
        return lonp,latp
        
    def get_url(self, starttime, endtime):
        print '''
        roms get url according to starttime and endtime.
        '''
        starttime = starttime
        self.hours = int((endtime-starttime).total_seconds()/60/60) # get total hours
        # time_r = datetime(year=2006,month=1,day=9,hour=1,minute=0)
        print self.hours
        sys.exit()
        url_oceantime = '''http://tds.marine.rutgers.edu:8080/thredds/dodsC/roms/espresso/2013_da/his/ESPRESSO_Real-Time_v2_History_Best?time'''
        url = """http://tds.marine.rutgers.edu:8080/thredds/dodsC/roms/espresso/2013_da/his/ESPRESSO_Real-Time_v2_History_Best?h[0:1:81][0:1:129],
        mask_rho[0:1:81][0:1:128],mask_u[0:1:81][0:1:128],mask_v[0:1:80][0:1:129],zeta[{0}:1:{1}][0:1:81][0:1:129],u[{0}:1:{1}][0:1:35][0:1:81][0:1:128],
        v[{0}:1:{1}][0:1:35][0:1:80][0:1:129],s_rho[0:1:35],lon_rho[0:1:81][0:1:129],lat_rho[0:1:81][0:1:129],lon_u[0:1:81][0:1:128],lat_u[0:1:81][0:1:128],
        lon_v[0:1:80][0:1:129],lat_v[0:1:80][0:1:129],time[0:1:19523]"""
        try:
            oceantime = netCDF4.Dataset(url_oceantime).variables['time'][:]
        except:
            print 'ROMS database is unavailable!'
            raise Exception()
        # get model works time horizon(UTC).
        fmodtime = datetime(2013,05,18) + timedelta(hours=float(oceantime[0]))
        emodtime = datetime(2013,05,18) + timedelta(hours=float(oceantime[-1]))
        mstt = fmodtime.strftime('%m/%d/%Y %H:%M') #model start time
        mett = emodtime.strftime('%m/%d/%Y %H:%M') #model end time
        # get number of hour from 05/18/2013
        t1 = (starttime - datetime(2013,05,18)).total_seconds()/3600 
        t2 = (endtime - datetime(2013,05,18)).total_seconds()/3600
        #t1 = int(round(t1)); t2 = int(round(t2))
        # judge if the starttime and endtime in the model time horizon
        print starttime, endtime
        if starttime<fmodtime or starttime>emodtime or endtime<fmodtime or endtime>emodtime:
            print 'Time: Error! Model(ROMS) only works between %s with %s.'%(mstt,mett)
            raise Exception()
        #index1 = np.where(oceantime==t1)[0][0]; #print index1
        #index2 = np.where(oceantime==t2)[0][0]; #print index2
        int1 = oceantime - t1; int2 = oceantime - t2
        index1 = np.argmin(abs(int1)); index2 = np.argmin(abs(int2))
        url = url.format(index1, index2)
        Times = []
        for i in range(self.hours+1):
            Times.append(fmodtime+timedelta(hours=i))
        self.mTime = Times
        self.url=url
        print url
        return url
    
    def shrink_data(self,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
        
    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','lon_u','lat_u','lon_v','lat_v','mask_rho','mask_u','mask_v','u','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'][:]#; mask_rho = data['mask_rho'][:]
        self.u = data['u']; self.v = data['v']; self.zeta = data['zeta']
        #return self.fmodtime, self.emodtime
        
    def get_track(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 = 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)

        # 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:
                lonrp,latrp = self.nearest_point(lon,lat,lonrho,latrho,radius)
                lonup,latup = self.nearest_point(lon,lat,lonu,latu,radius)
                lonvp,latvp = self.nearest_point(lon,lat,lonv,latv,radius)

                indexu = np.where(self.lon_u[0][:]==lonup)
                indeyu = np.where(self.lat_u[:,0]==latup)
                indexv = np.where(self.lon_v[0][:]==lonvp)
                indeyv = np.where(self.lat_v[:,0]==latvp)
                indexr = np.where(self.lon_rho[0][:]==lonrp)
                indeyr = np.where(self.lat_rho[:,0]==latrp)

                if not self.mask_u[(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."})
                
                if track_way=='backward': # backwards case
                    waterdepth = self.h[(indeyr,indexr)]+self.zeta[(i)][(indeyr,indexr)][0]
                    # nodes['time'].append(self.mTime[t-i-1])
                else:
                    waterdepth = self.h[(indeyr,indexr)]+self.zeta[(i)][(indeyr,indexr)][0]
                    #nodes['time'].append(self.mTime[i+1])

                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))
            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

            if track_way=='backward': # backwards case
                u_t = -1*self.u[i,layer][(indeyu,indexu)][0][0] 
                v_t = -1*self.v[i,layer][(indeyv,indexv)][0][0]
            else:
                u_t = self.u[i,layer][(indeyu,indexu)][0][0] 
                v_t = self.v[i,layer][(indeyv,indexv)][0][0]
            # 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)
            
            
        return nodes
        
class get_fvcom(track):
    def __init__(self, mod):
        self.modelname = mod
            
    def nearest_point(self, lon, lat, lons, lats, length):  #0.3/5==0.06
        print '''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 Exception()
        #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
        mindex = np.argmin(distancelist)
        # location the point
        lonp = insidep[mindex][0]; latp = insidep[mindex][1]
        print lonp, latp
        return lonp,latp
        
    def uvreflection(self,f1,f2,b1,b2,v):
        '''
        Return reflection point.
        f2 is a forcast point from f1. 
        b1,b2 are boundary point.
        v=(u,v)
        Format (lon,lat)
        '''
        def rotate_vector(v, angle):
            """Rotate a vector `v` by the given angle(radian)."""
            x, y = v
        
            cos_theta = math.cos(angle)
            sin_theta = math.sin(angle)
        
            nx = x*cos_theta - y*sin_theta
            ny = x*sin_theta + y*cos_theta
        
            return [nx, ny]
            
        fl = Line(f1,f2)
        bl = Line(b1,b2)
        anchor = intersection(fl,bl)[0]
        
        anchor = str(anchor)[6:-1]
        anchor = anchor.split(','); #print 'anchor:', anchor
        interpx = Fraction(anchor[0]); interpy = Fraction(anchor[1])    
        anchor = [float(interpx),float(interpy)]
        #angle : angle in radians'''
        ag1 = Line(anchor,b1); ag2 = Line(anchor,b2); ab = Line(anchor,f2)
        ang = ab.angle_between(ag1)
        if  ang > pi/2:
            reag = ab.angle_between(ag2)
            bline = ag2
        else :
            reag = ang
            bline = ag1
        #print reag
        v1 = rotate_vector(v, reag*2)
        dx = 60*60*v1[0]; dy = 60*60*v1[1]
        temlon = anchor[0] + (dx/(111111*np.cos(anchor[1]*np.pi/180)))
        temlat = anchor[1] + dy/111111 #'''
        fp = (temlon,temlat) #reflection point
        
        angle2 = 2*pi-reag*2 ; #print 'angle2:'
        v2 = rotate_vector(v, angle2)
        dx = 60*60*v2[0]; dy = 60*60*v2[1]
        temlon2 = anchor[0] + (dx/(111111*np.cos(anchor[1]*np.pi/180)))
        temlat2 = anchor[1] + dy/111111 #'''
        fp2 = (temlon2,temlat2) #reflection point
        
        tg = Line(anchor,fp)
        tg2 = Line(anchor,fp2)
        #print tg.angle_between(bline)-reag
        if abs(tg.angle_between(bline)-reag)<abs(tg2.angle_between(bline)-reag) :       
            fp = (round(fp[0],6),round(fp[1],6))
        else:
            fp = (round(fp2[0],6),round(fp2[1],6))
            
        return fp #,anchor
            
    def get_url(self, starttime, endtime):
        print '''
        get different url according to starttime and endtime.
        urls are monthly.
        '''
        self.hours = int(round((endtime-starttime).total_seconds()/60/60))
        #print self.hours
                
        if self.modelname == "GOM3":
            timeurl = '''http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc?Times[0:1:144]'''
            url = '''http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc?
            lon[0:1:51215],lat[0:1:51215],lonc[0:1:95721],latc[0:1:95721],siglay[0:1:39][0:1:51215],h[0:1:51215],nbe[0:1:2][0:1:95721],
            u[{0}:1:{1}][0:1:39][0:1:95721],v[{0}:1:{1}][0:1:39][0:1:95721],zeta[{0}:1:{1}][0:1:51215]'''
            '''urll = http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc?
            u[{0}:1:{1}][0:1:39][0:1:95721],v[{0}:1:{1}][0:1:39][0:1:95721],zeta[{0}:1:{1}][0:1:51215]'''
            try:
                mTime = netCDF4.Dataset(timeurl).variables['Times'][:]              
            except:
                print '"GOM3" database is unavailable!'
                raise Exception()
            Times = []
            for i in mTime:
                strt = '201'+i[3]+'-'+i[5]+i[6]+'-'+i[8]+i[9]+' '+i[11]+i[12]+':'+i[14]+i[15]
                Times.append(datetime.strptime(strt,'%Y-%m-%d %H:%M'))
            # print Times
            fmodtime = Times[0]; emodtime = Times[-1]         
            if starttime<fmodtime or starttime>emodtime or endtime<fmodtime or endtime>emodtime:
                print 'Time: Error! Model(GOM3) only works between %s with %s(UTC).'%(fmodtime,emodtime)
                raise Exception()
            npTimes = np.array(Times)
            tm1 = npTimes-starttime; #tm2 = mtime-t2
            index1 = np.argmin(abs(tm1))
            index2 = index1 + self.hours#'''
            print '================'
            print index1, index2
            url = url.format(index1, index2)
            print url
            self.mTime = Times[index1:index2+1]
            
            self.url = url
            
        elif self.modelname == "massbay":
            timeurl = '''http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc?Times[0:1:144]'''
            url = """http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc?
            lon[0:1:98431],lat[0:1:98431],lonc[0:1:165094],latc[0:1:165094],siglay[0:1:9][0:1:98431],h[0:1:98431],
            nbe[0:1:2][0:1:165094],u[{0}:1:{1}][0:1:9][0:1:165094],v[{0}:1:{1}][0:1:9][0:1:165094],zeta[{0}:1:{1}][0:1:98431]"""
            
            try:
                mTime = netCDF4.Dataset(timeurl).variables['Times'][:]              
            except:
                print '"massbay" database is unavailable!'
                raise Exception()
            Times = []
            for i in mTime:
                strt = '201'+i[3]+'-'+i[5]+i[6]+'-'+i[8]+i[9]+' '+i[11]+i[12]+':'+i[14]+i[15]
                Times.append(datetime.strptime(strt,'%Y-%m-%d %H:%M'))
            fmodtime = Times[0]; emodtime = Times[-1]         
            if starttime<fmodtime or starttime>emodtime or endtime<fmodtime or endtime>emodtime:
                print 'Time: Error! Model(massbay) only works between %s with %s(UTC).'%(fmodtime,emodtime)
                raise Exception()
            npTimes = np.array(Times)
            tm1 = npTimes-starttime; #tm2 = mtime-t2
            index1 = np.argmin(abs(tm1))
            index2 = index1 + self.hours#'''
            url = url.format(index1, index2)
            self.mTime = Times[index1:index2+1]
            
            self.url = url

        elif self.modelname == "30yr": #start at 1977/12/31 23:00, end at 2014/1/1 0:0, time units:hours
            timeurl = """http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3?time[0:1:316008]"""
            url = '''http://www.smast.umassd.edu:8080/thredds/dodsC/fvcom/hindcasts/30yr_gom3?h[0:1:48450],
            lat[0:1:48450],latc[0:1:90414],lon[0:1:48450],lonc[0:1:90414],nbe[0:1:2][0:1:90414],siglay[0:1:44][0:1:48450],
            u[{0}:1:{1}][0:1:44][0:1:90414],v[{0}:1:{1}][0:1:44][0:1:90414],zeta[{0}:1:{1}][0:1:48450]'''
            
            try:
                mtime = netCDF4.Dataset(timeurl).variables['time'][:]
            except:
                print '"30yr" database is unavailable!'
                raise Exception
            # get model's time horizon(UTC).
            '''fmodtime = datetime(1858,11,17) + timedelta(float(mtime[0]))
            emodtime = datetime(1858,11,17) + timedelta(float(mtime[-1]))
            mstt = fmodtime.strftime('%m/%d/%Y %H:%M')
            mett = emodtime.strftime('%m/%d/%Y %H:%M') #'''
            # get number of days from 11/17/1858
            t1 = (starttime - datetime(1858,11,17)).total_seconds()/86400 
            t2 = (endtime - datetime(1858,11,17)).total_seconds()/86400
            if not mtime[0]<t1<mtime[-1] or not mtime[0]<t2<mtime[-1]:
                #print 'Time: Error! Model(massbay) only works between %s with %s(UTC).'%(mstt,mett)
                print 'Time: Error! Model(massbay) only works between 1978-1-1 with 2014-1-1(UTC).'
                raise Exception()
            
            tm1 = mtime-t1; #tm2 = mtime-t2
            index1 = np.argmin(abs(tm1)); #index2 = np.argmin(abs(tm2)); print index1,index2
            index2 = index1 + self.hours
            url = url.format(index1, index2)
            Times = []
            for i in range(self.hours+1):
                Times.append(starttime+timedelta(hours=i))
            self.mTime = Times
            self.url = url
        #print url
        return url

    def get_data(self,url):
        print '''
        "get_data" not only returns boundary points but defines global attributes to the object
        '''
        self.data = get_nc_data(url,'lat','lon','latc','lonc','siglay','h','nbe','u','v','zeta')#,'nv'
        self.lonc, self.latc = self.data['lonc'][:], self.data['latc'][:]  #quantity:165095
        self.lons, self.lats = self.data['lon'][:], self.data['lat'][:]
        self.h = self.data['h'][:]; self.siglay = self.data['siglay'][:]; #nv = self.data['nv'][:]
        self.u = self.data['u']; self.v = self.data['v']; self.zeta = self.data['zeta']

        nbe1=self.data['nbe'][0];nbe2=self.data['nbe'][1];
        nbe3=self.data['nbe'][2]
        pointt = np.vstack((nbe1,nbe2,nbe3)).T; self.pointt = pointt
        wl=[]
        for i in pointt:
            if 0 in i: 
                wl.append(1)
            else:
                wl.append(0)
        self.wl = wl
        tf = np.array(wl)
        inde = np.where(tf==True)
        #b_index = inde[0]
        lonb = self.lonc[inde]; latb = self.latc[inde]        
        self.b_points = np.vstack((lonb,latb)).T#'''
        #self.b_points = b_points
        return self.b_points #,nv lons,lats,lonc,latc,,h,siglay
        
    def shrink_data(self,lon,lat,lons,lats,rad):
        print '''shrink data'''
        print lon, lat
        print len(lons), len(lats)
        lont = []; latt = []
        p = Path.circle((lon,lat),radius=rad)
        pints = np.vstack((lons,lats)).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'
            sys.exit()
        print len(lonl), len(latl)
        return lonl,latl
    
    def eline_path(self,lon,lat):
        print '''
        When drifter close to boundary(less than 0.1),find one nearest point to drifter from boundary points, 
        then find two nearest boundary points to previous boundary point, create a boundary path using that 
        three boundary points.
        '''
        def boundary_location(locindex,pointt,wl):
            '''
            Return the index of boundary points nearest to 'locindex'.
            '''
            loca = []
            dx = pointt[locindex]; #print 'func',dx 
            for i in dx: # i is a number.
                #print i  
                if i ==0 :
                    continue
                dx1 = pointt[i-1]; #print dx1
                if 0 in dx1:
                    loca.append(i-1)
                else:
                    for j in dx1:
                        if j != locindex+1:
                            if wl[j-1] == 1:
                                loca.append(j-1)
            return loca
        
        p = Path.circle((lon,lat),radius=0.02) #0.06
        dis = []; bps = []; pa = []
        tlons = []; tlats = []; loca = []
        for i in self.b_points:
            if p.contains_point(i):
                bps.append((i[0],i[1]))
                d = math.sqrt((lon-i[0])**2+(lat-i[1])**2)
                dis.append(d)
        bps = np.array(bps)
        if not dis:
            return None
        else:
            print "Close to boundary."
            dnp = np.array(dis)
            dmin = np.argmin(dnp)
            lonp = bps[dmin][0]; latp = bps[dmin][1]
            index1 = np.where(self.lonc==lonp)
            index2 = np.where(self.latc==latp)
            elementindex = np.intersect1d(index1,index2)[0] # location 753'''
            #print index1,index2,elementindex  
            loc1 = boundary_location(elementindex,self.pointt,self.wl) ; #print 'loc1',loc1
            loca.extend(loc1)
            loca.insert(1,elementindex)               
            for i in range(len(loc1)):
                loc2 = boundary_location(loc1[i],self.pointt,self.wl); #print 'loc2',loc2
                if len(loc2)==1:
                    continue
                for j in loc2:
                    if j != elementindex:
                        if i ==0:
                            loca.insert(0,j)
                        else:
                            loca.append(j)
            
            for i in loca:
                tlons.append(self.lonc[i]); tlats.append(self.latc[i])
                       
            for i in xrange(len(tlons)):
                pa.append((tlons[i],tlats[i]))
            #path = Path(pa)#,codes
            return pa
            
    def streamlinedata(self,nodes,depth,track_way):
        
        lonpps,latpps,US,VS,speeds = [],[],[],[],[]
        #uvz = netCDF4.Dataset(self.url)
        #u = uvz.variables['u']; v = uvz.variables['v']; zeta = uvz.variables['zeta']
        lons = nodes['lons']
        lats = nodes['lats']
        size = max((max(lons)-min(lons)),(max(lats)-min(lats)))/1
        ilon,mlon = min(lons)-size,max(lons)+size
        ilat,mlat = min(lats)-size,max(lats)+size
        msize = max((mlon-ilon),(mlat-mlat))*0.71  # shrink model radius      
        
        lonmax = max(lons); lonmin = min(lons)
        latmax = max(lats); latmin = min(lats)
        lon1 = (lonmax+lonmin)/2; lat1 = (latmax+latmin)/2
        
        
        latps, lonps = np.mgrid[ilat:mlat:30j, ilon:mlon:30j]#x = np.linspace(-3,3,100)       
        latss = latps.flatten();lonss = lonps.flatten()
        
        lonll,latll = self.shrink_data(lon1,lat1,self.lonc,self.latc,msize)
        lonk,latk = self.shrink_data(lon1,lat1,self.lons,self.lats,msize)
        
        nodeindexs = []; elementindexs = []
        for j in xrange(len(latss)):
            try:
                if self.modelname == "GOM3" or self.modelname == "30yr":
                    lonpp,latpp = self.nearest_point(lonss[j],latss[j],lonll,latll,0.2)
                    lonn,latn = self.nearest_point(lonss[j],latss[j],lonk,latk,0.3)
                if self.modelname == "massbay":
                    lonpp,latpp = self.nearest_point(lonss[j],latss[j],lonll,latll,0.03)
                    lonn,latn = self.nearest_point(lonss[j],latss[j],lonk,latk,0.05)
                
                index11 = np.where(self.lonc==lonpp)
                index22 = np.where(self.latc==latpp)
                elementindex = np.intersect1d(index11,index22); #print index00
                index3 = np.where(self.lons==lonn)
                index4 = np.where(self.lats==latn)
                nodeindex = np.intersect1d(index3,index4)
                nodeindexs.append(nodeindex); elementindexs.append(elementindex)
            except:
                nodeindexs.append(0); elementindexs.append(0)
        t = abs(self.hours)+1
        for i in xrange(t):
            us,vs = [],[]
            print 'Streamline the %dth of %d points.'%(i+1,t)
            for j in xrange(len(latss)):
                if nodeindexs[j]==0 or elementindex==0 :
                    mi = 0.000000000000000
                    us.append(mi); vs.append(mi)
                else :             
                    try:
                        ni = nodeindexs[j]; ei = elementindexs[j]
                        if track_way=='backward' : # backwards case
                            waterdepth = self.h[ni]+self.zeta[t-i-1,ni]
                        else:
                            waterdepth = self.h[ni]+self.zeta[i,ni]
                        
                        depth_total = self.siglay[:,ni]*waterdepth  
                        layer = np.argmin(abs(depth_total-depth))                    
                        
                        if track_way=='backward' : # backwards case
                            us.append(self.u[t-1-i,layer,ei][0]); vs.append(self.v[t-1-i,layer,ei][0])
                        else:
                            us.append(self.u[i,layer,ei][0]); vs.append(self.v[i,layer,ei][0])                        
                    except:
                        #print 'out of area'
                        mi = 0.000000000000000
                        us.append(mi); vs.append(mi)
            sh = np.shape(lonps); #print sh,len(us)
            U = np.array(us).reshape(sh); V = np.array(vs).reshape(sh)
            speed = np.sqrt(U*U + V*V);print 'max speed',max(speed.flatten())
            
            lonpps.append(lonps); latpps.append(latps);
            US.append(U); VS.append(V); speeds.append(speed)
        
        return lonpps,latpps,US,VS,speeds
        
    def uvt(self,u1,v1,u2,v2):
        t = 2
        a=0; b=0
        if u1==u2:
            a = u1
        else:
            ut = np.arange(u1,u2,float(u2-u1)/t)
            for i in ut:
                a += i
            a = a/t  
        
        if v1==v2:
            b = v1
        else:
            c = float(v2-v1)/t
            vt = np.arange(v1,v2,c)
            for i in vt:
                b += i
            b = b/t
               
        return a, b
        
    def get_track(self,lon,lat,depth,track_way,bcon): #,b_index,nvdepth, 
        print '''
        Get forecast points start at lon,lat
        '''
        print lon, lat, depth, track_way, bcon
        modpts = dict(lon=[lon], lat=[lat], layer=[], time=[]) #model forecast points
        #uvz = netCDF4.Dataset(self.url)
        #u = uvz.variables['u']; v = uvz.variables['v']; zeta = uvz.variables['zeta']
        #print 'len u',len(u)
        if lon>90:
            lon, lat = dm2dd(lon, lat)
        lonl,latl = self.shrink_data(lon,lat,self.lonc,self.latc,0.5)
        lonk,latk = self.shrink_data(lon,lat,self.lons,self.lats,0.5)
        try:
            if self.modelname == "GOM3" or self.modelname == "30yr":
                lonp,latp = self.nearest_point(lon, lat, lonl, latl,0.2)
                lonn,latn = self.nearest_point(lon,lat,lonk,latk,0.3)
            if self.modelname == "massbay":
                lonp,latp = self.nearest_point(lon, lat, lonl, latl,0.03)
                lonn,latn = self.nearest_point(lon,lat,lonk,latk,0.05)        
            index1 = np.where(self.lonc==lonp)
            index2 = np.where(self.latc==latp)
            print index1, index2
            elementindex = np.intersect1d(index1,index2)
            print elementindex
            index3 = np.where(self.lons==lonn)
            index4 = np.where(self.lats==latn)
            print index3, index4
            nodeindex = np.intersect1d(index3,index4); #print 'here index'
            print nodeindex
            ################## boundary 1 ####################
            pa = self.eline_path(lon,lat); #print 'here boundary_path'
            print 'pa: ', pa
            if track_way=='backward' : # backwards case
                if self.modelname == "30yr":
                    waterdepth = self.h[nodeindex]
                else:                
                    waterdepth = self.h[nodeindex]+self.zeta[-1,nodeindex]
                modpts['time'].append(self.mTime[-1])
            else:
                #print 'h,zeta',self.h[nodeindex],self.zeta[0,nodeindex]
                if self.modelname == "30yr":
                    waterdepth = self.h[nodeindex]
                else:
                    waterdepth = self.h[nodeindex]+self.zeta[0,nodeindex]
                modpts['time'].append(self.mTime[0])
            depth_total = self.siglay[:,nodeindex]*waterdepth; #print 'Here one' 
            layer = np.argmin(abs(depth_total+depth)); #print 'layer',layer
            modpts['layer'].append(layer); 
            if waterdepth<(abs(depth)): 
                print 'This point is too shallow.Less than %d meter.'%abs(depth)
                raise Exception()
        except:
            return modpts,0  
            
        t = abs(self.hours) 
        print 't: ',t
        #mapx = Basemap(projection='ortho',lat_0=lat,lon_0=lon,resolution='l')
        for i in xrange(t):            
            if i!=0 and i%24==0 :
                #print 'layer,lon,lat,i',layer,lon,lat,i
                lonl,latl = self.shrink_data(lon,lat,self.lonc,self.latc,0.5)
                lonk,latk = self.shrink_data(lon,lat,self.lons,self.lats,0.5)
            if track_way=='backward' : # backwards case
                u_t1 = self.u[t-i,layer,elementindex][0]*(-1); v_t1 = self.v[t-i,layer,elementindex][0]*(-1)
                u_t2 = self.u[t-i-1,layer,elementindex][0]*(-1); v_t2 = self.v[t-i-1,layer,elementindex][0]*(-1)
            else:
                u_t1 = self.u[i,layer,elementindex][0]; v_t1 = self.v[i,layer,elementindex][0]
                u_t2 = self.u[i+1,layer,elementindex][0]; v_t2 = self.v[i+1,layer,elementindex][0]
            u_t,v_t = self.uvt(u_t1,v_t1,u_t2,v_t2)
            #u_t = (u_t1+u_t2)/2; v_t = (v_t1+v_t2)/2
            '''if u_t==0 and v_t==0: #There is no water
                print 'Sorry, hits the land,u,v==0'
                return modpts,1 #'''
            #print "u[i,layer,elementindex]",u[i,layer,elementindex]
            dx = 60*60*u_t; dy = 60*60*v_t
                                    
            #x,y = mapx(lon,lat)
            #temlon,temlat = mapx(x+dx,y+dy,inverse=True)            
            temlon = lon + (dx/(111111*np.cos(lat*np.pi/180)))
            temlat = lat + dy/111111 #'''
            #########case for boundary 1 #############
            if pa:
                pat = Path(pa)
                teml = [(lon,lat),(temlon,temlat)]
                #print temlon,temlat
                tempa = Path(teml)
                if pat.intersects_path(tempa):
                    if bcon == 'stop':
                        modpts['lon'].append(temlon); modpts['lat'].append(temlat); modpts['layer'].append(0)
                        print 'Sorry, point hits land here. Path. Last point:',temlon,temlat
                        return modpts,1 #'''
                    if bcon == 'reflection':
                        f1 = (lon,lat); f2 = (temlon,temlat); 
                        v = (u_t1,v_t1)
                        #fl = Path([f1,f2])
                        for k in range(len(pa)-1):
                            kl = Path([pa[k],pa[k+1]])
                            if tempa.intersects_path(kl):
                                print 'reflection',f1,f2,pa[k],pa[k+1],v
                                (temlon,temlat) = self.uvreflection(f1,f2,pa[k],pa[k+1],v)
                                break
                        
                
            #########case for boundary 2 #############
            '''if pa :
                if not pa.contains_point([temlon,temlat]):
                    print 'Sorry, point hits land here.path'
                    return modpts,1 #'''
            #########################
            lon = temlon; lat = temlat
            #if i!=(t-1):                
            try:
                if self.modelname == "GOM3" or self.modelname == "30yr":
                    lonp,latp = self.nearest_point(lon, lat, lonl, latl,0.2)
                    lonn,latn = self.nearest_point(lon,lat,lonk,latk,0.3)
                if self.modelname == "massbay":
                    lonp,latp = self.nearest_point(lon, lat, lonl, latl,0.03)
                    lonn,latn = self.nearest_point(lon,lat,lonk,latk,0.05)
                index1 = np.where(self.lonc==lonp)
                index2 = np.where(self.latc==latp)
                elementindex = np.intersect1d(index1,index2); #print 'elementindex',elementindex
                index3 = np.where(self.lons==lonn)
                index4 = np.where(self.lats==latn)
                nodeindex = np.intersect1d(index3,index4)
                
                ################## boundary 1 ####################
        
                pa = self.eline_path(lon,lat)
                
                if track_way=='backward' : # backwards case
                    if self.modelname == "30yr":
                        waterdepth = self.h[nodeindex]
                    else:                
                        waterdepth = self.h[nodeindex]+self.zeta[t-i-1,nodeindex]
                    modpts['time'].append(self.mTime[t-i-1])
                else:
                    #print 'h,zeta',self.h[nodeindex],self.zeta[0,nodeindex]
                    if self.modelname == "30yr":
                        waterdepth = self.h[nodeindex]
                    else:
                        waterdepth = self.h[nodeindex]+self.zeta[i+1,nodeindex]
                    modpts['time'].append(self.mTime[i+1])
                
                depth_total = self.siglay[:,nodeindex]*waterdepth  
                layer = np.argmin(abs(depth_total+depth)) 
                modpts['lon'].append(lon); modpts['lat'].append(lat); modpts['layer'].append(layer); 
                print '%d,lat,lon,layer'%(i+1),temlat,temlon,layer
                if waterdepth<(abs(depth)): 
                    print 'This point hits the land here.Less than %d meter.'%abs(depth)
                    raise Exception()
            except:
                return modpts,1
        print 'modpts: ',modpts
        return modpts,2

class get_wind(track):
    '''
    Take a point,start time, end time, return lists of lons,lats,us,vs in the circle area with the center of given point
    '''
    
    def __init__(self):
        pass
    
    def nearest_point(self, lon, lat, lons, lats, length=0.06):  #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 Exception()
        #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
        mindex = np.argmin(distancelist)
        # location the point
        lonp = insidep[mindex][0]; latp = insidep[mindex][1]
        
        return lonp,latp
        
    def get_url(self, starttime, endtime):
        '''
        get url according to starttime and endtime.
        '''
        starttime = starttime
        self.hours = int((endtime-starttime).total_seconds()/3600) # get total hours
        # time_r = datetime(year=2006,month=1,day=9,hour=1,minute=0)
        
        url_oceantime = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_MET_FORECAST.nc?Times[0:1:72]'
        url = """http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_MET_FORECAST.nc?XLAT[0:1:140][0:1:182],
                 XLONG[0:1:140][0:1:182],U10[{0}:1:{1}][0:1:140][0:1:182],V10[{0}:1:{1}][0:1:140][0:1:182]"""
        try:
            mTime = netCDF4.Dataset(url_oceantime).variables['Times'][:]
        except:
            print 'Wind database is unavailable!'
            raise Exception()
        # get model works time horizon(UTC).
        Times = []
        for i in mTime:
            strt = '201'+i[3]+'-'+i[5]+i[6]+'-'+i[8]+i[9]+' '+i[11]+i[12]+':'+i[14]+i[15]
            Times.append(datetime.strptime(strt,'%Y-%m-%d %H:%M'))
        fmodtime = Times[0]; emodtime = Times[-1]         
        if starttime<fmodtime or starttime>emodtime or endtime<fmodtime or endtime>emodtime:
            print 'Time: Error! Model(massbay) only works between %s with %s(UTC).'%(fmodtime,emodtime)
            raise Exception()
        npTimes = np.array(Times)
        tm1 = npTimes-starttime; #tm2 = mtime-t2
        index1 = np.argmin(abs(tm1))
        index2 = index1 + self.hours#'''
        url = url.format(index1, index2)
        self.mTime = Times[index1:index2+1]
        
        self.url = url
        
        return url
    
    def shrink_data(self,lon,lat,lons,lats,radius):
        lont = []; latt = []
        p = Path.circle((lon,lat),radius)
        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'
            sys.exit()
        return lonl,latl
        
    def get_data(self, url):
        '''
        return the data needed.
        url is from get_roms.get_url(starttime, endtime)
        '''
        
        data = get_nc_data(url,'XLAT','XLONG','Times','U10','V10')
        self.xlat = data['XLAT'][:]; self.xlong = data['XLONG'][:] 
        #self.times = data['Times'][:]
        self.u10,self.v10 = data['U10'], data['V10']
        
        return self.mTime
        
    def get_uv(self,points):#, depth
        '''
        get the nodes of specific time period
        lon, lat: start point
        depth: 0~35, the 0th is the bottom.
        '''
        
        lons = points['lons']
        lats = points['lats']
        mlon,mlat = (max(lons)+min(lons))/2,(max(lats)+min(lats))/2
        radius = max((max(lons)-min(lons))/2,(max(lats)-min(lats))/2)
                
        
        ixlong,ixlat = self.shrink_data(mlon,mlat,self.xlong,self.xlat, 2*radius)
        
        #nodes = dict(lon=[lon], lat=[lat], time=[])
        uss,vss = [],[]
        t = len(self.mTime)
        for j in xrange(t):         
            us,vs = [],[]
            for i in range(len(ixlong)):
                indexr = np.where(self.xlong==ixlong[i]); #print indexr
                #indexlat = np.where(self.xlat==ixlat[i])
                #indexr = np.intersect1d(indexlon,indexlat)
                #print indexlon,indexlat,indexr
                if len(indexr[0]) > 1:
                    print 'Repeat error.',len(indexr),indexr
                    raise Exception()#'''
                u,v = self.u10[j][indexr][0],self.v10[j][indexr][0]
                us.append(u); vs.append(v)
            uss.append(us); vss.append(vs)
        
        return ixlong,ixlat,uss,vss
class get_drifter(track):

    def __init__(self, drifter_id, filename=None):
        self.drifter_id = drifter_id
        self.filename = filename
    def get_track(self, starttime=None, days=None):
        '''
        return drifter nodes
        if starttime is given, return nodes started from starttime
        if both starttime and days are given, return nodes of the specific time period
        '''
        if self.filename:
            temp=getrawdrift(self.drifter_id,self.filename)
        else:
            temp=getdrift(self.drifter_id)
        nodes = {}
        nodes['lon'] = np.array(temp[1])
        nodes['lat'] = np.array(temp[0])
        nodes['time'] = np.array(temp[2])
        #starttime = np.array(temp[2][0])
        if not starttime:
            starttime = np.array(temp[2][0])
        if days:
            endtime = starttime + timedelta(days=days)
            i = self.__cmptime(starttime, nodes['time'])
            j = self.__cmptime(endtime, nodes['time'])
            
            if i == j :
                print 'No data in given duration.'
                raise Exception()
                #return nodes
            nodes['lon'] = nodes['lon'][i:j+1]
            nodes['lat'] = nodes['lat'][i:j+1]
            nodes['time'] = nodes['time'][i:j+1]
        else:
            i = self.__cmptime(starttime, nodes['time'])
            nodes['lon'] = nodes['lon'][i:-1]
            nodes['lat'] = nodes['lat'][i:-1]
            nodes['time'] = nodes['time'][i:-1]
        return nodes
        
    def __cmptime(self, time, times):
        '''
        return indies of specific or nearest time in times.
        '''
        tdelta = []
        #print len(times)
        for t in times:
            tdelta.append(abs((time-t).total_seconds()))
            
        index = tdelta.index(min(tdelta))
        
        return index
class get_roms_rk4(get_roms):
    '''
    model roms using Runge Kutta
    '''
    def get_track(self, lon, lat, depth, url):
        '''
        get the nodes of specific time period
        lon, lat: start point
        url: get from get_url(starttime, endtime)
        depth: 0~35, the 36th is the bottom.
        '''
        self.startpoint = lon, lat
        
        if type(url) is str:
            nodes = self.__get_track(lon, lat, depth, url)
            
        else: # case where there are two urls, one for start and one for stop time
            nodes = dict(lon=[self.startpoint[0]],lat=[self.startpoint[1]])
            
            for i in url:
                temp = self.__get_track(nodes['lon'][-1], nodes['lat'][-1], depth, i)
                nodes['lon'].extend(temp['lon'][1:])
                nodes['lat'].extend(temp['lat'][1:])
                
        return nodes # dictionary of lat and lon
        
    def __get_track(self, lon, lat, depth, url):
        '''
        ???? ????
        '''
        data = self.get_data(url)
        nodes = dict(lon=lon, lat=lat)
        mask = data['mask_rho'][:]
        lon_rho = data['lon_rho'][:]
        lat_rho = data['lat_rho'][:]
        index, nearestdistance = self.nearest_point_index(lon,lat,lons,lats)
        depth_layers = data['h'][index[0][0]][index[0][1]]*data['s_rho']
        layer = np.argmin(abs(depth_layers+depth))
        u = data['u'][:,layer]
        v = data['v'][:,layer]
        
        for i in range(0, len(data['u'][:])):
            u_t = u[i, :-2, :]
            v_t = v[i, :, :-2]
            lon, lat, u_p, v_p = self.RungeKutta4_lonlat(lon,lat,lons,lats,u_t,v_t)
            
            if not u_p:
                print 'point hit the land'
                break
            nodes['lon'] = np.append(nodes['lon'],lon)
            nodes['lat'] = np.append(nodes['lat'],lat)
            
        return nodes
        
    def polygonal_barycentric_coordinates(self,xp,yp,xv,yv):
        '''
        ??? how is this one solved???
        '''
        N=len(xv)   
        j=np.arange(N)
        ja=(j+1)%N
        jb=(j-1)%N
        Ajab=np.cross(np.array([xv[ja]-xv[j],yv[ja]-yv[j]]).T,
                      np.array([xv[jb]-xv[j],yv[jb]-yv[j]]).T)
        Aj=np.cross(np.array([xv[j]-xp,yv[j]-yp]).T,
                    np.array([xv[ja]-xp,yv[ja]-yp]).T)
        Aj=abs(Aj)
        Ajab=abs(Ajab)
        Aj=Aj/max(abs(Aj))
        Ajab=Ajab/max(abs(Ajab))    
        w=xv*0.
        j2=np.arange(N-2)
        
        for j in range(N):
            
            w[j]=Ajab[j]*Aj[(j2+j+1)%N].prod()
          
        w=w/w.sum()
        
        return w
        
    def VelInterp_lonlat(self,lonp,latp,lons,lats,u,v):
        index, distance = self.nearest_point_index(lonp,latp,lons,lats)
        lonv,latv = lons[index[0],index[1]], lats[index[0],index[1]]
        w = self.polygonal_barycentric_coordinates(lonp,latp,lonv,latv)
        uf = (u[index[0],index[1]]/np.cos(lats[index[0],index[1]]*np.pi/180)*w).sum()
        vf = (v[index[0],index[1]]*w).sum()
        
        return uf, vf
        
    def RungeKutta4_lonlat(self,lon,lat,lons,lats,u,v):
        '''
        ?????????????
        '''
        tau = 60*60/111111.
        lon1=lon*1.;          lat1=lat*1.;        urc1,v1=self.VelInterp_lonlat(lon1,lat1,lons,lats,u,v);  
        lon2=lon+0.5*tau*urc1;lat2=lat+0.5*tau*v1;urc2,v2=self.VelInterp_lonlat(lon2,lat2,lons,lats,u,v);
        lon3=lon+0.5*tau*urc2;lat3=lat+0.5*tau*v2;urc3,v3=self.VelInterp_lonlat(lon3,lat3,lons,lats,u,v);
        lon4=lon+    tau*urc3;lat4=lat+    tau*v3;urc4,v4=self.VelInterp_lonlat(lon4,lat4,lons,lats,u,v);
        lon=lon+tau/6.*(urc1+2.*urc2+2.*urc3+urc4);
        lat=lat+tau/6.*(v1+2.*v2+2.*v3+v4); 
        uinterplation=  (urc1+2.*urc2+2.*urc3+urc4)/6    
        vinterplation= (v1+2.*v2+2.*v3+v4)/6
       
        return lon,lat,uinterplation,vinterplation

def min_data(*args):
    '''
    return the minimum of several lists
    '''
    data = []
    for i in range(len(args)):    
        data.append(min(args[i]))
    return min(data)
    
def max_data(*args):
    '''
    return the maximum of several lists
    '''
    data = []   
    for i in range(len(args)):
        data.append(max(args[i]))
    return max(data)
    
def angle_conversion(a):
    '''
    converts the angle into radians
    '''
    a = np.array(a)
    
    return a/180*np.pi
    
def dist(lon1, lat1, lon2, lat2):
    '''
    calculate the distance of points
    '''
    R = 6371.004
    lon1, lat1 = angle_conversion(lon1), angle_conversion(lat1)
    lon2, lat2 = angle_conversion(lon2), angle_conversion(lat2)
    l = R*np.arccos(np.cos(lat1)*np.cos(lat2)*np.cos(lon1-lon2)+
                    np.sin(lat1)*np.sin(lat2))
                    
    return l

    
def draw_basemap(ax, points, interval_lon=0.3, interval_lat=0.3):
    '''
    draw the basemap?
    '''
    
    lons = points['lons']
    lats = points['lats']
    size = max((max(lons)-min(lons)),(max(lats)-min(lats)))/2
    map_lon = [min(lons)-size,max(lons)+size]
    map_lat = [min(lats)-size,max(lats)+size]
    
    #ax = fig.sca(ax)
    dmap = Basemap(projection='cyl',
                   llcrnrlat=map_lat[0], llcrnrlon=map_lon[0],
                   urcrnrlat=map_lat[1], urcrnrlon=map_lon[1],
                   resolution='h',ax=ax)# resolution: c,l,i,h,f.
    dmap.drawparallels(np.arange(int(map_lat[0])-1,
                                 int(map_lat[1])+1,interval_lat),
                       labels=[1,0,0,0])
    dmap.drawmeridians(np.arange(int(map_lon[0])-1,
                                 int(map_lon[1])+1,interval_lon),
                       labels=[0,0,0,1])
    dmap.drawcoastlines()
    dmap.fillcontinents(color='grey')
    dmap.drawmapboundary()
    #dmap.etopo()

def uniquecolors(N):
    """
    Generate unique RGB colors
    input: number of RGB tuples to generate
    output: list of RGB tuples
    """
    HSV_tuples = [(x*1.0/N, 0.5, 0.5) for x in range(N)]
    RGB_tuples = map(lambda x: colorsys.hsv_to_rgb(*x), HSV_tuples)
    colors =  [colorsys.hsv_to_rgb(x*1.0/N, 0.5, 0.5) for x in range(N)]
    return colors

def basemap_region(region):
    path="" # Y:/bathy/"#give the path if these data files are store elsewhere
    #if give the region, choose the filename
    if region=='sne':
        filename='sne_coast.dat'
    if region=='cc':
        filename='capecod_outline.dat'
    if region=='bh':
        filename='bostonharbor_coast.dat'
    if region=='cb':
        filename='cascobay_coast.dat'
    if region=='pb':
        filename='penbay_coast.dat'
    if region=='ma': # mid-atlantic
        filename='necscoast_noaa.dat'
    if region=='ne': # northeast
        filename='necoast_noaa.dat'   
    if region=='wv': # world vec
        filename='necscoast_worldvec.dat'        
    
    #open the data
    f=open(path+filename)

    lon,lat=[],[]
    for line in f:#read the lat, lon
        lon.append(line.split()[0])
        lat.append(line.split()[1])
    nan_location=[]
    # plot the lat,lon between the "nan"
    for i in range(len(lon)):#find "nan" location
        if lon[i]=="nan":
            nan_location.append(i)

    for m in range(1,len(nan_location)):#plot the lat,lon between nan
        lon_plot,lat_plot=[],[]
        for k in range(nan_location[m-1],nan_location[m]):
            lat_plot.append(lat[k])
            lon_plot.append(lon[k])
        plt.plot(lon_plot,lat_plot,'r') 

def clickmap(n):
   # this allows users to click on a rough map and define lat/lon points
   # where "n" is the number of points
   fig=plt.figure()
   basemap_region('cc')
   pt=fig.ginput(n)
   plt.close('all')
   lon=list(zip(*pt)[0])
   lat=list(zip(*pt)[1])
   return lon,lat

def points_between(st_point,en_point,x):
    """ 
    For 2 positions, interpolate X number of points between them
    where "lat" and "lon" are two element list
    "x" is the number of points wanted between them
    returns lat0,lono
    """
    
    lato=[]
    lono=[]
    if not st_point: 
        lato.append(en_point[0]); lono.append(en_point[1])
        return lato,lono
    if not en_point: 
        lato.append(st_point[0]); lono.append(st_point[1])
        return lato,lono
    lati=(en_point[0]-st_point[0])/float(x+1)
    loni=(en_point[1]-st_point[1])/float(x+1)
    for j in range(x+2):
        lato.append(st_point[0]+lati*j)
        lono.append(st_point[1]+loni*j)
    
    return lato,lono
    
def points_square(point, hside_length):
    '''point = (lat,lon); length: units is decimal degrees.
       return a squre points(lats,lons) on center point'''
    (lat,lon) = point; length =float(hside_length)
    lats=[lat]; lons=[lon]
    bbox = [lon-length, lon+length, lat-length, lat+length]
    bbox = np.array(bbox)
    points = np.array([bbox[[0,1,1,0]],bbox[[2,2,3,3]]])
    lats.extend(points[1]); lons.extend(points[0])
    
    return lats,lons

def extend_square(point, radius,num):
    '''point = (lat,lon); length: units is decimal degrees.
       return a squre points(lats,lons) on center point'''
    (lat,lon) = point; 
    lats=[lat]; lons=[lon]
    length =float(radius)/(num+1)
    for i in xrange(num+1):
        lth = length*(i+1)
        bbox = [lon,lon-lth,lon+lth,lat,lat-lth,lat+lth]
        bbox = np.array(bbox)
        points = np.array([bbox[[1,2,2,1,1,2,0,0]],bbox[[4,4,5,5,3,3,4,5]]])
        lats.extend(points[1]); lons.extend(points[0])

    return lats,lons
    
def extend_units(point, unit, num):
    '''point = (lat,lon); length: units is decimal degrees.
       return a squre points(lats,lons) on center point'''
    (lat,lon) = point; 
    lats=[]; lons=[]
    #unit = unit
    leh = unit*num
    
    #ps = np.mgrid[lat-leh:lat+leh+1:2j, lon-leh:lon+leh+1:2j]; print ps
    #ps = ps*unit
    lts = np.arange(lat-leh,lat+leh+unit,unit); lns = np.arange(lon-leh,lon+leh+unit,unit)
    llp = np.meshgrid(lts,lns); #print llp
    lats.extend(llp[0].flatten()); lons.extend(llp[1].flatten())

    return lats,lons

def totdis(lons,lats):
    "return path length of list points" 
    ts = 0
    lp = len(lons)-1
    for i in xrange(lp):
        dx = lons[i+1]-lons[i]
        dy = lats[i+1]-lats[i]
        ds = math.sqrt(dx**2+dy**2)
        ts += ds
    return ts
