%matplotlib inline
from __future__ import division

import numpy as np
import healpy as hp
import astroML

#Esto suprime algunos deprecation warnings que molestan pero no interfieren con la ejecución
import warnings
warnings.filterwarnings('ignore')

# Generar un cuadro con versiones de las librerías utilizadas en este notebook
#https://github.com/jrjohansson/version_information
%load_ext version_information
%version_information numpy, healpy, astroML, astroPy

for NSIDE in 2.**np.arange(6):
    print u'El número de pixels para NSIDE = %2d es %5d' %(NSIDE, hp.nside2npix(NSIDE))

NSIDE = 1
ipix = np.arange(hp.nside2npix(NSIDE))
print u'La numeración de los pixels para NSIDE = %d es: \n' %NSIDE, ipix

theta, phi = hp.pix2ang(NSIDE,ipix)
print '\n array theta = ', theta/np.pi, 'radianes'
print '\n array phi = ', phi/np.pi, 'radianes'

NSIDE = 1
hp.ang2pix(NSIDE, np.pi/2, np.pi)

# Caso de arrays de la misma longitud: se emparejan elemento a elemento
theta = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4])
phi = np.array([0, np.pi/2, np.pi, 3*np.pi/2])
hp.ang2pix(NSIDE, theta, phi)

# Caso de arrays de diferente longitud: da un error
theta = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4])
phi2 = np.array([0, np.pi/2, np.pi, 3*np.pi/2, 3.5*np.pi/2])
try:
    hp.ang2pix(NSIDE, theta, phi2)
except Exception as e:
    print u"Atención: ",e

# Pero si por ej. cambiamos el array theta a un array columna
# El broadcasting de numpy es posible y se devuelve un array 2D
# con los números de pixels
theta2 = np.array([0, np.pi/4, np.pi/2, 3*np.pi/4]).reshape(-1,1)
phi2 = np.array([0, np.pi/2, np.pi, 3*np.pi/2, 3.5*np.pi/2])
hp.ang2pix(NSIDE, theta2, phi2)

npixels = hp.nside2npix(32)
np.random.seed(0)
test = np.random.rand(npixels)
hp.mollview(test, min=0, max=1,
            title = u'Prueba de proyección de Mollweide',
            unit=u'Descripción de la unidad de medida')

NSIDE = 2
npixels = hp.nside2npix(NSIDE)
img = np.linspace(0, 255, num=npixels)
index = np.arange(npixels)
theta, phi= hp.pix2ang(NSIDE,index)
hp.mollview(img, min=0, max = 255, unit=u'Rango de valores')
hp.projscatter(theta, phi)
for i in index:
    hp.projtext(theta[i]-0.05, phi[i], i)

lon = [30, 105, 270]
lat = [45, -30, 60]
hp.mollview(title="Mollweide en coordenadas (lon, lat)")
hp.graticule()    # añade retícula
hp.projscatter(lon, lat, lonlat=True)
hp.projtext(30, 45,'(30,45)', lonlat=True)
hp.projtext(105, -30,'(105,-30)', lonlat=True)
hp.projtext(270, 60,'(270,60)', lonlat=True);

from astroML.datasets import fetch_wmap_temperatures
wmap_unmasked = fetch_wmap_temperatures(masked=False)
wmap_unmasked.shape

hp.nside2npix(512)

hp.mollview(wmap_unmasked, min=-1, max=1, 
            title=u'Datos WMAP sin enmascarar \n en coordenadas galácticas',
            unit=r'$\Delta$T (mK)', xsize = 200)

hp.mollview(wmap_unmasked, min=-1, max=1, 
            title=u'Datos WMAP sin enmascarar \n en coordenadas ecuatoriales',
            coord='GC', rot=(180,0),
            unit=r'$\Delta$T (mK)', xsize = 200)
hp.graticule()

hp.projtext(180, 2,'RA=180', lonlat=True)
hp.projtext(90, 2, 'RA=90', lonlat=True)
hp.projtext(270, 2, 'RA=270', lonlat=True);