Healpy query_polygon returns incorrect area

I'm using healpy.query_polygon to find the pixel indices corresponding to a rectangular FOV. However, the indices which are returned do not match my input- creating a polygon in healpy space which is ~100x larger than the FOV (or returns ~14000 pixels instead of the expected ~30 pixels).

The query_disc function works as I expect, however, this is not the function I'm looking to use.

The corresponding outputs: Top(disc), bottom (polygon)

For hp.query_disc:

disc_center = astropy.coordinates.spherical_to_cartesian(1, np.deg2rad(47.3901), np.deg2rad(319.6428))
#(<Quantity 0.5158909828269898>, <Quantity -0.4383926760027675>, <Quantity 0.7359812195055898>)
radius = 0.06208
qd = hp.query_disc(128,disc_center,radius)
#plotting
for ind in range(len(prob)):
        if ind in qd:
            prob[ind] = 1
        else:
            prob[ind] = 0

For hp.query_polygon:

ra_poly, dec_poly = (array([ 48.51458308,  48.51458308,  46.20781856,  46.20781856]), array([ 317.00403838,  322.28167591,  317.11703852,  322.16867577]))
xyzpoly = astropy.coordinates.spherical_to_cartesian(1, np.deg2rad(dec_poly), np.deg2rad(ra_poly))
#xyzpoly = array([[ 0.48450204,  0.52400015,  0.50709248,  0.5465906 ],
   [-0.45174162, -0.40526109, -0.47093848, -0.42445795],
   [ 0.74912435,  0.74912435,  0.72185467,  0.72185467]])

qp = hp.query_polygon(128,xyzpoly)
#plotting
for ind in range(len(prob)):
        if ind in qp:
            prob[ind] = 1
        else:
            prob[ind] = 0

Can anyone explain this discrepancy? By all accounts I don't see any errors, unless the vertices are being implemented incorrectly into the function.

Solution

Your problem is easily solved by considering what query_polygon accepts as input: it wants an (N, 3) array of coordinates, not an (3, N). The value you're passing it is a 3-tuple of 4 elements each, which is transformed into an array of shape (3, 4). Not sure what query_polygon does with the fourth component (possibly ignoring it), but otherwise it will match a triangle (on the sphere).

The solution is straightforward: transpose the coordinate array first.

>>> qp = hp.query_polygon(128, array(xyzpoly).T)
>>> len(qp)
36

There are two problems with your example code:

swap RA and dec. spherical_to_cartesian will complain
swap the last two values for the RA. Healpy will crash Python(!) with the message that your polygon is not convex (connect the dots and you'll see).

An additional problem is that your polygon doesn't define what is inside, and what is outside. It looks like healpy does the "right" thing, and assumes the smallest area is "inside". Defining the polygon the "other way around" (reverse the spherical coordinates) yields the same result.