def interpolate_to_grid(x, y, z, interp_type='linear', hres=50000,
                        minimum_neighbors=3, gamma=0.25, kappa_star=5.052,
                        search_radius=None, rbf_func='linear', rbf_smooth=0,
                        boundary_coords=None):
    r"""Interpolate given (x,y), observation (z) pairs to a grid based on given parameters.

    Parameters
    ----------
    x: array_like
        x coordinate
    y: array_like
        y coordinate
    z: array_like
        observation value
    interp_type: str
        What type of interpolation to use. Available options include:
        1) "linear", "nearest", "cubic", or "rbf" from `scipy.interpolate`.
        2) "natural_neighbor", "barnes", or "cressman" from `metpy.interpolate`.
        Default "linear".
    hres: float
        The horizontal resolution of the generated grid, given in the same units as the
        x and y parameters. Default 50000.
    minimum_neighbors: int
        Minimum number of neighbors needed to perform barnes or cressman interpolation for a
        point. Default is 3.
    gamma: float
        Adjustable smoothing parameter for the barnes interpolation. Default 0.25.
    kappa_star: float
        Response parameter for barnes interpolation, specified nondimensionally
        in terms of the Nyquist. Default 5.052
    search_radius: float
        A search radius to use for the barnes and cressman interpolation schemes.
        If search_radius is not specified, it will default to the average spacing of
        observations.
    rbf_func: str
        Specifies which function to use for Rbf interpolation.
        Options include: 'multiquadric', 'inverse', 'gaussian', 'linear', 'cubic',
        'quintic', and 'thin_plate'. Defualt 'linear'. See `scipy.interpolate.Rbf` for more
        information.
    rbf_smooth: float
        Smoothing value applied to rbf interpolation.  Higher values result in more smoothing.
    boundary_coords: dictionary
        Optional dictionary containing coordinates of the study area boundary. Dictionary
        should be in format: {'west': west, 'south': south, 'east': east, 'north': north}

    Returns
    -------
    grid_x: (N, 2) ndarray
        Meshgrid for the resulting interpolation in the x dimension
    grid_y: (N, 2) ndarray
        Meshgrid for the resulting interpolation in the y dimension ndarray
    img: (M, N) ndarray
        2-dimensional array representing the interpolated values for each grid.

    Notes
    -----
    This function acts as a wrapper for `interpolate_points` to allow it to generate a regular
    grid.

    See Also
    --------
    interpolate_to_points

    """
    # Generate the grid
    if boundary_coords is None:
        boundary_coords = get_boundary_coords(x, y)
    grid_x, grid_y = generate_grid(hres, boundary_coords)

    # Handle grid-to-points conversion, and use function from `interpolation`
    points_obs = np.array(list(zip(x, y)))
    points_grid = generate_grid_coords(grid_x, grid_y)
    img = interpolate_to_points(points_obs, z, points_grid, interp_type=interp_type,
                                minimum_neighbors=minimum_neighbors, gamma=gamma,
                                kappa_star=kappa_star, search_radius=search_radius,
                                rbf_func=rbf_func, rbf_smooth=rbf_smooth)

    return grid_x, grid_y, img.reshape(grid_x.shape)
r"""Interpolate given (x,y), observation (z) pairs to a grid based on given parameters.

    Parameters
    ----------
    x: array_like
        x coordinate
    y: array_like
        y coordinate
    z: array_like
        observation value
    interp_type: str
        What type of interpolation to use. Available options include:
        1) "linear", "nearest", "cubic", or "rbf" from `scipy.interpolate`.
        2) "natural_neighbor", "barnes", or "cressman" from `metpy.interpolate`.
        Default "linear".
    hres: float
        The horizontal resolution of the generated grid, given in the same units as the
        x and y parameters. Default 50000.
    minimum_neighbors: int
        Minimum number of neighbors needed to perform barnes or cressman interpolation for a
        point. Default is 3.
    gamma: float
        Adjustable smoothing parameter for the barnes interpolation. Default 0.25.
    kappa_star: float
        Response parameter for barnes interpolation, specified nondimensionally
        in terms of the Nyquist. Default 5.052
    search_radius: float
        A search radius to use for the barnes and cressman interpolation schemes.
        If search_radius is not specified, it will default to the average spacing of
        observations.
    rbf_func: str
        Specifies which function to use for Rbf interpolation.
        Options include: 'multiquadric', 'inverse', 'gaussian', 'linear', 'cubic',
        'quintic', and 'thin_plate'. Defualt 'linear'. See `scipy.interpolate.Rbf` for more
        information.
    rbf_smooth: float
        Smoothing value applied to rbf interpolation.  Higher values result in more smoothing.
    boundary_coords: dictionary
        Optional dictionary containing coordinates of the study area boundary. Dictionary
        should be in format: {'west': west, 'south': south, 'east': east, 'north': north}

    Returns
    -------
    grid_x: (N, 2) ndarray
        Meshgrid for the resulting interpolation in the x dimension
    grid_y: (N, 2) ndarray
        Meshgrid for the resulting interpolation in the y dimension ndarray
    img: (M, N) ndarray
        2-dimensional array representing the interpolated values for each grid.

    Notes
    -----
    This function acts as a wrapper for `interpolate_points` to allow it to generate a regular
    grid.

    See Also
    --------
    interpolate_to_points