我有一个数据网格,其中的行表示theta(0,pi),列表示phi(0,2 * pi),其中f(theta,phi)是该位置暗物质的密度。我想为此计算功率谱,并决定使用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的条件)的好的文档/教程,那将非常有用。

最佳答案

您去了那里,希望它是您想要的。随意发表评论:)

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)

python - Healpy:从数据到Healpix map-LMLPHP

由于上面的 map 仅包含噪声,因此功率谱应仅包含散粒噪声,即是平坦的。
# Get the power spectrum
Cl = hp.anafast(hpxmap)
plt.figure()
plt.plot(Cl)

python - Healpy:从数据到Healpix map-LMLPHP

关于python - Healpy:从数据到Healpix map ,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/31573572/

10-10 18:21
查看更多