Python scipy.interpolate.griddata() Examples

def plot_solution(X_star, u_star, index):
    
    lb = X_star.min(0)
    ub = X_star.max(0)
    nn = 200
    x = np.linspace(lb[0], ub[0], nn)
    y = np.linspace(lb[1], ub[1], nn)
    X, Y = np.meshgrid(x,y)
    
    U_star = griddata(X_star, u_star.flatten(), (X, Y), method='cubic')
    
    plt.figure(index)
    plt.pcolor(X,Y,U_star, cmap = 'jet')
    plt.colorbar() 

def test_complex_2d(self):
        x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
                     dtype=np.double)
        y = np.arange(x.shape[0], dtype=np.double)
        y = y - 2j*y[::-1]

        xi = x[:,None,:] + np.array([0,0,0])[None,:,None]

        for method in ('nearest', 'linear', 'cubic'):
            for rescale in (True, False):
                msg = repr((method, rescale))
                yi = griddata(x, y, xi, method=method, rescale=rescale)

                assert_equal(yi.shape, (5, 3), err_msg=msg)
                assert_allclose(yi, np.tile(y[:,None], (1, 3)),
                                atol=1e-14, err_msg=msg) 

def _view_griddata(cls, data, data_view, target_view, method='nearest',
                       x_tolerance=1e-3, y_tolerance=1e-3):
        t_xmin, t_ymin, t_xmax, t_ymax = target_view.bbox
        d_xmin, d_ymin, d_xmax, d_ymax = data_view.bbox
        nodata = target_view.nodata
        view = np.full(target_view.shape, nodata)
        viewcoords = target_view.coords
        datacoords = data_view.coords
        yx_tolerance = np.sqrt(x_tolerance**2 + y_tolerance**2)
        row_bool = ((datacoords[:,0] <= t_ymax + y_tolerance) &
                    (datacoords[:,0] >= t_ymin - y_tolerance))
        col_bool = ((datacoords[:,1] <= t_xmax + x_tolerance) &
                    (datacoords[:,1] >= t_xmin - x_tolerance))
        yx_grid = datacoords[row_bool & col_bool]
        view = interpolate.griddata(yx_grid,
                                    data.flat[row_bool & col_bool],
                                    viewcoords, method=method,
                                    fill_value=nodata)
        view = view.reshape(target_view.shape)
        return view 

def interpolation_near(x1, y1, x2, y2, x1grd, y1grd, method='linear', **kwargs):
    ''' Interpolate values of x2/y2 onto full-res grids of x1/y1 using
    linear interpolation of nearest points
    Parameters
    ----------
        x1 : 1D vector - X coordinates of keypoints on image 1
        y1 : 1D vector - Y coordinates of keypoints on image 1
        x1 : 1D vector - X coordinates of keypoints on image 2
        y1 : 1D vector - Y coordinates of keypoints on image 2
        x1grd : 1D vector - source X coordinate on img1
        y1grd : 1D vector - source Y coordinate on img2
        method : str - parameter for SciPy griddata
    Returns
    -------
        x2grd : 1D vector - destination X coordinate on img1
        y2grd : 1D vector - destination Y coordinate on img2
    '''
    src = np.array([y1, x1]).T
    dst = np.array([y1grd, x1grd]).T
    x2grd = griddata(src, x2, dst, method=method).T
    y2grd = griddata(src, y2, dst, method=method).T

    return x2grd, y2grd 

def interpj(self, key_or_var, lay, lat, angle):
        from scipy.interpolate import griddata
        lays = self.variables['LAY']
        lats = self.variables['LAT']
        angles = self.variables['ANGLE']
        points = np.array([[[(lay_, lat_, angle_) for angle_ in angles]
                            for lat_ in lats] for lay_ in lays]).reshape(-1, 3)
        if isinstance(key_or_var, str):
            v = self.variables[key_or_var].ravel()
        else:
            v = np.array(key_or_var)
            assert(key_or_var.shape == (self.NLAYS, self.NLATS, self.NANGLES))
        if isinstance(lay, (int, float)):
            idx = np.array((lay, lat, angle), ndmin=2)
        else:
            idx = np.array((lay, lat, angle)).swapaxes(0, 1)
            assert(idx.shape[0] == 2)
        return griddata(points, v, idx)[0] 

def regrid(self):
        ''' Regrid the non-uniform data on a regular mesh '''

        if self.values is None:
            warnings.warn('Cannont regrid: data missing.')
            return

        # First we need to interpolate the non-uniformly sampled data
        # onto a uniform grid
        from scipy.interpolate import griddata
        
        x = int(np.sqrt(self.n_points / 2))

        G = np.mgrid[-np.pi:np.pi:2j*x, 0:np.pi:1j*x]
        spherical_grid = np.squeeze(G.reshape((2,1,-1)))

        gridded = griddata(self.spherical.T, self.values, spherical_grid.T, method='nearest', fill_value=0.)
        dirty_img = gridded.reshape((2*x, x))

        return G[0], G[1], gridded.reshape((2*x, x)) 

def get_reference_elevation(input, recGrid, elevation):
    """
    The following function define the elevation from the TIN to a regular grid...

    Args:
        input: class containing XML input file parameters.
        recGrid: class describing the regular grid characteristics.
        elevation: TIN elevation mesh.

    Returns:
        - ref_elev - interpolated elevation on the regular grid
    """
    if input.searef:
        x_ref, y_ref = input.searef
        pts = recGrid.tinMesh["vertices"]
        ref_elev = griddata(
            points=pts, values=elevation, xi=[x_ref, y_ref], method="nearest"
        )
    else:
        ref_elev = 0.0

    return ref_elev 

def plot_landscape(data):
        """
        Plot the data as an energy landscape.

        Args:
            data: (x, y, xx, yy, z, xlim, ylim). x, y, z represent the \
                  coordinates of the points that will be interpolated. xx, yy \
                  represent the meshgrid used to interpolate the points. xlim, \
                  ylim are tuples containing the limits of the x and y axes.

        Returns:
            Plot representing the interpolated energy landscape of the target points.

        """
        x, y, xx, yy, z, xlim, ylim = data
        zz = griddata((x, y), z, (xx, yy), method="linear")
        mesh = holoviews.QuadMesh((xx, yy, zz))
        contour = holoviews.operation.contours(mesh, levels=8)
        scatter = holoviews.Scatter((x, y))
        contour_mesh = mesh * contour * scatter
        return contour_mesh.redim(
            x=holoviews.Dimension("x", range=xlim), y=holoviews.Dimension("y", range=ylim),
        ) 

def test_xi_1d(self):
        # Check that 1-D xi is interpreted as a coordinate
        x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
                     dtype=np.double)
        y = np.arange(x.shape[0], dtype=np.double)
        y = y - 2j*y[::-1]

        xi = np.array([0.5, 0.5])

        for method in ('nearest', 'linear', 'cubic'):
            p1 = griddata(x, y, xi, method=method)
            p2 = griddata(x, y, xi[None,:], method=method)
            assert_allclose(p1, p2, err_msg=method)

            xi1 = np.array([0.5])
            xi3 = np.array([0.5, 0.5, 0.5])
            assert_raises(ValueError, griddata, x, y, xi1,
                          method=method)
            assert_raises(ValueError, griddata, x, y, xi3,
                          method=method) 

def test_square_rescale_manual(self):
        points = np.array([(0,0), (0,100), (10,100), (10,0), (1, 5)], dtype=np.double)
        points_rescaled = np.array([(0,0), (0,1), (1,1), (1,0), (0.1, 0.05)], dtype=np.double)
        values = np.array([1., 2., -3., 5., 9.], dtype=np.double)

        xx, yy = np.broadcast_arrays(np.linspace(0, 10, 14)[:,None],
                                     np.linspace(0, 100, 14)[None,:])
        xx = xx.ravel()
        yy = yy.ravel()
        xi = np.array([xx, yy]).T.copy()

        for method in ('nearest', 'linear', 'cubic'):
            msg = method
            zi = griddata(points_rescaled, values, xi/np.array([10, 100.]),
                          method=method)
            zi_rescaled = griddata(points, values, xi, method=method,
                                   rescale=True)
            assert_allclose(zi, zi_rescaled, err_msg=msg,
                            atol=1e-12) 

def test_1d_borders(self):
        # Test for nearest neighbor case with xi outside
        # the range of the values.
        x = np.array([1, 2.5, 3, 4.5, 5, 6])
        y = np.array([1, 2, 0, 3.9, 2, 1])
        xi = np.array([0.9, 6.5])
        yi_should = np.array([1.0, 1.0])

        method = 'nearest'
        assert_allclose(griddata(x, y, xi,
                                 method=method), yi_should,
                        err_msg=method,
                        atol=1e-14)
        assert_allclose(griddata(x.reshape(6, 1), y, xi,
                                 method=method), yi_should,
                        err_msg=method,
                        atol=1e-14)
        assert_allclose(griddata((x, ), y, (xi, ),
                                 method=method), yi_should,
                        err_msg=method,
                        atol=1e-14) 

def warp_image(im, flow):
    """
    Use optical flow to warp image to the next
    :param im: image to warp
    :param flow: optical flow
    :return: warped image
    """
    from scipy import interpolate
    image_height = im.shape[0]
    image_width = im.shape[1]
    flow_height = flow.shape[0]
    flow_width = flow.shape[1]
    n = image_height * image_width
    (iy, ix) = np.mgrid[0:image_height, 0:image_width]
    (fy, fx) = np.mgrid[0:flow_height, 0:flow_width]
    fx += flow[:,:,0]
    fy += flow[:,:,1]
    mask = np.logical_or(fx <0 , fx > flow_width)
    mask = np.logical_or(mask, fy < 0)
    mask = np.logical_or(mask, fy > flow_height)
    fx = np.minimum(np.maximum(fx, 0), flow_width)
    fy = np.minimum(np.maximum(fy, 0), flow_height)
    points = np.concatenate((ix.reshape(n,1), iy.reshape(n,1)), axis=1)
    xi = np.concatenate((fx.reshape(n, 1), fy.reshape(n,1)), axis=1)
    warp = np.zeros((image_height, image_width, im.shape[2]))
    for i in range(im.shape[2]):
        channel = im[:, :, i]
        plt.imshow(channel, cmap='gray')
        values = channel.reshape(n, 1)
        new_channel = interpolate.griddata(points, values, xi, method='cubic')
        new_channel = np.reshape(new_channel, [flow_height, flow_width])
        new_channel[mask] = 1
        warp[:, :, i] = new_channel.astype(np.uint8)

    return warp.astype(np.uint8) 

def sarpy2ortho(ro, pix, decimation=10):

    nx = ro.sicdmeta.ImageData.FullImage.NumCols
    ny = ro.sicdmeta.ImageData.FullImage.NumRows

    nx_dec = round(nx / decimation)
    ny_dec = round(ny / decimation)

    xv, yv = np.meshgrid(range(nx), range(ny), indexing='xy')
    xv = xv[::decimation, ::decimation]
    yv = yv[::decimation, ::decimation]
    npix = xv.size

    xv = np.reshape(xv, (npix,1))
    yv = np.reshape(yv, (npix,1))
    im_points = np.concatenate([yv,xv], axis=1)

    ground_coords = image_to_ground(im_points, ro.sicdmeta)
    ground_coords = ecf_to_geodetic(ground_coords)

    minx = np.min(ground_coords[:,1])
    maxx = np.max(ground_coords[:,1])
    miny = np.min(ground_coords[:,0])
    maxy = np.max(ground_coords[:,0])

    xi, yi = create_ground_grid(minx, maxx, miny, maxy, nx_dec, ny_dec)

    ground_coords[:,[0,1]] = ground_coords[:,[1,0]]
    pix = np.reshape(pix, npix)
    gridded = griddata(ground_coords[:,0:2], pix, (xi,yi), method='nearest').astype(np.uint8)

    ul = [maxy, minx]
    lr = [miny, maxx]

    extent = [ul,lr]

    extent_poly = bounds_2_shapely_polygon(min_x=minx, max_x=maxx, min_y=miny, max_y=maxy)
    geot = world_poly_to_geo_t(extent_poly, npix_x=nx_dec, npix_y=ny_dec)

    return gridded, extent, geot 

def sarpy2ortho(ro, pix, decimation=10):

    nx = ro.sicdmeta.ImageData.FullImage.NumCols
    ny = ro.sicdmeta.ImageData.FullImage.NumRows

    nx_dec = round(nx / decimation)
    ny_dec = round(ny / decimation)

    xv, yv = np.meshgrid(range(nx), range(ny))
    xv = xv[::decimation, ::decimation]
    yv = yv[::decimation, ::decimation]
    npix = xv.size

    xv = np.reshape(xv, (npix,1))
    yv = np.reshape(yv, (npix,1))
    im_points = np.concatenate([xv,yv], axis=1)

    ground_coords = image_to_ground(im_points, ro.sicdmeta)
    ground_coords = ecf_to_geodetic(ground_coords)

    minx = np.min(ground_coords[:,1])
    maxx = np.max(ground_coords[:,1])
    miny = np.min(ground_coords[:,0])
    maxy = np.max(ground_coords[:,0])

    xi, yi = create_ground_grid(minx, maxx, miny, maxy, nx_dec, ny_dec)

    ground_coords[:,[0,1]] = ground_coords[:,[1,0]]
    pix = np.reshape(pix, npix)
    gridded = griddata(ground_coords[:,0:2], pix, (xi,yi), method='nearest').astype(np.uint8)

    ul = [maxy, minx]
    lr = [miny, maxx]

    extent = [ul,lr]

    extent_poly = bounds_2_shapely_polygon(min_x=minx, max_x=maxx, min_y=miny, max_y=maxy)
    geot = world_poly_to_geo_t(extent_poly, npix_x=nx_dec, npix_y=ny_dec)

    return gridded, extent, geot 

def sarpy2ortho(ro, pix, decimation=10):

    nx = ro.sicdmeta.ImageData.FullImage.NumCols
    ny = ro.sicdmeta.ImageData.FullImage.NumRows

    nx_dec = round(nx / decimation)
    ny_dec = round(ny / decimation)

    xv, yv = np.meshgrid(range(nx), range(ny), indexing='xy')
    xv = xv[::decimation, ::decimation]
    yv = yv[::decimation, ::decimation]
    npix = xv.size

    xv = np.reshape(xv, (npix,1))
    yv = np.reshape(yv, (npix,1))
    im_points = np.concatenate([yv,xv], axis=1)

    ground_coords = image_to_ground(im_points, ro.sicdmeta)
    ground_coords = ecf_to_geodetic(ground_coords)

    minx = np.min(ground_coords[:,1])
    maxx = np.max(ground_coords[:,1])
    miny = np.min(ground_coords[:,0])
    maxy = np.max(ground_coords[:,0])

    xi, yi = create_ground_grid(minx, maxx, miny, maxy, nx_dec, ny_dec)

    ground_coords[:,[0,1]] = ground_coords[:,[1,0]]
    pix = np.reshape(pix, npix)
    gridded = griddata(ground_coords[:,0:2], pix, (xi,yi), method='nearest').astype(np.uint8)

    ul = [maxy, minx]
    lr = [miny, maxx]

    extent = [ul,lr]

    extent_poly = bounds_2_shapely_polygon(min_x=minx, max_x=maxx, min_y=miny, max_y=maxy)
    geot = world_poly_to_geo_t(extent_poly, npix_x=nx_dec, npix_y=ny_dec)

    return gridded, extent, geot 

def interp_map_irregular_grid(a, x, y, x_i, y_i, method="linear", debug=False):
    """Interpolates fields from any grid to another grid
    !!! Careful when using this on regular grids.
    Results are not unique and it takes forever.
    Use interp_map_regular_grid instead!!!
    """
    xx, yy = np.meshgrid(x, y)
    xx_i, yy_i = np.meshgrid(x_i, y_i)

    # pad margins to avoid nans in the interpolation
    xx = xx[:, [-1] + list(range(xx.shape[1])) + [0]]
    xx = xx[[-1] + list(range(xx.shape[0])) + [0], :]

    xx[:, 0] = xx[:, 0] - 360
    xx[:, -1] = xx[:, -1] + 360

    yy = yy[:, [-1] + list(range(yy.shape[1])) + [0]]
    yy = yy[[-1] + list(range(yy.shape[0])) + [0], :]

    yy[0, :] = yy[0, :] - 180
    yy[-1, :] = yy[-1, :] + 180

    a = a[:, [-1] + list(range(a.shape[1])) + [0]]
    a = a[[-1] + list(range(a.shape[0])) + [0], :]

    if debug:
        print("a shape", a.shape)
        print("x shape", xx.shape)
        print("y shape", yy.shape)
        print("x values", xx[0, :])
        print("y values", yy[:, 0])

    points = np.vstack((xx.flatten(), yy.flatten())).T
    values = a.flatten()
    int_points = np.vstack((xx_i.flatten(), yy_i.flatten())).T
    a_new = spi.griddata(points, values, int_points, method=method)

    return a_new.reshape(xx_i.shape) 

def interp_file(Lat,Long,file):
    # Import the DOVPM file
    f = gdal.Open(file)
    # load band (akin to a variable in dataset)
    band = f.GetRasterBand(1)   
    # get the pixel width, height, etc.
    transform = f.GetGeoTransform()
    # Grid resolution (known)
    res=transform[1]
    # convert lat,lon to row,col
    column = (Long - transform[0]) / transform[1]
    row = (Lat - transform[3]) / transform[5]
    # get pixel values surrounding data point
    Surrounding_data=(band.ReadAsArray(np.floor(column-2), np.floor(row-2), 5, 5))
    # convert row,col back to north,east
    Long_c = transform[0] + np.floor(column) * res
    Lat_c = transform[3] - np.floor(row) * res
    # set up matrices for interpolation
    count=-1
    pos=np.zeros((25,2))
    Surrounding_data_v=np.zeros((25,1))
    for k in range(-2,3):
        for j in range(-2,3):
            count=count+1
            pos[count]=(Long_c+j*res,Lat_c-k*res)
            Surrounding_data_v[count]=Surrounding_data[k+2,j+2]          
    interp_val=griddata(pos,Surrounding_data_v,(Long,Lat),method='cubic')
    return interp_val

#___________________________________________________________________________#
# Functions to handle the conversions from one height to another 

def fill_depth(depth):
    """ Fill in the holes in the depth map """
    ht, wd = depth.shape
    x, y = np.meshgrid(np.arange(wd), np.arange(ht))
    xx = x[depth > 0].astype(np.float32)
    yy = y[depth > 0].astype(np.float32)
    zz = depth[depth > 0].ravel()
    return interpolate.griddata((xx, yy), zz, (x, y), method='nearest') 

def load_RSM(filename):
    om, tt, psd = xu.io.getxrdml_map(filename)
    om = np.deg2rad(om)
    tt = np.deg2rad(tt)
    wavelength = 1.54056

    q_y = (1 / wavelength) * (np.cos(tt) - np.cos(2 * om - tt))
    q_x = (1 / wavelength) * (np.sin(tt) - np.sin(2 * om - tt))

    xi = np.linspace(np.min(q_x), np.max(q_x), 100)
    yi = np.linspace(np.min(q_y), np.max(q_y), 100)
    psd[psd < 1] = 1
    data_grid = griddata(
        (q_x, q_y), psd, (xi[None, :], yi[:, None]), fill_value=1, method="cubic"
    )
    nx, ny = data_grid.shape

    range_values = [np.min(q_x), np.max(q_x), np.min(q_y), np.max(q_y)]
    output_data = (
        Panel(np.log(data_grid).reshape(nx, ny, 1), minor_axis=["RSM"])
        .transpose(2, 0, 1)
        .to_frame()
    )

    return range_values, output_data 

def test_contour_heatmap():
    from scipy.interpolate import griddata
    from matplotlib import pyplot as plt

    mesh3D = mesh(200)
    mesh2D = proj_to_2D(mesh3D)

    data = np.zeros((3,3))
    data[0,1] += 2

    vals = np.exp(log_dirichlet_density(mesh3D,2.,data=data.sum(0)))
    temp = log_censored_dirichlet_density(mesh3D,2.,data=data)
    censored_vals = np.exp(temp - temp.max())

    xi = np.linspace(-1,1,1000)
    yi = np.linspace(-0.5,1,1000)

    plt.figure()
    plt.contour(griddata((mesh2D[:,0],mesh2D[:,1]),vals,(xi[None,:],yi[:,None]), method='cubic'))
    plt.axis('off')
    plt.title('uncensored likelihood')

    plt.figure()
    plt.contour(griddata((mesh2D[:,0],mesh2D[:,1]),censored_vals,(xi[None,:],yi[:,None]),method='cubic'))
    plt.axis('off')
    plt.title('censored likelihood') 

def fill_depth(depth):
    x, y = np.meshgrid(np.arange(depth.shape[1]).astype("float32"),
                       np.arange(depth.shape[0]).astype("float32"))
    xx = x[depth > 0]
    yy = y[depth > 0]
    zz = depth[depth > 0]

    grid = interpolate.griddata((xx, yy), zz.ravel(),
                                (x, y), method='nearest')
    return grid 

def fill_depth(depth):
    """ Fill in the holes in the depth map """
    ht, wd = depth.shape
    x, y = np.meshgrid(np.arange(wd), np.arange(ht))
    xx = x[depth > 0].astype(np.float32)
    yy = y[depth > 0].astype(np.float32)
    zz = depth[depth > 0].ravel()
    return interpolate.griddata((xx, yy), zz, (x, y), method='nearest') 

def fill_depth(depth):
    """ Fill in the holes in the depth map """
    ht, wd = depth.shape
    x, y = np.meshgrid(np.arange(wd), np.arange(ht))
    xx = x[depth > 0].astype(np.float32)
    yy = y[depth > 0].astype(np.float32)
    zz = depth[depth > 0].ravel()
    return interpolate.griddata((xx, yy), zz, (x, y), method='linear') 

def nearest_griddata(x, y, z, xi, yi):
    """
    Nearest Neighbor Interpolation Method.

    Nearest-neighbor interpolation (also known as proximal interpolation or, in some contexts, point sampling) is a simple method of multivariate interpolation in one or more dimensions.<br/>

    Interpolation is the problem of approximating the value of a function for a non-given point in some space when given the value of that function in points around (neighboring) that point.<br/>
    The nearest neighbor algorithm selects the value of the nearest point and does not consider the values of neighboring points at all, yielding a piecewise-constant interpolant. <br/>
    The algorithm is very simple to implement and is commonly used (usually along with mipmapping) in real-time 3D rendering to select color values for a textured surface.<br/>

    zi = nearest_griddata(x,y,z,xi,yi) fits a surface of the form z = f*(*x, y) to the data in the (usually) nonuniformly spaced vectors (x, y, z).<br/>
    griddata() interpolates this surface at the points specified by (xi, yi) to produce zi. xi and yi must describe a regular grid.<br/>

    Parameters
    ----------
    x:  array-like
        x-coord [1D array]
    y:  array-like
        y-coord [1D array]
    z:  array-like
        z-coord [1D array]
    xi: array-like
        meshgrid for x-coords [2D array] see <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.meshgrid.html">numpy.meshgrid</a>
    yi: array-like
        meshgrid for y-coords [2D array] see <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.meshgrid.html">numpy.meshgrid</a>

    Returns
    -------
    Interpolated Zi Coord

    zi: array-like
        zi interpolated-value [2D array]  for (xi,yi)
    """
    zi = griddata(zip(x,y), z, (xi, yi), method='nearest')
    return zi 

def warp_image(im, flow):
    """
    Use optical flow to warp image to the next
    :param im: image to warp
    :param flow: optical flow
    :return: warped image
    """
    from scipy import interpolate
    image_height = im.shape[0]
    image_width = im.shape[1]
    flow_height = flow.shape[0]
    flow_width = flow.shape[1]
    n = image_height * image_width
    (iy, ix) = np.mgrid[0:image_height, 0:image_width]
    (fy, fx) = np.mgrid[0:flow_height, 0:flow_width].astype(float)
    fx += flow[:,:,0]
    fy += flow[:,:,1]
    mask = np.logical_or(fx <0 , fx > flow_width)
    mask = np.logical_or(mask, fy < 0)
    mask = np.logical_or(mask, fy > flow_height)
    fx = np.minimum(np.maximum(fx, 0), flow_width)
    fy = np.minimum(np.maximum(fy, 0), flow_height)
    points = np.concatenate((ix.reshape(n,1), iy.reshape(n,1)), axis=1)
    xi = np.concatenate((fx.reshape(n, 1), fy.reshape(n,1)), axis=1)
    warp = np.zeros((image_height, image_width, im.shape[2]))
    for i in range(im.shape[2]):
        channel = im[:, :, i]
        plt.imshow(channel, cmap='gray')
        values = channel.reshape(n, 1)
        new_channel = interpolate.griddata(points, values, xi, method='linear')
        new_channel = np.reshape(new_channel, [flow_height, flow_width])
        new_channel[mask] = 1
        warp[:, :, i] = new_channel.astype(np.uint8)

    return warp.astype(np.uint8) 

def warp_image(im, flow):
    """
    Use optical flow to warp image to the next
    :param im: image to warp
    :param flow: optical flow
    :return: warped image
    """
    from scipy import interpolate
    image_height = im.shape[0]
    image_width = im.shape[1]
    flow_height = flow.shape[0]
    flow_width = flow.shape[1]
    n = image_height * image_width
    (iy, ix) = np.mgrid[0:image_height, 0:image_width]
    (fy, fx) = np.mgrid[0:flow_height, 0:flow_width]
    fx += flow[:,:,0]
    fy += flow[:,:,1]
    mask = np.logical_or(fx <0 , fx > flow_width)
    mask = np.logical_or(mask, fy < 0)
    mask = np.logical_or(mask, fy > flow_height)
    fx = np.minimum(np.maximum(fx, 0), flow_width)
    fy = np.minimum(np.maximum(fy, 0), flow_height)
    points = np.concatenate((ix.reshape(n,1), iy.reshape(n,1)), axis=1)
    xi = np.concatenate((fx.reshape(n, 1), fy.reshape(n,1)), axis=1)
    warp = np.zeros((image_height, image_width, im.shape[2]))
    for i in range(im.shape[2]):
        channel = im[:, :, i]
        plt.imshow(channel, cmap='gray')
        values = channel.reshape(n, 1)
        new_channel = interpolate.griddata(points, values, xi, method='cubic')
        new_channel = np.reshape(new_channel, [flow_height, flow_width])
        new_channel[mask] = 1
        warp[:, :, i] = new_channel.astype(np.uint8)

    return warp.astype(np.uint8) 

def nlloc_project_grid(self):
        """
        Projecting the grid to the new coordinate system.

        This function also determines the 3D grid from the 2D grids from
        NonLinLoc

        """

        # Generating the correct NonLinLoc Formatted Grid
        if self.NLLoc_proj == "NONE":
            GRID_NLLOC = Grid3D(cell_count=self.NLLoc_n,
                                cell_size=self.NLLoc_size,
                                azimuth=0.0,
                                dip=0.0)
            GRID_NLLOC.nlloc_grid_centre(self.NLLoc_org[0], self.NLLoc_org[1])
        else:
            GRID_NLLOC = Grid3D(cell_count=self.NLLoc_n,
                                cell_size=self.NLLoc_size,
                                azimuth=self.NLLoc_MapOrg[2],
                                dip=0.0)
            GRID_NLLOC.lonlat_centre(self.NLLoc_MapOrg[0],
                                     self.NLLoc_MapOrg[1])

        # TO-DO: What is the text in NLLoc_MapOrg[3]?
        if self.NLLoc_proj == "LAMBERT":
            GRID_NLLOC.projections(grid_proj_type=self.NLLoc_MapOrg[3])

        if self.NLLoc_proj == "TRANS_MERC":
            GRID_NLLOC.projections(grid_proj_type=self.NLLoc_MapOrg[3])

        OrgX, OrgY, OrgZ = GRID_NLLOC.grid_xyz
        NewX, NewY, NewZ = self.grid_xyz

        self.NLLoc_data = griddata((OrgX.flatten(),
                                    OrgY.flatten(),
                                    OrgZ.flatten()),
                                   self.NLLoc_data.flatten(),
                                   (NewX, NewY, NewZ),
                                   method="nearest") 

def warp_image(im, flow):
    """
    Use optical flow to warp image to the next
    :param im: image to warp
    :param flow: optical flow
    :return: warped image
    """
    from scipy import interpolate
    image_height = im.shape[0]
    image_width = im.shape[1]
    flow_height = flow.shape[0]
    flow_width = flow.shape[1]
    n = image_height * image_width
    (iy, ix) = np.mgrid[0:image_height, 0:image_width]
    (fy, fx) = np.mgrid[0:flow_height, 0:flow_width]
    fx += flow[:,:,0]
    fy += flow[:,:,1]
    mask = np.logical_or(fx <0 , fx > flow_width)
    mask = np.logical_or(mask, fy < 0)
    mask = np.logical_or(mask, fy > flow_height)
    fx = np.minimum(np.maximum(fx, 0), flow_width)
    fy = np.minimum(np.maximum(fy, 0), flow_height)
    points = np.concatenate((ix.reshape(n,1), iy.reshape(n,1)), axis=1)
    xi = np.concatenate((fx.reshape(n, 1), fy.reshape(n,1)), axis=1)
    warp = np.zeros((image_height, image_width, im.shape[2]))
    for i in range(im.shape[2]):
        channel = im[:, :, i]
        plt.imshow(channel, cmap='gray')
        values = channel.reshape(n, 1)
        new_channel = interpolate.griddata(points, values, xi, method='cubic')
        new_channel = np.reshape(new_channel, [flow_height, flow_width])
        new_channel[mask] = 1
        warp[:, :, i] = new_channel.astype(np.uint8)

    return warp.astype(np.uint8) 

def warp_image(im, flow):
    """
    Use optical flow to warp image to the next
    :param im: image to warp
    :param flow: optical flow
    :return: warped image
    """
    from scipy import interpolate
    image_height = im.shape[0]
    image_width = im.shape[1]
    flow_height = flow.shape[0]
    flow_width = flow.shape[1]
    n = image_height * image_width
    (iy, ix) = np.mgrid[0:image_height, 0:image_width]
    (fy, fx) = np.mgrid[0:flow_height, 0:flow_width]
    # print(flow.shape, flow.shape)
    fx = fx + flow[:,:,0]
    fy = fy + flow[:,:,1]
    mask = np.logical_or(fx <0 , fx > flow_width)
    mask = np.logical_or(mask, fy < 0)
    mask = np.logical_or(mask, fy > flow_height)
    fx = np.minimum(np.maximum(fx, 0), flow_width)
    fy = np.minimum(np.maximum(fy, 0), flow_height)
    points = np.concatenate((ix.reshape(n,1), iy.reshape(n,1)), axis=1)
    xi = np.concatenate((fx.reshape(n, 1), fy.reshape(n,1)), axis=1)
    warp = np.zeros((image_height, image_width, im.shape[2]))
    for i in range(im.shape[2]):
        channel = im[:, :, i]
        plt.imshow(channel, cmap='gray')
        values = channel.reshape(n, 1)
        new_channel = interpolate.griddata(points, values, xi, method='cubic')
        new_channel = np.reshape(new_channel, [flow_height, flow_width])
        new_channel[mask] = 1
        warp[:, :, i] = new_channel.astype(np.uint8)

    return warp.astype(np.uint8) 

def test_1d_unsorted(self):
        x = np.array([2.5, 1, 4.5, 5, 6, 3])
        y = np.array([1, 2, 0, 3.9, 2, 1])

        for method in ('nearest', 'linear', 'cubic'):
            assert_allclose(griddata(x, y, x, method=method), y,
                            err_msg=method, atol=1e-10)
            assert_allclose(griddata(x.reshape(6, 1), y, x, method=method), y,
                            err_msg=method, atol=1e-10)
            assert_allclose(griddata((x,), y, (x,), method=method), y,
                            err_msg=method, atol=1e-10)