query_disc misses pixels

[talensgj]

Hi,

I've encountered an issue where using query_disc with inclusive=True does not return all required pixels. I encountered this issue near the pole. Increasing the fact parameter does not seem to solve this issue. The code example below illustrates that when querying a 9.9 degree disc centered on lon = 45 deg and lat = -84 deg in a grid with nside=4, pixel 190 is not returned when it should be.

import numpy as np
import healpy as hp

import matplotlib.pyplot as plt
from matplotlib.patches import Circle

print(hp.__version__)

nside = 4

# Query pixels contributing to a 9.9 deg disk at lon,lat = 45,-84
vec = hp.ang2vec(45, -84, lonlat=True)
pixels = hp.query_disc(nside, vec, np.radians(9.9), inclusive=True)
print(pixels)

# Plot the returned pixels.
npix = hp.nside2npix(nside)
array = np.zeros(npix)
array[pixels] = 1

hp.orthview(array, rot=(45, -84, 0), half_sky=True)
c = Circle((0, 0), np.sin(np.radians(9.9)), facecolor='None', edgecolor='k')
plt.gca().add_patch(c)
plt.show()

# Try with a high fact parameter.
pixels = hp.query_disc(nside, vec, np.radians(9.9), inclusive=True, fact=1024)
print(pixels)

npix = hp.nside2npix(nside)
array = np.zeros(npix)
array[pixels] = 1

hp.orthview(array, rot=(45, -84, 0), half_sky=True)
c = Circle((0, 0), np.sin(np.radians(9.9)), facecolor='None', edgecolor='k')
plt.gca().add_patch(c)
plt.show()

Output: 1.17.3 [180 181 182 187 188 189 191] [180 181 182 187 188 189 191]

[mreineck]

This reminds me of #904. I'm travelling at the moment, but will have a look as soon as possible.

[mreineck]

I have a tentative fix, but I'm not yet terribly happy with it and hope to come up with something more elegant soon(ish). The fix will have to happen in healpix_cxx, so it's not something that we can address with a quick patch to the healpy sources. Is this a time-critical issue for you, @talensgj?

[talensgj]

I've managed to work around it by increasing the search radius for now, so it's not urgent. Thanks for looking into it so promptly.