<p>您必须避免循环,并在底层的C代码中、在Numpy或Astropy中进行繁重的计算,以便进行天空/像素转换。有几个选项可以用<code>astropy.wcs</code>来实现这一点。你知道吗</p>
<p>第一个是用<code>SkyCoord</code>。我们首先为像素索引创建一个值网格:</p>
<pre class="lang-py prettyprint-override"><code>In [30]: xx, yy = np.mgrid[:5, :5]
...: xx, yy
Out[30]:
(array([[0, 0, 0, 0, 0],
[1, 1, 1, 1, 1],
[2, 2, 2, 2, 2],
[3, 3, 3, 3, 3],
[4, 4, 4, 4, 4]]), array([[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4]]))
</code></pre>
<p>现在我们可以从像素索引创建<code>SkyCoord</code>对象(这是Numpy数组子类),并使用wcs:</p>
<pre class="lang-py prettyprint-override"><code>In [33]: from astropy.coordinates import SkyCoord
...: sky = SkyCoord.from_pixel(xx, yy, wcs)
...: sky
Out[33]:
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
[[(53.17127889, -27.78771333), (53.17127889, -27.78765778),
(53.17127889, -27.78760222), (53.17127889, -27.78754667),
(53.17127889, -27.78749111)],
....
</code></pre>
<p>注意,这是使用<code>wcs.utils.skycoord_to_pixel</code>。这个对象还有一个方法,可以用wcs投影到像素。出于实际目的,我在这里也将这样做:</p>
<pre class="lang-py prettyprint-override"><code>In [34]: sky.to_pixel(wcs)
Out[34]:
(array([[ 0.00000000e+00, -1.11022302e-16, -2.22044605e-16,
-3.33066907e-16, 1.13149046e-10],
...
[ 4.00000000e+00, 4.00000000e+00, 4.00000000e+00,
4.00000000e+00, 4.00000000e+00]]),
array([[-6.31503738e-11, 1.00000000e+00, 2.00000000e+00,
3.00000000e+00, 4.00000000e+00],
...
[-1.11457732e-10, 1.00000000e+00, 2.00000000e+00,
3.00000000e+00, 4.00000000e+00]]))
</code></pre>
<p>我们得到一个新的x和y索引的浮点值元组。因此,您需要对这些值进行舍入,并将其转换为int以用作数组索引。你知道吗</p>
<p>第二个选项是使用较低级别的函数,例如<code>wcs.pixel_to_world_values</code>和<code>wcs.world_to_pixel_values</code>,它接受Nx2数组并返回:</p>
<pre class="lang-py prettyprint-override"><code>In [37]: wcs.pixel_to_world_values(np.array([xx.ravel(), yy.ravel()]).T)
Out[37]:
array([[ 53.17127889, -27.78771333],
[ 53.17127889, -27.78765778],
[ 53.17127889, -27.78760222],
[ 53.17127889, -27.78754667],
...
</code></pre>