code: def natural_neighbor_to_grid(xp, yp, variable, grid_x, grid_y):
      r
      # 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)
docs: 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