Healpy: Von Daten zu Healpix Karte

Question

Ich habe ein Datenraster, wo die Zeilen Theta (0, Pi) und die Spalten repräsentieren Phi (0, 2 * Pi) und wo f (Theta, Phi) ist die Dichte der dunklen Materie an diesem Ort. Ich wollte das Leistungsspektrum dafür berechnen und habe mich dazu entschieden, healty zu verwenden.

Was ich nicht verstehen kann ist, wie meine Daten für healpy zu verwenden formatieren. Wenn jemand Code zur Verfügung stellen könnte (aus offensichtlichen Gründen in Python) oder mich auf ein Tutorial verweisen würde, wäre das großartig! Ich habe meine Hand versucht, es mit dem folgenden Code zu tun:

#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 

Aber, wenn ich versuche, den Code auszuführen, um das Leistungsspektrum zu tun. Insbesondere:

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

ich diesen Fehler:

Typeerror: schlecht Anzahl der Pixel

Wenn jemand mir ein gutes Tutorial zeigen könnte oder helfen, meinen Code zu bearbeiten, das wäre toll.

Weitere Spezifikationen: meine Daten sind in einem 2d NP-Array und ich kann die numRows/numCols ändern, wenn ich muss.

Edit:

Ich habe dieses Problem, indem zuerst die args von anafast healpix_map Ändern gelöst. Ich verbesserte auch den Abstand, indem ich meine Nows * Ncols = 12 * nside * nside machte. Aber mein Leistungsspektrum gibt immer noch Fehler. Wenn jemand Links zu einer guten Dokumentation/Anleitung zur Berechnung des Leistungsspektrums (Zustand von Theta/Phi Args) hat, wäre das unglaublich hilfreich.

Answer 1

Da gehen Sie, hoffe es ist, was Sie suchen. Fühlen Sie sich frei mit Fragen Kommentar :)

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) 

Da die Karte oben nichts als Rauschen enthält, sollte das Leistungsspektrum nur Schrotrauschen enthalten, das heißt flach sein.

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