def interpolate_bilinear_lonlat(lon, lat, values, order='ring'):
    """
    Interpolate values at specific longitudes/latitudes using bilinear interpolation

    Parameters
    ----------
    lon, lat : :class:`~astropy.units.Quantity`
        The longitude and latitude values as :class:`~astropy.units.Quantity` instances
        with angle units.
    values : `~numpy.ndarray`
        Array with the values in each HEALPix pixel. The first dimension should
        have length 12 * nside ** 2 (and nside is determined automatically from
        this).
    order : { 'nested' | 'ring' }
        Order of HEALPix pixels

    Returns
    -------
    result : float `~numpy.ndarray`
        The interpolated values
    """
    nside = npix_to_nside(values.shape[0])
    indices, weights = bilinear_interpolation_weights(lon, lat, nside, order=order)
    values = values[indices]
    # At this point values has shape (N, M) where both N and M might be several
    # dimensions, and weights has shape (N,), so we need to transpose in order
    # to benefit from broadcasting, then transpose back so that the dimension
    # with length 4 is at the start again, ready for summing.
    result = (values.T * weights.T).T
    return result.sum(axis=0)