def natural_neighbor_to_grid(xp, yp, variable, grid_x, grid_y):
    r"""Generate a natural neighbor interpolation of the given points to a regular grid.

    This assigns values to the given grid using the Liang and Hale [Liang2010]_.
    approach.

    Parameters
    ----------
    xp: (N, ) ndarray
        x-coordinates of observations
    yp: (N, ) ndarray
        y-coordinates of observations
    variable: (N, ) ndarray
        observation values associated with (xp, yp) pairs.
        IE, variable[i] is a unique observation at (xp[i], yp[i])
    grid_x: (M, 2) ndarray
        Meshgrid associated with x dimension
    grid_y: (M, 2) ndarray
        Meshgrid associated with y dimension

    Returns
    -------
    img: (M, N) ndarray
        Interpolated values on a 2-dimensional grid

    See Also
    --------
    natural_neighbor_to_points

    """
    # Handle grid-to-points conversion, and use function from `interpolation`
    points_obs = list(zip(xp, yp))
    points_grid = generate_grid_coords(grid_x, grid_y)
    img = natural_neighbor_to_points(points_obs, variable, points_grid)
    return img.reshape(grid_x.shape)
r"""Generate a natural neighbor interpolation of the given points to a regular grid.

    This assigns values to the given grid using the Liang and Hale [Liang2010]_.
    approach.

    Parameters
    ----------
    xp: (N, ) ndarray
        x-coordinates of observations
    yp: (N, ) ndarray
        y-coordinates of observations
    variable: (N, ) ndarray
        observation values associated with (xp, yp) pairs.
        IE, variable[i] is a unique observation at (xp[i], yp[i])
    grid_x: (M, 2) ndarray
        Meshgrid associated with x dimension
    grid_y: (M, 2) ndarray
        Meshgrid associated with y dimension

    Returns
    -------
    img: (M, N) ndarray
        Interpolated values on a 2-dimensional grid

    See Also
    --------
    natural_neighbor_to_points