python - Healpy: From Data to Healpix map -

i have data grid rows represent theta (0, pi) , columns represent phi (0, 2*pi) , f(theta,phi) density of dark matter @ location. wanted calculate power spectrum , have decided use healpy.

what can not understand how format data healpy use. if provide code (in python obvious reasons) or point me tutorial, great! have tried hand @ doing following code:

#grid dimensions nrows*ncols (subject 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 try execute code power spectrum. specifically, :

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

i error:

typeerror: bad number of pixels

if point me tutorial or edit code great.

more specifications: data in 2d np array , can change numrows/numcols if need to.

edit:

i have solved problem first changing args of anafast healpix_map. improved spacing making nrows*ncols=12*nside*nside. but, power spectrum still giving errors. if has links documentation/tutorial on how calculate power spectrum (condition of theta/phi args), incredibly helpful.

there go, hope it's you're looking for. feel free comment questions :)

import healpy hp import numpy np import matplotlib.pyplot plt  # set number of sources , coordinates input nsources = int(1.e4) nside = 16 npix = hp.nside2npix(nside)  # coordinates , density field f thetas = np.random.random(nsources) * np.pi phis = np.random.random(nsources) * np.pi * 2. fs = np.random.randn(nsources)  # go healpix coordinates indices indices = hp.ang2pix(nside, thetas, phis)  # initate map , fill values hpxmap = np.zeros(npix, dtype=np.float) hpxmap[indices] += fs[indices]  # inspect map hp.mollview(hpxmap) 

since map above contains nothing noise, power spectrum should contain shot noise, i.e. flat.

# power spectrum cl = hp.anafast(hpxmap) plt.figure() plt.plot(cl)