Tengo una cuadrícula de datos donde las filas representan theta (0, pi) y las columnas representan phi (0, 2 * pi) y donde f (theta, phi) es la densidad de la materia oscura en esa ubicación. Quería calcular el espectro de potencia para esto y he decidido usar healpy.

Lo que no puedo entender es cómo formatear mis datos para que healpy los use. Si alguien pudiera proporcionar código (en python por razones obvias) o señalarme un tutorial, ¡sería genial! He intentado hacerlo con el siguiente código:

#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 

Pero, cuando trato de ejecutar el código para hacer el espectro de potencia. Específicamente:

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

Me sale este error:

TypeError: mala cantidad de píxeles

Si alguien pudiera señalarme un buen tutorial o ayudarme a editar mi código, sería genial.

Más especificaciones: mis datos están en una matriz de 2d np y puedo cambiar numRows / numCols si es necesario.

Editar:

He resuelto este problema cambiando primero los argumentos de anafast a healpix_map. También mejoré el espacio haciendo mis Nrows * Ncols = 12 * nside * nside. Pero, mi espectro de potencia sigue dando errores. Si alguien tiene enlaces a una buena documentación / tutorial sobre cómo calcular el espectro de potencia (condición de theta / phi args), eso sería increíblemente útil.