在healpy中对HEALPix贴图应用旋转

Question

我有一张HEALPix地图，我已经用healpy读过了，但是它是在银河系的坐标系中，我需要它在天体/赤道坐标系下。有人知道转换地图的简单方法吗？在

我尝试过使用healpy.Rotator将（l，b）转换为（phi，theta），然后使用healpy.ang2pix对像素进行重新排序，但是地图看起来还是很奇怪。在

如果有一个类似于Rotator的函数，你可以这样调用它：map = AnotherRotator(map,coord=['G','C'])。有人知道这种功能吗？？在

谢谢

亚历克斯

Answer 1

我知道这个问题很久以前就有人问过了，但我本周也遇到了同样的问题，找到了你的帖子。我已经找到了几个潜在的解决方案，所以我将分享，以防有人发现它有用。在

解决方案1：这取决于数据的格式。我的是一个（θ，φ）网格。在

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：还没有完全完成测试，只是在完成了上面这些烦人的工作之后才发现的。。。在

^{pr2}$

它提供了一个2D numpy数组。我很想知道如何将这个转换回healpix地图。。。在

Answer 2

这个函数似乎能做到这一点（相当慢，但应该比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)))

Answer 3

我断断续续地找了几个月，找到了另一个可能的解决办法。我还没有做过很多测试，所以请小心！在

索尔的建议是很好的解决办法在

基本上，您将healpy.mollview（gnomview，cartview，和{}的功能）与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标题（如下例所示）。我写我的有相同的平均像素比例作为我想要的HEALPix地图。在

第3步：使用reproject.transform_to_healpix

reproject包含一个函数，用于将“normal”数组（或FITS文件）映射到HEALPix投影中。再加上返回healpy.mollview/cartview/orthview/gnomview，你可以将一个坐标系（天体）的HEALPix地图旋转到另一个坐标系（银河系）。在

^{pr2}$

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

-完整的工作示例（iPython笔记本格式+FITS示例数据）

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

那里的代码应该运行得很快，但那是因为映射严重退化。我对我的nside1024和2048地图也做了同样的工作，而且 大约一个小时。在

图像前后

{a4}