BUG: neighbours() produces warnings for many cases

[taldcroft]

In [1]: import astropy_healpix as hp

In [2]: hp.neighbours(480, nside=16, order="ring")
/Users/aldcroft/miniconda3-arm/envs/ska3/lib/python3.11/site-packages/astropy_healpix/core.py:659: RuntimeWarning: invalid value encountered in neighbours_ring
  return np.stack(func(healpix_index, nside))
Out[2]: array([544,  -1, 543, 479, 420, 481, 545, 608])

In [3]: hp.__version__
Out[3]: '1.0.2'

More generally:

In [8]: import astropy_healpix as hp

In [9]: npix = hp.nside_to_npix(16)

In [10]: npix
Out[10]: 3072

In [12]: import warnings

In [13]: for idx in range(npix):
    ...:     with warnings.catch_warnings(record=True) as warns:
    ...:         hp.neighbours(idx, nside=16)
    ...:         if len(warns) > 0:
    ...:             print(f"warning for {idx}")
warning for 480
warning for 495
warning for 496
warning for 511
warning for 512
warning for 527
warning for 528
warning for 543
warning for 544
warning for 560
warning for 576
warning for 592
warning for 2464
warning for 2480
warning for 2496
warning for 2512
warning for 2528
warning for 2543
warning for 2544
warning for 2559
warning for 2560
warning for 2575
warning for 2576
warning for 2591

In [14]: for idx in range(npix):
    ...:     with warnings.catch_warnings(record=True) as warns:
    ...:         hp.neighbours(idx, nside=16, order="nested")
    ...:         if len(warns) > 0:
    ...:             print(f"warning for {idx}")
warning for 85
warning for 170
warning for 341
warning for 426
warning for 597
warning for 682
warning for 853
warning for 938
warning for 1024
warning for 1279
warning for 1280
warning for 1535
warning for 1536
warning for 1791
warning for 1792
warning for 2047
warning for 2133
warning for 2218
warning for 2389
warning for 2474
warning for 2645
warning for 2730
warning for 2901
warning for 2986

[lpsinger]

Is this a bug? The return value does really contain invalid pixel indices (-1) for pixels that are near the poles and do not have the full complement of neighbors.

[mreineck]

To be precise, this happens at the locations where the 12 base pixels meet and that are not on the poles or the equator. For Nside>1, at most one of the neighbours can be -1 (i.e. nonexistant), for Nside==1, two of the neighbours may be absent.

[Edit: by "meet" I mean at the corners of the base pixels.]

[mreineck]

An easier interpretation: these warnings happen at pixels that touch corners where exactly 3 base pixels meet. Every Healpix map has exactly 8 of those locations, and each is touched by three pixels,so you get 24 warnings.

[taldcroft]

I would call this a bug and a documentation issue.

Bug - there is no reason this code should be issuing a scary sounding warning for a situation that is expected and normal. This warning ended up making me not trust the results and I spent an hour tracking down what was happening and submitting this issue.

Documentation - the neighbours() docstring and any examples should highlight that -1 is an expected return value in this case, effectively "Not a Neighbour". Basically I didn't understand or trust astropy_healpix until I read the corresponding healpy documentation and confirmed the -1 return value with that code.

[lpsinger]

I would call this a bug and a documentation issue.

I agree. Proposed fix in #220.

Bug - there is no reason this code should be issuing a scary sounding warning for a situation that is expected and normal. This warning ended up making me not trust the results and I spent an hour tracking down what was happening and submitting this issue.

It would be easy enough to remove the warning; we would just remove this line: https://github.com/astropy/astropy-healpix/blob/v1.0.3/astropy_healpix/_core.c#L302.

However, that would make this method inconsistent with the others; all of the functions set the invalid status if they return any invalid pixel values.