Healpy 教程

Persus & Xie

简介

healpy是一个Python包，用于处理球体上的像素化数据。它基于分层等面积孤立像素化(HEALPix)方案，并捆绑了HEALPix c++库。

HEALPix的开发是为了有效地处理来自BOOMERANG和WMAP等宇宙学实验的宇宙微波背景数据，但它现在被用于天体物理学的其他分支，以存储来自全天观测的数据。目标受众过去主要是宇宙学科学界，但目前任何对处理球面像素化数据感兴趣的人都非常欢迎提出新的特性。

Healpy提供实用程序:

• 在HEALPix嵌套和环形方案中天空坐标和像素索引之间的转换
• 在天空中寻找磁盘、多边形或条带中的像素
• 在银河参考系、黄道参考系和赤道参考系之间应用坐标变换
• 应用自定义旋转的矢量或完整的地图
• 读取和写入HEALPix映射到磁盘的FITS格式
• 升级和降低现有HEALPix映射的分辨率
• 在Mollweide、Gnomonic和Cartographic投影中可视化地图
• 使用多线程C++例程将映射转换为球面调和空间
• 从地图中计算自动和交叉功率谱，并从光谱创建地图实现

所需其他模块

在数据处理以及可视化时需要用到两个模块：matplotlib和numpy。

import matplotlib.pyplot as plt
import numpy as np
import healpy as hp

NSIDE和排序

地图就是简单的numpy数组，其中每个数组元素都指向由Healpix像素化方案(参见Healpix网站)定义的天空中的一个位置。

映射的分辨率由NSIDE参数定义，该参数通常是2的幂，比如：$2^2=4,2^3=8,2^4=16,2^5=32,\cdot \cdot \cdot$。

程序实例1

import healpy as hp

NSIDE = 32
print(
    "Approximate resolution at NSIDE {} is {:.2} deg".format(
        NSIDE, hp.nside2resol(NSIDE, arcmin=True) / 60
    )
)

输出

Approximate resolution at NSIDE 32 is 1.8 deg

程序实例2

函数healpy.pixelfunc.nside2npix可以给出map中的像素数目NPIX。

import healpy as hp

NSIDE = 32
NPIX = hp.nside2npix(NSIDE)
print(NPIX)

输出

12288

程序实例3

import matplotlib.pyplot as plt
import numpy as np
import healpy as hp

NSIDE = 1
NPIX = hp.nside2npix(NSIDE)
m = np.arange(NPIX)
hp.mollview(m, title="Mollview image RING")
hp.graticule()
plt.show()

输出为

0.0 180.0 -180.0 180.0

为了方便了解数据的排序方式，对图中的颜色块进行标识序号：

其坐标系中的纬度为余纬$\theta$。余纬指“纬度的余”，即：余纬+纬度=$90^{\circ }$。北极纬度为0，$\pi /2$为赤道，$\pi$为南极。经度$\varphi$由东向西从0递增到$2\pi$。在Mollview投影中，$\varphi$在地图的中心，向左为东。

程序实例4

也可以用矢量表示坐标，例如vec是一个归一化矢量，其指向：$\theta =\pi /2,\varphi =3/4\pi$

import healpy as hp
vec = hp.ang2vec(np.pi / 2, np.pi * 3 / 4)
print(vec)

输出为：

[-7.07106781e-01 7.07106781e-01 6.12323400e-17]

程序实例5

在里面找到某个点10度以内所有像素的索引，然后在map上更改这些索引的值：

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt

NSIDE = 16
NPIX = hp.nside2npix(NSIDE)
vec = hp.ang2vec(np.pi / 2, np.pi * 3 / 4)
ipix_disc = hp.query_disc(nside=NSIDE, vec=vec, radius=np.radians(10))
m = np.arange(NPIX)
m[ipix_disc] = m.max()
hp.mollview(m, title="Mollview image RING")
plt.show()

结果图例

程序实例6

我们可以使用pix2ang检索每个像素的纬度和经度，在这种情况下，可以注意到前4个像素中心覆盖了北极~1.46度，这四个点都在同一纬度。第五个像素已经是另一个像素环的一部分

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt

NSIDE = 16
NPIX = hp.nside2npix(NSIDE)
vec = hp.ang2vec(np.pi / 2, np.pi * 3 / 4)
ipix_disc = hp.query_disc(nside=NSIDE, vec=vec, radius=np.radians(10))
m = np.arange(NPIX)
m[ipix_disc] = m.max()
theta, phi = np.degrees(hp.pix2ang(nside=32, ipix=[0, 1, 2, 3, 4]))
print("the latitude:",theta)
print("the longitude:",phi)

输出为

the latitude: [1.46197116 1.46197116 1.46197116 1.46197116 2.92418036]

the longitude: [ 45. 135. 225. 315. 22.5]

程序实例7

RING排序对于球面谐波变换是必要的，另一个选项是嵌套排序，它对于映射域操作非常有效，因为向上和向下缩放映射只需要乘以和舍入像素索引。看看下面的像素是如何在嵌套方案中排序的，它是在12个HEALPix基本像素的结构(NSIDE 1)中再分割生成的:

NSIDE = 2
NPIX = hp.nside2npix(NSIDE)
m = np.arange(NPIX)
hp.mollview(m, nest=True, title="Mollview image NESTED")
hp.graticule()
plt.savefig("test3.png",bbox_inches='tight',dpi=300,format="png")
plt.show()

结果显示

为了方便了解数据的排序方式，对部分图中的颜色块进行标识序号：

所有healpy例程都假定RING顺序，事实上，只要你用read_map读取一个map，即使它被存储为NESTED，它也会被转换为RING。但是，对于大多数healpy routines可以按照嵌套顺序将nest=True`参数传递。

未完待续······