4

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

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

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

谢谢，

亚历克斯

4

3 回答 3

3

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

解决方案 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)

它给出了一个 2D numpy 数组。我很想知道如何将其转换回 healpix 地图...

于 2014-10-31T18:37:04.473 回答
3

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

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

基本上，您将healpy.mollview( gnomview,cartview和orthviewwork as well ) 的reproject_to_healpix功能与reproject包中的功能 ( http://reproject.readthedocs.org/en/stable/ ) 结合起来。

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

-----基本大纲----------

步骤 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 步：编写模板全天 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 notebook 格式 + FITS 样本数据）------

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

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

------之前和之后的图像------

NSIDE 128，天体坐标：*之前*

NSIDE 128，银河坐标：*之后*

于 2016-04-25T05:26:30.227 回答
1

这个函数似乎可以解决问题（相当慢，但应该比 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)))

在此处输入图像描述

于 2017-08-08T04:02:34.827 回答