Neighbours in interpolate_bilinear_lonlat/skycoord methods

[mroldrobot]

Hi all,
could somebody please give me a hint with how interpolate_bilinear_... methods
actually work.  The comments say that it takes 4 neighbours to make interpolation, but it's not clear for me which pixels are the neighbours.  From  _core.c  I found that the four nearest neighbours are probably meant.  So, I tried to do some straightforward test:

from astropy.coordinates import ICRS
from astropy_healpix import HEALPix
from astropy_healpix.core import bilinear_interpolation_weights

NSIDE = 8

hp = HEALPix(nside=NSIDE, frame=ICRS()) # ring order
test_coord = SkyCoord(214*u.deg, 49*u.deg, frame=ICRS())
pix_num =  hp.skycoord_to_healpix(test_coord)

print('Test point falls inside:', pix_num)

pix_neighbours = hp.neighbours(pix_num)
print('All neighbours of the test pixel:', pix_neighbours)

sorted_neighbours = np.argsort(test_coord.separation(hp.healpix_to_skycoord(pix_neighbours)))
print('Sorted by distance:', pix_neighbours[sorted_neighbours])
print('Four nearest by distance:', pix_num, *pix_neighbours[sorted_neighbours][:3])

bilinear_neighbours, _ = bilinear_interpolation_weights(test_coord.ra, test_coord.dec, NSIDE) # ring order
print('Neighbours by bilinear interpolator:', *bilinear_neighbours)

it prints:

Test point falls inside: 100
All neighbours of the test pixel: [130  99  73  51  74 101 131 163]
Sorted by distance: [ 74 101 131 130  73  99  51 163]
Four nearest by distance: 100 74 101 131
Neighbours by bilinear interpolator: 100 101 74 73

so  the nearest pixels by means of distance are not equal to nearest pixels, returned by  bilinear_interpolation method. So what are these neighbours? May be these are the nearest
pixels by means of their boundaries?

Another yet not clear statement is "If a position does not have four neighbours, this currently returns NaN."  I guess it means if one (or more) of the four nearest neighbours has an invalid value. What are invalid values? e.g. in healpy, an invalid value is set by UNSEEN constant. I think it should be np.NaN, but it would be nice to document this somewhere...

Thanks in advance for any help!