Healpy: From Data to Healpix Card

Question

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 this point. I wanted to calculate the power spectrum for this and decided to use healpy.

I cannot figure out how to format my data to use healpy. If someone can provide the code (in python for obvious reasons) or point me to a tutorial, that would be great! I tried my hand by doing this 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 make a power spectrum. In particular:

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

I get this error:

TypeError: bad pixel count

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

Additional features: my data is in a 2d np array, and I can change numRows / numCols if I need to.

Edit:

I solved this problem by first changing the anaphast arguments to healpix_map. I also improved the interval by making my Nrows * Ncols = 12 * nside * nside. But my power spectrum still gives errors. If anyone has links to good documentation / tutorial on how to calculate the power spectrum (theta / phi args state), this would be incredibly useful.

Answer 1

That is all, hopefully this is what you are looking for. Feel free to comment with questions :)

 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) 

Since there is nothing but noise on the map above, the power spectrum should only contain a shot, i.e. be flat.

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