2

我healpy.query_polygon用来获取多边形内的 healpix 索引列表。根据文档：

vertices：顶点数组，包含多边形的顶点，形状为 (N, 3)。

但是当我尝试从以下多边形中获取所有索引时，会RuntimeError: Unknown exception出现一个：

在[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

输出[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]])

在[2]：

healpix_indexes_test = hp.query_polygon(4, vertex_array)

输出[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

在这里，您可以看到这些点位于球体上的可视化效果。

只是为了好玩，我尝试转置输入数组，所以它的形状变成了 (3, 4)。Unknown exception问题消失了。但是这样的输入与文档相矛盾，所以我不相信。

在[1]：

print(vertex_array.T.shape)
vertex_array.T

输出[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]])

在[2]：

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

healpix_indexes_test_1

输出[2]：

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

我将不胜感激任何建议。

4

1 回答 1

3

在 Github 上 Healpy 成员的帮助下，找到了解决方案。正确定义顶点的顺序很重要。在我的情况下，这意味着应该交换最后两个点的坐标以使矩形简单，而不是自交叉：

在[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

输出[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]])

在[2]：

healpix_indexes_test = hp.query_polygon(4, vertex_array_fixed)

healpix_indexes_test

输出[2]：

array([24, 40, 71])

这是可视化：

在[3]中：

Nside = 2048

healpix_indexes_test = hp.query_polygon(Nside, vertex_array_fixed)

healpix_indexes_test

输出[3]：

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

在[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]： 输出图

于 2019-04-08T13:24:13.970 回答