Healpy query_polygon 返回不正确的区域答案

【问题标题】：Healpy query_polygon returns incorrect areaHealpy query_polygon 返回不正确的区域
【发布时间】：2016-06-13 22:39:31
【问题描述】：

我正在使用healpy.query_polygon 来查找对应于矩形FOV 的像素索引。但是，返回的索引与我的输入不匹配 - 在 healpy 空间中创建一个比 FOV 大约 100 倍的多边形（或返回约 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

谁能解释这种差异？总的来说，我没有看到任何错误，除非顶点被错误地实现到函数中。

【问题讨论】：

如果您不混淆 RA 和 dec，它可能会有所帮助；复制粘贴第二个代码块的前两行会导致错误（虽然很容易修复）。请在发布时验证您的示例代码，以避免混淆。

标签： python healpy

【解决方案1】：

您的问题很容易通过考虑query_polygon 接受的输入来解决：它需要(N, 3) 坐标数组，而不是(3, N)。您传递给它的值是一个包含 4 个元素的 3 元组，它被转换为一个形状为 (3, 4) 的数组。不确定query_polygon 对第四个组件做了什么（可能忽略它），否则它将匹配一个三角形（在球体上）。

解决方法很简单：先转置坐标数组。

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

您的示例代码有两个问题：

交换 RA 和 dec。 spherical_to_cartesian会抱怨
交换 RA 的最后两个值。 Healpy 将导致 Python(!) 崩溃，并显示您的多边形不是凸面的消息（连接点，您会看到）。

另一个问题是您的多边形没有定义内部和外部。看起来 healpy 做了“正确”的事情，并假设最小的区域是“内部”。以“其他方式”（反转球坐标）定义多边形会产生相同的结果。

【讨论】：