Healpy query_polygon возвращает неправильную область

Question

Я использую healpy.query_polygon для поиска индексов пикселей, соответствующих прямоугольному FOV. Однако возвращаемые индексы не совпадают с моими вводами - создание многоугольника в healpy space, которое на 100x больше, чем FOV (или возвращает ~ 14000 пикселей вместо ожидаемых ~ 30 пикселей).

Функция query_disc работает так, как я ожидаю, однако это не функция, которую я ищу использовать.

Соответствующие выходы: Top(disc), bottom (polygon)

Для 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 

Для 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 

Может кто-нибудь объяснить это несоответствие? По всем учетным записям я не вижу никаких ошибок, если вершины не будут внедрены неправильно в функцию.

Answer 1

Ваша проблема легко решена путем учета того, что query_polygon принимает в качестве входных данных: ему нужен массив координат , а не (3, N). Значение, которое вы передаете, представляет собой 3-элементный 4 элемента каждый, который преобразуется в массив формы (3, 4). Не уверен, что query_polygon делает с четвертым компонентом (возможно, игнорируя его), но в противном случае он будет соответствовать треугольнику (на сфере).

Решение прост: сначала перенесите сначала координатный массив.

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

---

Есть две проблемы с вашим примером кода:

• своп РА и Декабре spherical_to_cartesian будет жаловаться

• обменять последние два значения для RA. Healpy приведет к сбою Python (!) С сообщением о том, что ваш многоугольник не выпуклый (соедините точки, и вы увидите).

Дополнительная проблема заключается в том, что ваш многоугольник не определяет, что внутри, и что снаружи. Похоже, что healpy делает «правильную» вещь и предполагает, что самая маленькая область «внутри». Определение полигона «наоборот» (обратное сферическим координатам) дает тот же результат.