在healpy中对HEALPix贴图应用旋转

Question

我有一个HEALPix地图，我在使用healpy时读到的，但是它是在银河系坐标中，我需要在天体/赤道坐标中使用它。有没有人知道转换地图的简单方法？

我尝试过使用healpy.Rotator将(l，b)转换为(phi，theta)，然后使用healpy.ang2pix重新排序像素，但地图看起来仍然很奇怪。

如果有一个类似于Rotator的函数，你可以像这样调用：map = AnotherRotator(map,coord=['G','C'])，那就太好了。有人知道有这样的函数吗？？

谢谢,

亚历克斯

Answer 1

我意识到这个问题很久以前就被问过了，但是这周我自己也遇到了同样的问题，我找到了你的帖子。我已经找到了几个潜在的解决方案，所以我会分享给其他人，以防其他人发现它很有用。

解决方案1：这在某种程度上取决于数据的格式。我的是一个(theta，phi)网格。

import numpy as np
import healpy as H
map = <your original map>
nside = <your map resolution, mine=256>
npix = H.nside2npix(nside)
pix = N.arange(npix)
t,p = H.pix2ang(nside,pix) #theta, phi

r = H.Rotator(deg=True, rot=[<THETA ROTATION>, <PHI ROTATION>])

map_rot = np.zeros(npix)

for i in pix:
    trot, prot = r(t[i],p[i])
    tpix = int(trot*180./np.pi) #my data came in a theta, phi grid -- this finds its location there
    ppix = int(prot*180./np.pi)
    map_rot[i] = map[ppix,tpix] #this being the rright way round may need double-checking

解决方案2：还没有完成测试，但在做了上面烦人的工作后才发现它……

map_rot = H.mollview(map,deg=True,rot=[<THETA>,<PHI>], return_projected_map=True)

这给出了一个二维numpy数组。我很有兴趣知道如何将它转换回一个健康地图...

Answer 2

在断断续续地寻找了几个月后，我找到了另一个潜在的解决方案。我还没有过多的测试，所以请多加小心！

上面索尔的Solution 2是关键(很棒的建议！)

基本上，您将healpy.mollview (也可以使用gnomview、cartview和orthview )的功能与reproject包(http://reproject.readthedocs.org/en/stable/)中的reproject_to_healpix函数的功能结合起来。

生成的地图适合我的角度比例，但我不能说与其他方法相比，转换有多精确。

-基本大纲-

步骤1:通过cartview读入映射并生成矩形数组。正如索尔上面指出的，这也是一种进行旋转的方法。如果你只是做一个标准的旋转/坐标变换，那么你所需要的就是coord关键字。从天体坐标到银河系坐标，设置coord = ['C','G']

map_Gal = hp.cartview(map_Cel, coord=['C','G'], return_projected_map=True, xsize=desired_xsize, norm='hist',nest=False)

第2步：编写一个模板all-sky FITS header (如下例所示)。我写的地图的平均像素比例与我想要的HEALPix地图相同。

第3步：使用reproject.transform_to_healpix

reproject包含一个将“普通”数组(或FITS文件)映射到HEALPix投影的函数。将它与返回由HEALPix py.Mollview/cartview/orthview/gnomview创建的数组的能力结合起来，您可以将一个坐标系(天球)的HEALPix地图旋转到另一个坐标系(银河)。

map_Gal_HP, footprint_Gal_HP = rp.reproject_to_healpix((map_Gal, target_header), coord_system_out= 'GALACTIC', nside=nside, nested=False)

从本质上讲，它归结为这两个命令。然而，你必须制作一个模板标题，给出与你想要制作的中间全天地图相对应的像素比例和大小。

-完整工作示例(iPython notebook格式+ FITS示例数据)-

https://github.com/aaroncnb/healpix_coordtrans_example/tree/master

那里的代码应该运行得非常快，但这是因为地图已经严重降级。我为我的NSIDE 1024和2048地图做了同样的事情，花了大约一个小时。

-前后镜像-

​

​

​

​

Answer 3

这个函数似乎做到了这一点(速度相当慢，但应该比for循环更好)：

def rotate_map(hmap, rot_theta, rot_phi):
    """
    Take hmap (a healpix map array) and return another healpix map array 
    which is ordered such that it has been rotated in (theta, phi) by the 
    amounts given.
    """
    nside = hp.npix2nside(len(hmap))

    # Get theta, phi for non-rotated map
    t,p = hp.pix2ang(nside, np.arange(hp.nside2npix(nside))) #theta, phi

    # Define a rotator
    r = hp.Rotator(deg=False, rot=[rot_phi,rot_theta])

    # Get theta, phi under rotated co-ordinates
    trot, prot = r(t,p)

    # Interpolate map onto these co-ordinates
    rot_map = hp.get_interp_val(hmap, trot, prot)

    return rot_map

对来自PyGSM的数据使用此方法可获得以下结果：

hp.mollview(np.log(rotate_map(gsm.generated_map_data, 0,0)))

​

​

在旋转phi时

hp.mollview(np.log(rotate_map(gsm.generated_map_data, 0,np.pi)))

​

​

或旋转theta

hp.mollview(np.log(rotate_map(gsm.generated_map_data, np.pi/4,0)))

​

​