Healpy: от Data to Healpix map

Question

У меня есть сетка данных, где строки представляют тета (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, если нужно.

Edit:

Я решил эту проблему, сначала изменив арг из anafast в healpix_map. Я также улучшил интервал, сделав свои Nrows * Ncols = 12 * nside * nside. Но мой спектр мощности все еще дает ошибки. Если у кого есть ссылки на хорошую документацию/учебник о том, как рассчитать спектр мощности (состояние тета/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) 
hpxmap[indices] += fs[indices] 

# Inspect the map 
hp.mollview(hpxmap) 

Поскольку выше карта не содержит ничего, кроме шума, спектр мощности должен содержать только дробовой шум, т.е. быть плоской.

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