astropy / astropy-healpix

BSD-licensed HEALPix for Astropy - maintained by @astrofrog and @lpsinger
https://astropy-healpix.readthedocs.io
BSD 3-Clause "New" or "Revised" License
50 stars 22 forks source link

Why the results of cone_search_lonlat() are different from healpy.query_disc()? #213

Closed StarCycle closed 4 months ago

StarCycle commented 4 months ago

Hi,

I use both cone_search_lonlat() and healpy.query_disc() to find which pixel is located in the Field of View of a telescope with the following code (the field angle is 1°):

import numpy as np
import healpy as hp
from astropy import units as u
from astropy_healpix import HEALPix
from utils import get_args

nside = 32
radius = 1

npix = hp.nside2npix(nside)
pixels = np.arange(npix)
ra, dec = hp.pix2ang(nside, pixels, lonlat=True)

pixels_in_FOV = []
probs_in_FOV = np.zeros(npix)

healpix = HEALPix(nside=nside, order='nested')
for i in range(npix):
    index_nested = healpix.cone_search_lonlat(ra[i] * u.deg, dec[i] * u.deg, radius=radius * u.deg)
    index_ring = hp.nest2ring(nside, index_nested)
    pixels_in_FOV.append(index_ring)

pixels_in_FOV_ = []
for i in range(npix):
    index = hp.query_disc(nside, hp.ang2vec(ra[i], dec[i], True), radius=radius * np.pi / 180, inclusive=True)
    pixels_in_FOV_.append(index)

But the results from these two functions are different:

pixels_in_FoV[0]: array([0,4,5])  # cone_search_lonlat()
pixels_in_FoV_[0]: array([0,1,3,4,5]) # query_disc()

Why are they different?

mreineck commented 4 months ago

My guess is that cone_search_lonlat only returns pixels whose centers lie inside the cone; query_disc will do exactly the same if you provide inclusive=False. For inclusive=True it will return (perhaps with a few false positives) all pixels that partially overlap the cone, even if their center is not in the cone.

lpsinger commented 4 months ago

Yes, @mreineck is correct.