Healpy:从数据到 Healpix 地图 [英] Healpy: From Data to Healpix map

问题描述

我有一个数据网格，其中行代表 theta (0, pi)，列代表 phi (0, 2*pi)，其中 f(theta,phi) 是该位置的暗物质密度.我想为此计算功率谱并决定使用healpy.

I have a data grid where the rows represent theta (0, pi) and the columns represent phi (0, 2*pi) and where f(theta,phi) is the density of dark matter at that location. I wanted to calculate the power spectrum for this and have decided to use healpy.

我无法理解的是如何格式化我的数据以供健康使用.如果有人可以提供代码(出于显而易见的原因在 python 中)或指向我的教程，那就太好了！我已尝试使用以下代码进行操作:

What I can not understand is how to format my data for healpy to use. If someone could provide code (in python for obvious reasons) or point me to a tutorial, that would be great! I have tried my hand at doing it with the following code:

#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 

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

But, when I try to execute the code to do the power spectrum. Specifically, :

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

我收到此错误:

类型错误:像素数错误

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

If anyone could point me to a good tutorial or help edit my code that would be great.

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

More specifications: my data is in a 2d np array and I can change the numRows/numCols if I need to.

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

I have solved this problem by first changing the args of anafast to healpix_map. I also improved the spacing by making my Nrows*Ncols=12*nside*nside. But, my power spectrum is still giving errors. If anyone has links to good documentation/tutorial on how to calculate the power spectrum (condition of theta/phi args), that would be incredibly helpful.

推荐答案

有一种使用 numpy.add.at 进行地图初始化的更快方法，遵循 这个答案.

There is a faster way to do the map initialization using numpy.add.at, following this answer.

与丹尼尔的优秀回答的第一部分相比，这在我的机器上快了数倍:

This is several times faster on my machine as compared to the first section of 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)

这篇关于Healpy:从数据到 Healpix 地图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！