9

我有一个数据网格，其中行代表 theta (0, pi)，列代表 phi (0, 2*pi)，其中 f(theta,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

但是，当我尝试执行代码来计算功率谱时。具体来说， ：

cl = hp.anafast(healpix_map[pix], lmax=1024)

我收到此错误：

TypeError：错误的像素数

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

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

编辑：

我通过首先将anafast的参数更改为healpix_map解决了这个问题。我还通过设置 Nrows*Ncols=12*nside*nside 来改进间距。但是，我的功率谱仍然出现错误。如果有人有关于如何计算功率谱（theta/phi args 的条件）的良好文档/教程的链接，那将非常有帮助。

4

2 回答 2

3

你去吧，希望它是你正在寻找的东西。随时评论问题:)

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)

在此处输入图像描述

由于上面的地图只包含噪声，因此功率谱应该只包含散粒噪声，即平坦。

# Get the power spectrum
Cl = hp.anafast(hpxmap)
plt.figure()
plt.plot(Cl)

在此处输入图像描述

于 2017-10-22T04:20:11.613 回答
3

numpy.add.at按照这个答案，有一种更快的方法可以使用 进行地图初始化。

与丹尼尔出色答案的第一部分相比，这在我的机器上要快几倍：

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)
于 2019-10-29T01:15:22.333 回答