Healpy：从数据到Healpix映射

Question

我有一个数据网格，其中的行表示θ（0，π），列表示φ（0，2*pi），其中f（θ，phi）是该位置暗物质的密度。我想计算这个的功率谱，并决定使用healpy。在

我不明白的是如何格式化我的数据供healpy使用。如果有人能提供代码（在Python中有明显的原因）或者指向我的教程，那就太好了！我试过用以下代码来实现：

#grid dimensions are Nrows*Ncols (subject to change)
theta = np.linspace(0, np.pi, num=grid.shape[0])[:, None]
phi = np.linspace(0, 2*np.pi, num=grid.shape[1])
nside = 512
print "Pixel area: %.2f square degrees" % hp.nside2pixarea(nside, degrees=True)
pix = hp.ang2pix(nside, theta, phi)
healpix_map = np.zeros(hp.nside2npix(nside), dtype=np.double)
healpix_map[pix] = grid 

但是，当我试图执行代码来做功率谱的时候。具体来说：

^{pr2}$

我得到这个错误：

类型错误：像素数错误

如果有人能给我指出一个好的教程或帮助编辑我的代码，那就太好了。在

更多规格： 我的数据在2dnp数组中，如果需要，我可以更改numRows/numCols。在

编辑：

我首先将anafast的参数改为healpix_映射，从而解决了这个问题。 我还通过使Nrows*Ncols=12*nside*nside来改进间距。 但是，我的功率谱仍然有误差。如果任何人有关于如何计算功率谱（theta/phi args的条件）的好的文档/教程的链接，那将是非常有用的。在

Answer 1

在this answer之后，有一种使用numpy.add.at进行映射初始化的更快方法。在

这在我的机器上比Daniel's excellent answer的第一部分快几倍：

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

# Set the number of sources and the coordinates for the input
nsources = int(1e7)
nside = 64
npix = hp.nside2npix(nside)

# Coordinates and the density field f
thetas = np.random.uniform(0, np.pi, nsources)
phis = np.random.uniform(0, 2*np.pi, nsources)
fs = np.random.randn(nsources)

# Go from HEALPix coordinates to indices
indices = hp.ang2pix(nside, thetas, phis)

# Baseline, from Daniel Lenz's answer:
# time: ~5 s
hpxmap1 = np.zeros(npix, dtype=np.float)
for i in range(nsources):
    hpxmap1[indices[i]] += fs[i]

# Using numpy.add.at
# time: ~0.6 ms
hpxmap2 = np.zeros(npix, dtype=np.float)
np.add.at(hpxmap2, indices, fs)

Answer 2

给你，希望这就是你要找的。欢迎提出问题：）

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

# Set the number of sources and the coordinates for the input
nsources = int(1.e4)
nside = 16
npix = hp.nside2npix(nside)

# Coordinates and the density field f
thetas = np.random.random(nsources) * np.pi
phis = np.random.random(nsources) * np.pi * 2.
fs = np.random.randn(nsources)

# Go from HEALPix coordinates to indices
indices = hp.ang2pix(nside, thetas, phis)

# Initate the map and fill it with the values
hpxmap = np.zeros(npix, dtype=np.float)
for i in range(nsources):
    hpxmap[indices[i]] += fs[i]

# Inspect the map
hp.mollview(hpxmap)

由于上面的图只包含噪声，所以功率谱应该只包含散粒噪声，即平坦的。在

^{pr2}$