Healpy:从Data到Healpix地图

Question

我有一个数据网格,其中行代表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的args更改为healpix_map来解决这个问题.我还通过使我的Nrows*Ncols = 12*nside*nside来改善间距.但是,我的功率谱仍然存在错误.如果有人有关于如何计算功率谱的良好文档/教程的链接(theta/phi args的条件),那将是非常有帮助的.

Answer 1

就这样，希望这是您正在寻找的。有问题请随时评论:)

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)