Preprocessing of the SLA data ::
       * process positions :
          * interpolate over gaps
          * find continents (extend the positions over continents to get the discontinuity)
          * find track edges
          * find gap lengths

       * clean SLA data::
          * Remove gaps greater than maximum allowed length over which interpolate is OK.
          * Remove time steps with not enough coverage
          * get sub-segments of valid data of a given length
          
    :parameter lon: longitude
    :parameter lat: longitude
    :parameter sla: data
    
    :keyword N_min: Length of subsegments (cf :func:`altimetry.tools.spectrum.get_segments`)
    :keyword per_min: Minimum percentage of valid data to allow.
    :keyword max_gap: Maximum gap length to interpolate over (interpolation is done 1st, THEN long gaps are eliminated)
    :keyword leave_gaps: Leave gaps (equivalent to setting max_gap to number of points in track).
    
    :keyword remove_edges: discard data at track edges.
    :keyword truncate_if_continents: Force truncating data if a continent is found within a segment of data.
    :keyword last: Get segments of data sticked to the last point in track (cf :func:`altimetry.tools.spectrum.get_segments`)
    :keyword first: Get segments of data sticked to the first point in track (cf :func:`altimetry.tools.spectrum.get_segments`)
    :keyword mid: Get segments of data sticked to the middle point in track (cf :func:`altimetry.tools.spectrum.get_segments`)

def preprocess(lat,lon,sla,
               N_min=None,
               per_min=15.0,
               max_gap=None,
               leave_gaps=False,
               remove_edges=True,
               interp_over_continents=False,
               truncate_if_continents=True,
               discard_continental_gaps = True, #do not consider gaps on continents to discard time steps using max_gap.
               flag_interp=False,
               return_lonlat=False,
               return_interpolated=False,
               last=True,mid=None,first=None, #for get_segments
               verbose=1):
    '''   
    Preprocessing of the SLA data ::
       * process positions :
          * interpolate over gaps
          * find continents (extend the positions over continents to get the discontinuity)
          * find track edges
          * find gap lengths

       * clean SLA data::
          * Remove gaps greater than maximum allowed length over which interpolate is OK.
          * Remove time steps with not enough coverage
          * get sub-segments of valid data of a given length
          
    :parameter lon: longitude
    :parameter lat: longitude
    :parameter sla: data
    
    :keyword N_min: Length of subsegments (cf :func:`altimetry.tools.spectrum.get_segments`)
    :keyword per_min: Minimum percentage of valid data to allow.
    :keyword max_gap: Maximum gap length to interpolate over (interpolation is done 1st, THEN long gaps are eliminated)
    :keyword leave_gaps: Leave gaps (equivalent to setting max_gap to number of points in track).
    
    :keyword remove_edges: discard data at track edges.
    :keyword truncate_if_continents: Force truncating data if a continent is found within a segment of data.
    :keyword last: Get segments of data sticked to the last point in track (cf :func:`altimetry.tools.spectrum.get_segments`)
    :keyword first: Get segments of data sticked to the first point in track (cf :func:`altimetry.tools.spectrum.get_segments`)
    :keyword mid: Get segments of data sticked to the middle point in track (cf :func:`altimetry.tools.spectrum.get_segments`)
    
    '''
    sh=sla.shape
    nt=sh[0]
    nx=sh[1]
    dumsla=sla.copy()
        
    #Remove profiles with less than 3 points
    ok=np.where(sla.mask.sum(axis=1) < (nx -3))[0]
    if (nt != len(ok)) : message(2, '%i time steps on %i removed: contain less than 3 valid data points' % (nt - len(ok),nt),verbose=verbose)
    dumsla=dumsla[ok,:]
    ntinit=nt
    nt=len(ok)
    
    # 1 - regrid tracks regularly
    #     get gap lengths
    #############################
    fg_dumlon=nx
    for i in np.arange(nt):
        
        #grid track regularly
        fg=~dumsla.mask[i,:]
        dst, dumlon, dumlat, dsla, lgaps, n, edges, inter = grid_track(lat[fg], lon[fg], dumsla[i,:][fg],remove_edges=False,backbone=[lon,lat],interp_over_continents=interp_over_continents)
        if isinstance(dumlon,np.ma.masked_array) : fg_dumlon_new = dumlon.mask.sum() 
        else : fg_dumlon_new = np.isfinite(dumlon).sum()
        if fg_dumlon_new < fg_dumlon :
            fg_dumlon = fg_dumlon_new
            lonout = dumlon
            latout = dumlat
        
        #extend matrix width if track has gone over any land (ie. any points not found in the backbone)
        if (len(dumlon) > len(lon)) & (i == 0) :            
            lendiff = len(dumlon) - len(lon)
            print '[WARNING] : Pass goes over a land mass, changing the track size from {0} to {1}'.format(nx,nx+lendiff)
            nx+=lendiff

        #~ if i==0:
            #~ dumint=inter.reshape((1,len(dsla)))
            #~ dumslaout=dsla.reshape((1,len(dsla)))
        #~ else:
            #~ dumint=np.ma.array(np.concatenate([dumint.data,inter.data.reshape((1,len(dsla)))],axis=0),
                               #~ mask=np.concatenate([dumint.mask,inter.mask.reshape((1,len(dsla)))],axis=0))
            #~ dumslaout=np.ma.array(np.concatenate([dumslaout.data,dsla.data.reshape((1,len(dsla)))],axis=0),
                                  #~ mask=np.concatenate([dumslaout.data,dsla.data.reshape((1,len(dsla)))],axis=0))
        dumint=inter.reshape((1,len(dsla))) if i == 0 else np.ma.concatenate([dumint,inter.reshape((1,len(dsla)))],axis=0)
        dumslaout=dsla.reshape((1,len(dsla))) if i == 0 else np.ma.concatenate([dumslaout,dsla.reshape((1,len(dsla)))],axis=0)
    
        if i == 0 :
            gaplen = [lgaps]
            gapedges= [edges]
        else :
            gaplen.append(lgaps)
            gapedges.append(edges)
        ngaps = n if i == 0 else np.append(ngaps,n) 
    dumsla=dumslaout.copy()
    
    #These points are the continental ones.
    dumint=dumint.astype(bool)
    continent=dumsla.mask & dumint
    flagged=dumint.astype(bool) & ~continent
    
    #check for continental gaps
    iscont=np.sum(continent,axis=0) == nt
    indcont=np.arange(nx)[iscont]
    
    #here we get the position of the gap from its egdes and intersect with the continental points
    # -> if any point is over continent, cont_gap is set to True
    cont_gap = [[len(set(indcont).intersection(range(gapedges[j][0][jj],gapedges[j][1][jj]))) > 0 for jj in xrange(ngaps[j])] for j in xrange(nt)]

    if discard_continental_gaps :
        gaplen=[ np.array([g[jj] for jj in xrange(len(g)) if not cont_gap[j][jj]]) for j,g in enumerate(gaplen)]
    
    # 2 - remove/subsample using stats from previous loop and keyword arguments 
    ###########################################################################
    if max_gap is not None:
        
        #Remove profiles with long gaps
        gapmax=np.array([np.max(g) if len(g) > 0 else 0 for g in gaplen])
        id1 = np.where(gapmax <= max_gap)[0] if not leave_gaps else ok
        
        if len(id1) == 0 : raise Exception('[ERROR] : All gaps in current track are greater than the maximum specified gap')
        if (len(id1) != nt) : message(2, '%i time steps on %i removed: gaps > %i point' %(nt - len(id1), ntinit, int(max_gap)),verbose=verbose)
        
        dumsla=dumsla[id1,:]
        
        #Remove profiles with not enough coverage :
        # -> based on points with interpolation flag from grid_track and not on continents
#        per=100 * dumsla.mask.sum(axis=0) / np.float(nt)
        per = 100. * flagged[id1,:].sum(axis=1)/np.float(nx)
        if N_min is None :
#             N_min = np.round((0.01* (100- per_min)) * nx)
            N_min = nx
#        if per_min is None : per_min = 100* (1 - N_min/np.float(nx))
        
        id2 = np.where( per <= per_min)[0]
        
        if len(id2) == 0 : raise Exception('[ERROR] : All time steps in current track have a percentage of invalid data > than the maximum allowed (%i)' % int(per_min))
        if (len(id2) != len(id1)) : message(2, '%i time steps on %i removed: exceed maximum allowed percentage of invalid data (%i)' %(len(id1) - len(id2), ntinit, int(per_min)),verbose=verbose)
        
        
        dumsla=dumsla[id2,:]
        
        
        
        #At this point track edges are removed
        dumsla,id3=get_segment(dumsla,N_min,remove_edges=remove_edges,truncate_if_continents=truncate_if_continents,last=last,mid=mid,first=first)
        
        if len(id3) == 0 : raise Exception('[ERROR] : Remaining time steps do not reach the minimum length of %i points' % int(N_min))
        if (len(id3) != len(id2)) : message(2, '%i time steps no reaching rhe minimum length of %i points have been removed)' %(len(id2) - len(id3), int(N_min)),verbose=verbose)
        
        res=(dumsla, ok[id1[id2[id3]]])
        
        
    else :
        res=(dumsla, ngaps, gaplen)
    
    nt=res[0].shape[0]
    if (nt != ntinit) : message(1, '%i time steps on %i removed by data pre-processing' %(ntinit - nt, ntinit),verbose=verbose)
    if return_lonlat : res+=(lonout, latout)
    if return_interpolated : res += (dumint,)
    return res