Autor: Eduardo Martín Calleja
En esta entrada, que puede considerarse continuación de la entrada anterior, dedicada a la proyección de Mollweide, vamos a ver una introducción a la utilización del paquete healpy de Python. Healpy proporciona el acceso desde Python al conjunto de funciones de HEALPix, que constituyen el estándar para la representación de los datos procedentes de las distintas misiones que miden la temperatura de la radiación de fondo de microondas del universo.
HEALPIx proporciona un algoritmo de subidivisión de la esfera en elementos de imagen (pixels) con la propiedad de que todos ellos representan exactamente el mismo área de la superficie esférica original, y además estos pixeles se distribuyen en líneas de igual latitud. A continuación esta subdivisión se puede trasladar a una proyección en el plano como la de Mollweide que conserva las áreas en la representación plana resultante.
Disponer de representaciones gráficas de la superficie esférica con pixels de igual área es crucial para la presentación y posterior análisis de la distribución de masas, energía, radiación, etc., permitiendo comparar densidades relativas y aplicar diversos algoritmos. En particular, la librería HEALPix se desarrolló inicialmente para representar la distribución de la radiación de fondo de microondas cuyos datos se han ido recopilando en diversas misiones, si bien su utilización se ha extendido posteriormente a otros campos.
Importaciones y referencias¶
%matplotlib inline
from __future__ import division
import numpy as np
import healpy as hp
import astroML
#Esto suprime algunos deprecation warnings que molestan
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
Uso de coordenadas y numeración de pixels¶
En healpy se utilizan coordenadas esféricas para posicionar los pixeles: el ángulo polar "theta" \(\theta\) medido desde el polo Norte (colatitud), el cual toma valores en el intervalo \([0,\pi]\), y el ángulo azimutal "phi" \(\varphi\) (longitud) que toma valores en el intervalo \([0,2 \pi]\), medido hacia el Este del globo.
La resolución de la subdivisión en pixels se define con el parámetro NSIDE.
La función nside2npix(NSIDE) devuelve el número de pixels en que se descompondrá la supericie esférica para un NSIDE dado. La resolución mínima son 12 pixels, y cada vez que se dobla NSIDE (que debe ser siempre una potencia de 2) se divide cada pixel en 4, como se ve en la siguiente tabla:
for NSIDE in 2.**np.arange(6):
print u'El número de pixels para NSIDE = %2d es %5d' %(NSIDE, hp.nside2npix(NSIDE))
Los pixels se numeran por defecto según el esquema RING, que consiste en irlos numerando de forma consecutiva a lo largo de anillos de latitud constante, comenzando por el por el anillo de pixels situado más al Norte.
La función pix2ang(nside,ipix) devuelve las coordenadas angulares para un número de pixel o array de números de pixels dado. Téngase en cuenta que theta está comprendido entre 0 y \(\pi\) radianes, y phi entre 0 y \(2\pi\) radianes. Veamos un ejemplo de utilización de esta función:
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'
Inversamente, La función ang2pix(nside, theta, phi) devuelve la numeración de los pixels para unas coordenadas angulares dadas \((\theta, \varphi)\), que al igual que antes pueden ser coordenadas individuales o sendos arrays
NSIDE = 1
hp.ang2pix(NSIDE, np.pi/2, np.pi)
Veamos que pasa ahora si utilizamos sendos arrays con las coordenadas angulares:
# 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)
Representación gráfica en proyección de Mollweide¶
La función para representar gráficamente un array de pixels en la proyección de Mollweide es mollview(). Para ser representado con mollview() el array debe ser de una dimensión, y su longitud debe coincidir con nside2npix(NSIDE). Cada elemento del array corresponderá a un pixel en la descomposición de HEALPix. Veamos un ejemplo:
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')
Podemos ver el esquema de numeración de los pixels en el ejemplo siguiente, lo que de paso nos permitirá introducir las funciones projscatter() y projtext(). Se utilizará el esquema por defecto: "RING":
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)
Observese que, en el esquema anterior, el orden de representación de los pixels es de arriba (Norte) abajo (Sur), y de derecha (Oeste) a izquierda (Este), como corresponde a una proyección de la bóveda celeste, en que el Este se situa a la izquierda. Si se desea representar un mapa geográfico, con el Este a la derecha, se deberá utilizar la opción flip="geo"
healpy trabaja internamente con coordendas angulares theta y phi en radianes, pero se le pueden proporcionar también coordenadas en otros sistemas, por ejemplo coordenadas galácticas (lon, lat) o coordenadas ecuatoriales (RA, Dec). Por este motivo, las funciones projscatter() y projtext() admiten la opción lonlat = True, que indica que las coordenadas que se proporcionan no son theta (colatitud) y phi (longitud) en radianes, sino una longitud y una latitud en grados. Por ejemplo, en coordenadas ecuatoriales las coordenadas angulares serían una ascensión recta y una declinación.
Veamos un ejemplo siguiendo este post
Recuerdese que el Este se representa a la izquierda de la figura
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);
Descarga y visualización de los datos de la radiación de fondo de microondas (CMB) de la misión WMAP¶
Para acceder a los datos del fóndo cósmico de microondas del WMAP utilizaremos el módulo astroML de Python. La función para obtener un array listo para su representación gráfica con healpy es:
astroML.datasets.fetch_wmap_temperatures()
Se pueden descargar dos versiones de estos datos: una versión "completa" que incluye las frecuencias contaminantes, principalmente radiación proveniente de nuestra propia galaxia u otras fuentes extragalácticas, o bien una versión en la que las bandas de frecuencias donde se concentra esta contaminación han sido enmascaradas.
Descarguemos por ejemplo los datos no enmascarados:
from astroML.datasets import fetch_wmap_temperatures
wmap_unmasked = fetch_wmap_temperatures(masked=False)
wmap_unmasked.shape
Estos datos están en forma de un array de Numpy ya preparado para su representación gráfica con healpy. Puede comprobarse que el nivel de resolución gráfica es el que corresponde a un NSIDE = 512
hp.nside2npix(512)
Procedamos ahora a su representación gráfica:
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)
La ordenación de los pixels en el array anterior corresponde al sistema de coordenadas galácticas, como podemos verificar por la presencia en el gráfico, en la posición lat = 0, de una banda de contaminación térmica procedente de nuestra galaxia. Esto es irrelevante en lo que respecta a la construcción por healpy de la representación gráfica, ya que se limita a leer el array y colocar cada pixel en la posición que le corresponda en el esquema de numeración adoptado (RING por defecto) y según el nivel de resolución indicado por la longitud del array.
No obstante, podemos indicar a healpy que estamos trabajando en un sistema de coordenadas determinado y deseamos hacer un cambio a otro sistema de coordenadas diferente. healpy puede trabajar con tres tipos de coordenadas: "G": Galácticas, "E": Eclipticas y "C": Celestes Ecuatoriales. Por ejemplo, si deseamos redibujar el gráfico convirtiendo las coordenadas galácticas en ecuatoriales utilizaremos la opción: coord='GC'.
Más aún, podemos incluso girar el gráfico anterior para por ejemplo centrar el origen en RA, dec = (180,0), obteniendo así un mapa en el que se aprecia la situación de la "vía láctea" en coordenadas ecuatoriales centradas en RA = 180 grados
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);
Este comentario ha sido eliminado por el autor.
ResponderEliminarEste comentario ha sido eliminado por el autor.
Eliminar