INTERP2D2D : Interpolate a 2D matrix into another 2D matrix
    
    @param x: 1st dimension vector of size NX
    @param y: 2nd dimension vector of size NY
    @param Z: Array to interpolate (NXxNY)
    @param xout: 1st dimension vector of size NXout
    @param yout: 2nd dimension vector of size NYout
    
    @keyword split_factor:Nummber of times to split arrays. Nb of threads is equal to split_factor**2. 
    
    @return: Interpolated array (NXoutxNYout)
    
    @author: Renaud DUSSURGET, LER/PAC, Ifremer La Seyne

def interp2d2d(x,y,Z,xout,yout,split_factor=1,**kwargs):
    """    
    INTERP2D2D : Interpolate a 2D matrix into another 2D matrix
    
    @param x: 1st dimension vector of size NX
    @param y: 2nd dimension vector of size NY
    @param Z: Array to interpolate (NXxNY)
    @param xout: 1st dimension vector of size NXout
    @param yout: 2nd dimension vector of size NYout
    
    @keyword split_factor:Nummber of times to split arrays. Nb of threads is equal to split_factor**2. 
    
    @return: Interpolated array (NXoutxNYout)
    
    @author: Renaud DUSSURGET, LER/PAC, Ifremer La Seyne
    """
    
    #Queued thread
    ##############
    class ThreadClass(threading.Thread):
        
        #We override the __init__ method
        def __init__(self, input_q,indices_q,result_q):
            threading.Thread.__init__(self)
            self.input = input_q
            self.indices = indices_q
            self.result = result_q
        
        def task(self,NaN=True):
            
            #grabs host from queue
            inargs = self.input.get()
            ind = inargs[0] #thread index
            x = inargs[1]
            y = inargs[2]
            Z = inargs[3]
            xout = inargs[4]
            yout = inargs[5]
            
            Zout=_interp2d2d(x,y,Z,xout,yout,**kwargs)
            
        #    Zout = sc.interpolate.griddata(points, gz, xi, **kwargs)
            
            self.indices.put(ind)
            self.result.put(Zout)
    
            del Zout
            
#            if verbose > 0 : print "%s ended at time: %s" % (self.getName(), datetime.datetime.now())
        
        def run(self):
            
                #Starts the queue
                self.task()
                #signals to queue job is done
                self.input.task_done()
       

    #Setup input and output queues
    input_q = Queue.Queue()
    indices_q = Queue.Queue()
    result_q = Queue.Queue()
    
    #Map the data along X axis
    N_threads=split_factor**2
#    over=np.ceil(cut/dx) #Overlay between each time series processed in parallel
    
    #Get dimensions to split matrix
    nxin=x.size
    nyin=y.size
    nxout=xout.size
    nyout=yout.size
    
    gx = np.reshape(np.repeat(x,nyin),(nxin,nyin))
    gy = np.reshape(np.repeat(y,nxin),(nyin,nxin)).transpose((1,0))

    gxout = np.reshape(np.repeat(xout,nyout),(nxout,nyout))
    gyout = np.reshape(np.repeat(yout,nxout),(nyout,nxout)).transpose((1,0))
    
    #Map output coordinates
    ind=[]
    
    #Map output coordinates
    xsplit=[]
    ysplit=[]
    for i in np.arange(split_factor) : 
        for j in np.arange(split_factor) :
            xsplit.append(i*(nxout/float(split_factor)))
            xsplit.append((i+1)*(nxout/float(split_factor))-1)
            ysplit.append(j*(nyout/float(split_factor)))
            ysplit.append((j+1)*(nyout/float(split_factor))-1)
    
    #Round
    xsplit=np.round(xsplit).astype(int)
    ysplit=np.round(ysplit).astype(int)
    
    N_threads = len(xsplit)/2
    
    
    
#    th_xout=gxout[0:nxout/split_factor,0:2]
#    th_yout=

#    gz = Z.flatten()
     
#    points = zip(*(gx.flatten(),gy.flatten())) 
    
    #spawn a pool of threads, and pass them queue instance 
    for i in np.arange(N_threads):
        t = ThreadClass(input_q,indices_q,result_q)
        t.setDaemon(True)
        t.start()
    
    #Feed threads with data
    for i in range(N_threads):
        xoind=xsplit[[i*2,(i*2)+1]]
        yoind=ysplit[[i*2,(i*2)+1]]
        xiind=x[xout[xoind].astype(int)].astype(int)
        yiind=y[yout[yoind].astype(int)].astype(int)
        th_x=x[xiind[0]:xiind[1]+1]
        th_y=y[yiind[0]:yiind[1]+1]
        th_xo=xout[xoind[0]:xoind[1]+1]
        th_yo=yout[yoind[0]:yoind[1]+1]
        th_z=Z[xiind[0]:xiind[1]+1,yiind[0]:yiind[1]+1]
#        input_q.put((i,np.copy(points),gz.copy(),xout[iind[i]].copy(),yout[iind[i]].copy()))
        input_q.put((i,th_x,th_y,th_z,th_xo,th_yo))

   
#    outvar=Z.copy()
#    outvar.data[:]=outvar.fill_value
    
#    for i in range(N_threads):
#        print "%s launched time: %s" % (i, datetime.datetime.now())
#        input_q.put((i,np.copy(points),gz.copy(),xout[iind[i]].copy(),yout[iind[i]].copy()))

    

    #wait on the queue until everything has been processed     
    input_q.join()
    
    #Sort threads
    tnb=[]
    for i in range(N_threads) : tnb.append(indices_q.get(i))
    tsort=np.argsort(tnb)
    
    #Get back the results for each thread in a list of results
    for i in np.arange(N_threads) :
        r=result_q.get(i)
        if i == 0 : dum = [r]
        else : dum.append(r)

    #Reorder data from each thread into output matrix
    for i in tsort :
        if i == tsort[0] : outmat = dum[i]
        else : outmat=np.ma.concatenate((outmat,dum[i]),0) if isinstance(outmat,np.ma.masked_array) else np.concatenate((outmat,dum[i]),0)
    
    if len(outmat) != len(outmat) :
        raise '[ERROR]Output array is not coherent with input array - check array reconstruction'
    return(outmat)