在healpy中对HEALPix地图应用旋转

5 投票

3 回答

4495 浏览

数据工程师

提问于 2025-04-18 12:31

我有一个HEALPix地图，是用healpy读取的，不过它是以银河坐标表示的，我需要把它转换成天球/赤道坐标。有没有人知道简单的方法可以转换这个地图？

我尝试过用healpy.Rotator把坐标从(l,b)转换成(phi,theta)，然后用healpy.ang2pix来重新排列像素，但地图看起来还是很奇怪。

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

谢谢，

Alex

3 个回答

这个函数看起来能解决问题（虽然速度比较慢，但应该比用循环要好）：

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)))

回答于 2025-04-18 由 Python大师

分享举报

我找到了一种可能的解决方案，经过几个月的反复搜索。虽然我还没有进行太多测试，所以请小心使用！

上面Saul提到的解决方案2是关键（真是个好建议！）

基本上，你需要把healpy.mollview（gnomview、cartview和orthview也可以）和reproject包里的reproject_to_healpix功能结合起来（http://reproject.readthedocs.org/en/stable/）。

这样得到的地图适合我的角度尺度，但我不能确定这个转换和其他方法相比有多准确。

-----基本步骤----------

步骤1：读取地图并通过cartview生成矩形数组。正如Saul上面提到的，这也是一种旋转的方法。如果你只是进行标准的旋转/坐标转换，那么只需要使用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：写一个模板的全天图FITS头文件（如下例所示）。我写的头文件使得像素尺度的平均值与我想要的HEALPix地图相同。

步骤3：使用reproject.transform_to_healpix

reproject包含一个将“普通”数组（或FITS文件）映射到HEALPix投影的功能。把这个功能和healpy.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笔记本格式 + FITS示例数据）------

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

那里的代码运行起来应该很快，但这是因为地图的质量被大幅降低。我对我的NSIDE 1024和2048的地图做了同样的处理，花了大约一个小时。

------前后对比图------

回答于 2025-04-18 由 Python大师

分享举报

我知道这个问题早就有人问过了，但我这周也遇到了同样的问题，刚好看到了你的帖子。我找到了一些可能的解决办法，所以分享出来，希望能帮到其他人。

解决办法 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数组。我想知道怎么把这个转换回healpix地图...

回答于 2025-04-18 由 Python大师

分享举报

在healpy中对HEALPix地图应用旋转

3 个回答

撰写回答