我有一个数据网格,其中的行表示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)
由于上面的 map 仅包含噪声,因此功率谱应仅包含散粒噪声,即是平坦的。
# Get the power spectrum
Cl = hp.anafast(hpxmap)
plt.figure()
plt.plot(Cl)
关于python - Healpy:从数据到Healpix map ,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/31573572/