# GRID_TRACK
    # @summary: This function allow detecting gaps in a set of altimetry data and rebin this data regularly, with informations on gaps.
    # @param lat {type:numeric} : latitude
    # @param lon {type:numeric} : longitude
    # @param sla {type:numeric} : data
    # @return:
    #    outdst : resampled distance
    #    outlon : resampled longitude
    #    outlat : resampled latitude
    #    outsla : resampled data
    #    gaplen : length of the longest gap in data
    #    dx : average spatial sampling
    #    interpolated : True when data was interpolated (empty bin)
    #
    # @author: Renaud DUSSURGET (RD) - LER/PAC, Ifremer
    # @change: Created by RD, July 2012
    #    29/08/2012 : Major change -> number of output variables changes (added INTERPOLATED), and rebinning modified
    #    06/11/2012 : Included in alti_tools lib
    #    19/12/2012 : Added backbone option (reproject along the backbone grid)

def grid_track(lat,lon,sla,remove_edges=None,backbone=None,interp_over_continents=True):
    """
    # GRID_TRACK
    # @summary: This function allow detecting gaps in a set of altimetry data and rebin this data regularly, with informations on gaps.
    # @param lat {type:numeric} : latitude
    # @param lon {type:numeric} : longitude
    # @param sla {type:numeric} : data
    # @return:
    #    outdst : resampled distance
    #    outlon : resampled longitude
    #    outlat : resampled latitude
    #    outsla : resampled data
    #    gaplen : length of the longest gap in data
    #    dx : average spatial sampling
    #    interpolated : True when data was interpolated (empty bin)
    #
    # @author: Renaud DUSSURGET (RD) - LER/PAC, Ifremer
    # @change: Created by RD, July 2012
    #    29/08/2012 : Major change -> number of output variables changes (added INTERPOLATED), and rebinning modified
    #    06/11/2012 : Included in alti_tools lib
    #    19/12/2012 : Added backbone option (reproject along the backbone grid)
    """
    
    
    
    #Find gaps in data
    if backbone is not None :
        backlon=backbone[0]
        backlat=backbone[1]
        ascending=track_orient(lon,lat)
        dst=calcul_distance(backlat[0],backlon[0],lat,lon)
        if ascending : dst[lat < backlat[0]]*=-1
        if not ascending : dst[lat > backlat[0]]*=-1
        dstback=calcul_distance(backlat,backlon)
        dx = dstback[1:] - dstback[:-1]
        mn_dx = np.median(dx)
        bins = np.round(dstback.max() / mn_dx)+1
        range=(0/2.,mn_dx * bins) - mn_dx/2
#        dfback=list(set(dstback).difference(set(dst)))
        bhist,bbin_edges=np.histogram(dstback, bins=bins, range=range)
        continent=np.where(bhist==0)[0]
        if remove_edges is None : remove_edges=False
    else :
        dst=calcul_distance(lat,lon)
        #Find gaps in data
        dx = dst[1:] - dst[:-1]
        mn_dx = np.median(dx)
        bins = np.ceil(dst.max() / mn_dx) + 1
        range=(0/2.,mn_dx * bins) - mn_dx/2
        
        if remove_edges is None : remove_edges=True
    
    
    hist,bin_edges=np.histogram(dst, bins=bins, range=range) #We have binned the data along a regular grid of size (bins) in the range (range)
                                                             #Missing data is thus represented by no data in a given bin
    
    #Remove leading and trailing edges (and synchronise bin_edges)
    if remove_edges == True :
        while hist[0] == 0 : 
            hist=np.delete(hist,[0])
            bin_edges=np.delete(bin_edges,[0])
        while hist[-1] == 0 :
            hist=np.delete(hist,[len(hist)-1])
            bin_edges=np.delete(bin_edges,[len(bin_edges)-1])
    
    nH =len(hist)
    
    #Get filled bins indices
    ok = np.arange(len(hist)).compress(np.logical_and(hist,True or False))
    empty = np.arange(len(hist)).compress(~np.logical_and(hist,True or False)) 
    
    if isinstance(sla,np.ma.masked_array) : outsla = np.ma.masked_array(np.repeat(sla.fill_value,nH),mask=np.ones(nH,dtype=bool),dtype=sla.dtype)
    else : outsla = np.ma.masked_array(np.repeat(np.ma.default_fill_value(1.0),nH),mask=np.ones(nH,dtype=bool),dtype=np.float32)
    if isinstance(sla,np.ma.masked_array) : outlon = np.ma.masked_array(np.repeat(lon.fill_value,nH),mask=np.ones(nH,dtype=bool),dtype=lon.dtype)
    else : outlon = np.ma.masked_array(np.repeat(np.ma.default_fill_value(1.0),nH),mask=np.ones(nH,dtype=bool),dtype=np.float32)
    if isinstance(sla,np.ma.masked_array) : outlat = np.ma.masked_array(np.repeat(lat.fill_value,nH),mask=np.ones(nH,dtype=bool),dtype=lat.dtype)
    else : outlat = np.ma.masked_array(np.repeat(np.ma.default_fill_value(1.0),nH),mask=np.ones(nH,dtype=bool),dtype=np.float32)
    
    outdst = bin_edges [:-1]+ mn_dx/2 #distances is taken at bins centers
    outsla[ok] = sla
    outlon[ok] = lon
    outlat[ok] = lat
    
    #Remove land mass point if asked
    if not interp_over_continents :
        sempty=np.sort(np.array(list(set(empty).difference(set(continent)))))
    else : sempty=empty.copy()
    
    #Fill the gaps if there are some
    if len(empty) > 0 : 
        #Interpolate lon,lat @ empty positions
        outlon[empty] = interp1d(ok, outlon[ok], empty, kind='cubic', fill_value=lon.fill_value)
        outlat[empty] = interp1d(ok, outlat[ok], empty, kind='cubic', fill_value=lat.fill_value)
    if len(sempty) > 0 :outsla[sempty] = interp1d(ok, outsla[ok], empty, kind=0, fill_value=sla.fill_value) #0-th order Spline interpolation
    
    outlon.mask[outlon.data == outlon.fill_value] = outlon.fill_value
    outlat.mask[outlat.data == outlat.fill_value] = outlat.fill_value
    outsla.mask[outsla.data == outsla.fill_value] = outsla.fill_value
    
    
    #Get gap properties
    ind=np.arange(len(hist))
    dhist=(hist[1:] - hist[:-1])
    
    #Remove edge gaps
    if (dhist!=0).sum()> 0 :
        if (dhist[dhist!=0])[0] == 1 : dhist[(np.arange(nH)[dhist!=0])[0]]=0 #This do not count start and end of track as gaps
    if (dhist!=0).sum()> 0 :
        if (dhist[dhist!=0])[-1] == -1 : dhist[(np.arange(nH)[dhist!=0])[-1]]=0
    st=ind[dhist==-1]+1
    en=ind[dhist==1]
    gaplen=(en-st) + 1
    ngaps=len(st)
    gapedges=np.array([st,en])
    
    #Get empty bin flag
    interpolated=~hist.astype('bool')
    
    return outdst, outlon, outlat, outsla, gaplen, ngaps, gapedges, interpolated