Healpy query_polygon RuntimeError: Unknown exception

Question

I'm using healpy.query_polygon to get a list of healpix indexes inside a polygon. According to the documentation:

vertices: Vertex array containing the vertices of the polygon, shape (N, 3).

But when I try to get all indexes from the following polygon, a RuntimeError: Unknown exception appears:

In[1]:

import healpy as hp

vertex_array = np.array([[0.65, -0.04, 0.76], [0.58, 0.38, 0.72], [0.91, -0.29, 0.31],[0.91, 0.18, 0.38]])

print(vertex_array.shape)
vertex_array

Out[1]:

(4, 3)
array([[ 0.65, -0.04,  0.76],
       [ 0.58,  0.38,  0.72],
       [ 0.91, -0.29,  0.31],
       [ 0.91,  0.18,  0.38]])

In[2]:

healpix_indexes_test = hp.query_polygon(4, vertex_array)

Out[2]:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-63-5a14f69cb078> in <module>
----> 1 healpix_indexes_test = hp.query_polygon(4, vertex_array)

healpy/src/_query_disc.pyx in healpy._query_disc.query_polygon()

RuntimeError: Unknown exception

Here you can see visualisation of these dots lying on a sphere.

Just for fun I tried to transpose the input array, so its shape became (3, 4). The Unknown exception problem disappeared. But input like this contradict the documentation, so I don't believe it.

In[1]:

print(vertex_array.T.shape)
vertex_array.T

Out[1]:

(3, 4)
array([[ 0.65,  0.58,  0.91,  0.91],
       [-0.04,  0.38, -0.29,  0.18],
       [ 0.76,  0.72,  0.31,  0.38]])

In[2]:

healpix_indexes_test_1 = hp.query_polygon(4, vertex_array.T)

healpix_indexes_test_1

Out[2]:

array([ 42,  58,  75, 107, 123, 140])

I'll be gratefull for any suggestions.

score 3 · Answer 1 · answered Apr 08 '19 at 13:24

With the help of Healpy members on Github, the solution was found. It's important to define the order of the vertices correctly. In my case it means that the last two dot's coordinates should be swopped to make the rectangle simple, not self-crossing:

In[1]:

vertex_array_fixed = np.array([[0.65, -0.04, 0.76], [0.58, 0.38, 0.72], [0.91, 0.18, 0.38], [0.91, -0.29, 0.31]])

print(vertex_array_fixed.shape)
vertex_array_fixed

Out[1]:

(4, 3)
array([[ 0.65, -0.04,  0.76],
       [ 0.58,  0.38,  0.72],
       [ 0.91,  0.18,  0.38],
       [ 0.91, -0.29,  0.31]])

In[2]:

healpix_indexes_test = hp.query_polygon(4, vertex_array_fixed)

healpix_indexes_test

Out[2]:

array([24, 40, 71])

Here comes the visualisation:

In[3]:

Nside = 2048

healpix_indexes_test = hp.query_polygon(Nside, vertex_array_fixed)

healpix_indexes_test

Out[3]:

array([ 5968575,  5968576,  5968577, ..., 17387119, 17387120, 17395310])

In[4]:

Npix = hp.nside2npix(Nside)

whole_map = np.arange(Npix, dtype=float)
whole_map[healpix_indexes_test] = hp.UNSEEN

hp.mollview(m, title="Fixed rectangle")

Out[4]: Output plot

Healpy query_polygon RuntimeError: Unknown exception

1 Answers1