El color de los objetos celestes: parte III¶

Autor: Eduardo Martín Calleja

Este post es el tercero que dedico al tema del color (espero que no os hayais aburrido aún del tema), los dos anteriores podeis encontrarlos en:

El color de los objetos celestes: parte I

El color de los objetos celestes: parte II

En este tercer post exploraremos la base de datos del Sloan Digital Sky Survey SDSS, cuyo sistema fotométrico es distinto del UBV, el cual es el que hasta ahora hemos estado utilizando. Echaremos mano de las técnicas para enviar queries SQL a la base de datos del SDSS desde Python que desarrollé en el post Como acceder con SQL a la base de datos del SDSS desde Python. Si vais a seguir este post conviene que le echeis un vistazo antes.

El propósito que perseguimos es acceder a todas las estrellas en la base de datos del SDSS que posean una clasificación espectral en el sistema MK (Morgan-Keenan). La razón de ello es que el catálogo que utiizábamos en el post anterior, el "Bright Star Catalogue" se limitaba a cubrir el conjunto de estrellas que pueden ser visibles a simple vista, por lo cual habia un gran número de estrellas que quedaban fuera del alcance de nuestro estudio, y la proporción de estrellas de cada clase resultaba sesgada. Por otro lado, la base de datos del SDSS es una fuente inagotable de información, por lo que estoy interesado en irme familiarizando con ella.

No obstante, para no hacer muy extensa esta entrada, de momento nos centraremos solo en las estrellas de la secuencia principal.

Esta entrada ha sido escrita íntegramente, como todas las anteriores, con el Notebook de IPython, y en ella iremos viendo con detalle el proceso de obtención, depuración y representación gráfica de los datos.

Importaciones y referencias¶

El SDSS emplea un sistema de cinco filtros cuyas características principales, que he encontrado resumidas en Tesis Doctroral son:

$$ \begin{array}{|l|c|c|} \text{Filtro} & \lambda_{eff} (A) & \Delta \lambda \\\\ \hline \\\\ \text{Ultraviolet (u)} & 3540 & 599 \\\\ \text{Green (g)} & 4770 & 1379 \\\\ \text{Red (r)} & 6222 & 1382 \\\\ \text{Near Infrared (i)} & 7632 & 1535 \\\\ \text{Infrared (z)} & 9049 & 1370 \\\\ \end{array} $$

En el sistema fotométrico del SDSS los filtros en la banda de la luz visible son g y r, por ello se utiliza a menudo la magnitud g-r para indicar el color en el espectro visible.

Estas son las librerías de Python que vamos a utilizar:

In [1]:

%matplotlib inline

from __future__ import division

import numpy as np
import matplotlib.pyplot as plt
import mechanize
import StringIO
import pandas as pd
import quantities as pq

# 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, matplotlib, pandas, quantities

Out[1]:

Software	Version
Python	2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython	2.3.1
OS	Linux 3.13.0 45 generic x86_64 with debian jessie sid
numpy	1.9.1
matplotlib	1.4.2
pandas	0.15.2
quantities	0.10.1
Mon Feb 23 23:34:25 2015 CET

Preparación del acceso a la base de datos del SDSS¶

In [2]:

# URL SQL Search DR10
url = "http://skyserver.sdss3.org/dr10/en/tools/search/sql.aspx"
# Abrimos una instancia de Browser para hacer los queries
br = mechanize.Browser()

Vamos a definir una función para hacer los queries SQL a la base de datos con más comodidad, obteniendo como resultado un dataframe de Pandas:

In [3]:

def SDSS_select(sql):
    br.open(url)
    br.select_form(name="sql")
    br['cmd'] = sql
    br['format']=['csv']
    response = br.submit()
    file_like = StringIO.StringIO(response.get_data())
    return pd.read_csv(file_like,  skiprows=1)

Creación de un dataframe con estrellas en la secuencia principal¶

Dado que queremos obtener datos tanto fotométricos como del espectro (la clasificación espectral del objeto), la tabla más apropiada de la base de datos del SDSS puede ser la vista SpecPhoto, de ella nos van a interesar en particular los datos relativos a las magnitudes en cada filtro, los campos "class" y "type" para seleccionar únicamente estrellas, y el campo "subClass" en el que se encuentra la información relativa a la clasificación espectral del objeto.

Vamos a preparar un query SQL para obtener un dataframe con las estrellas que de acuerdo a los datos de su espectro han sido clasificadas en la base de datos del SDSS como pertenecientes a la secuencia principal.

Hay que tener en cuenta que el campo "type" establece una clasificación de los objetos atendiendo a su morfología: type = 6 significa "STAR". En cambio, "class" clasifica los objetos atendiendo a su espectro en tres categorías: "GALAXY", "QSO", o "STAR". Por tanto, objetos con un "type" igual a 6 en la base de datos pueden ser clasificados espectrográficamente como estrellas, glaxias o quasares. Recíprocamente, aunque la clasificación espectral es más segura, en la base de datos se encuentran también algunos objetos con "class" con el valor "STAR" y con tipos diferentes de 6, como por ejemplo 0 (UNKNOWN) o 3 (GALAXY). Por ello lo mejor para asegurarnos de que estamos seleccionando únicamente estrellas es buscar con type=6 y class='STAR'.

Sabemos por el post anterior, que en el esquema Morgan-Keenan las estrellas de la secuencia principal son las que tienen una clase de luminosidad de tipo "V". Prepararemos el SQL para obtener precisamente las estrellas con esta clasificación:

In [4]:

SQL_V = (
'SELECT '
 'objID, subclass, flags, zwarning, '
 'modelMag_u, modelMag_g, modelMag_r, '
 'modelMag_i, modelMag_z '
'FROM SpecPhoto '
'WHERE '
 'class = "STAR" '
'AND '
 'subclass LIKE "%V%" '
'AND '
 'subclass NOT LIKE "%IV%" '
'AND '
 'subclass NOT LIKE "%VI%" '
'AND '
 'type = 6')

Y lanzamos el query a la base de datos (paciencia, puede tardar unos minutos):

In [5]:

df_V = SDSS_select(SQL_V)
df_V.shape

Out[5]:

(93061, 9)

Cambiamos el nombre de las columnas y echamos un vistazo

In [6]:

df_V.columns = ['objID','subclass','flags','zwarning','u','g','r','i','z']
df_V.head()

Out[6]:

	objID	subclass	flags	u	g	r	i	z
0	1237651752387346781	M2V	35253360001040	25.06914	20.94270	19.25908	18.39424	17.87121
1	1237651752387674306	M2V	68987912448	23.42549	19.97286	18.28300	17.27243	16.69180
2	1237661382768591228	M0V	158398662443776	23.37826	20.52139	18.69852	17.81060	17.31269
3	1237661125617123535	M0V	68987912448	23.77004	21.41212	19.70245	19.00505	18.44237
4	1237661383843053809	M0V	68987912960	22.95323	21.83584	20.00628	18.98125	18.34794

Se ha incluido en la selección, como puede verse, el campo "zwarning" por si interesa filtrar por este campo. En la descripción de la tabla en el Schema Browser se indica al respecto de este campo que un valor igual a cero indica ausencia de problemas. Vamos a ver cuantas entradas hay con este campo diferente de cero:

In [7]:

df_V[df_V['zwarning']!=0].shape

Out[7]:

(13668, 9)

Hay una cantidad respetable. Para estar seguros de que los datos no tienen ningún problema, vamos eliminar de nuestro dataframe todos estos casos.

In [8]:

df_V = df_V[df_V['zwarning']==0]

Y quitamos la columna correspondiente, que ya no tiene utilidad

In [9]:

df_V = df_V.drop('zwarning', axis=1)
df_V.shape

Out[9]:

(79393, 8)

Asignaremos como índice la primera columna:

In [10]:

df_V = df_V.set_index('objID')
df_V.head()

Out[10]:

	subclass	flags	u	g	r	i	z
objID
1237651752387346781	M2V	35253360001040	25.06914	20.94270	19.25908	18.39424	17.87121
1237651752387674306	M2V	68987912448	23.42549	19.97286	18.28300	17.27243	16.69180
1237661382768591228	M0V	158398662443776	23.37826	20.52139	18.69852	17.81060	17.31269
1237661125617123535	M0V	68987912448	23.77004	21.41212	19.70245	19.00505	18.44237
1237661383843053809	M0V	68987912960	22.95323	21.83584	20.00628	18.98125	18.34794

A continuación, vamos a examinar los flags, para ver si tenemos que excluir más entradas del dataframe

Depuración utilizando los flags¶

A título de ejemplo vamos a analizar los flags del primer objeto en el dataframe

In [11]:

flg = df_V.loc[1237662640123281773,'flags']
hex(flg)

Out[11]:

'0x1010000200L'

Vamos a ver su significado. Para ello lo mejor es hacer un query directamente en la base de datos. Ir a la página SQL Search e introducir el siguiente query:

Obtendremos como resultado los nombres de los tres flags que corresponden a los tres bits con valor 1 en el campo "flags" del objeto:

1237662640123281773,STATIONARY BINNED1 MANYPETRO

Vamos ahora a la descripción de los flags: para ello volvamos a la página SQL Search y escribamos:

Como respuesta se obtiene una descripción de cada flag. En particular:

STATIONARY,0x0000001000000000,This object was consistent with being stationary

BINNED1,0x0000000010000000,Object was detected in 1x1 binned image

MANYPETRO,0x0000000000000200,More than one Petrosian radius was found.

Se puede verificar que coinciden los caracteres hexadecimales con el valor obtenido del campo "flags" del primer objeto en el dataframe

Los flags están también explicados en la página Summary table, agrupados por conceptos, lo cual resulta bastante útil. Para interpretar el bit que corresponde a cada flag téngase en cuenta que hay dos campos de 32 bits de flags, que en nuestra tabla SpecPhoto vienen agrupados en un único valor de 64 bits en el campo flags". El campo de 32 bits al que pertenece lo indicaremos con 1/2 o 2/2. Según la página referenciada, los flags que indican posibles problemas con la fotometría son:

CANONICAL_CENTER: bit 0 (1/2)
PEAKCENTER: bit 5 (1/2)
DEBLEND_NOPEAK: bit 14 (2/2)
NOPROFILE: bit 7 (1/2)
NOTCHECKED: bit 19 (1/2)
NOTCHECKED_CENTER: bit 26 (2/2)
TOO_LARGE: bit 24 (1/2)
BADSKY: bit 22 (1/2)

De modo que vamos a construir una máscara para filtrar todos esos casos:

In [12]:

mask = 2**0 + 2**5 + 2**(14+32) + 2**7 + 2**19 + 2**(26+32) + 2**24 + 2**22

Vamos a ver cuantos casos hay en nuestro dataframe que tengan alguno de estos masks activados:

In [13]:

f = lambda s: (s & mask) == 0
b =df_V['flags'].map(f)
b.value_counts()

Out[13]:

True     74733
False     4660
dtype: int64

Es decir, hay 4660 casos con alguno de estos flags activo. Vamos a eliminarlos:

In [14]:

df_V_clean =df_V[b]
df_V_clean.shape

Out[14]:

(74733, 7)

Representación color-color g-r vs. u-g de la secuencia principal¶¶

Veamos en primer lugar cuales son las clases espectrales representadas en el dataframe, para lo cual basta coger el primer carácter del campo "subclass"

In [15]:

f = lambda s: s[0]
clases = df_V_clean['subclass'].map(f)
clases.value_counts()

Out[15]:

F    40908
K    11501
G     8485
A     6184
C     4006
M     2421
B     1228
dtype: int64

Se ve que falta las clases "O" y hay aparentemente una clase que empieza por "C". Veamos que clase espectral es esta:

In [16]:

b = df_V_clean['subclass'].apply(lambda s : s[0] == 'C')
df_V_clean['subclass'][b].value_counts()

Out[16]:

CV    4006
dtype: int64

¡Ah!: Hay 4006 objetos cuyo espectro aparece clasificado como "CV": se trata de "estrellas variables cataclísmicas", las cuales deben ser excluidas de nuestro dataframe que solo queremos que contenga estrellas de la secuancia principal:

In [17]:

df_V_clean = df_V_clean[~b]
df_V_clean.shape

Out[17]:

(70727, 7)

Por otro lado, con el fin de disponer de unos diagramas color-color más completos, vamos a completar el dataframe incorporando estrellas de la clase "O"

In [18]:

# Preparamos la query
#(este sería otro ejemplo de como preparar el string con el SQL)
SQL_O = '                            \
SELECT                               \
 objID, subclass, flags,             \
 modelMag_u, modelMag_g, modelMag_r, \
 modelMag_i, modelMag_z              \
FROM SpecPhoto                       \
WHERE                                \
 class = "STAR"                      \
AND                                  \
 subclass LIKE "O%"                  \
AND                                  \
 type = 6                            \
AND                                  \
 zwarning = 0'

In [19]:

# Y la ejecutamos
df_O = SDSS_select(SQL_O)
df_O.columns = ['objID','subclass','flags','u','g','r','i','z']
df_O = df_O.set_index('objID')
df_O.shape

Out[19]:

(1748, 7)

In [20]:

# depuramos los objetos con flags "malos":
f = lambda s: (s & mask) == 0
b =df_O['flags'].map(f)
df_O_clean =df_O[b]
df_O_clean.shape

Out[20]:

(1694, 7)

Finalmente concatenamos los dos dataframes, creando un nuevo dataframe "df" que es con el que vamos a trabajar:

In [21]:

df = pd.concat([df_V_clean, df_O_clean])
df.shape

Out[21]:

(72421, 7)

Ampliemos el dataframe con nuevas columnas con los índices de color:

In [22]:

df['u-g'] = df['u'] - df['g']
df['g-r'] = df['g'] - df['r']
df['r-i'] = df['r'] - df['i']
df['i-z'] = df['i'] - df['z']

Para distinguir en el diagrama color-color las distintas clases espectrales, vamos a crear un diccionario de colores:

In [23]:

colors = {'O':'#0000FF', 'B':'#CCCCFF', 'A':'#FFFFFF',
          'F':'#FFFFB2', 'G':'#FFFF00', 'K':'#FFA500',
          'M':'#FF0000'}

Y generamos el diagrama color-color:

In [24]:

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, axisbg='0.6')

ax.set_xlim(-1, 6)
ax.set_ylim(-1, 2.5)
ax.grid()
ax.set_title('Secuencia principal: g-r vs. u-g')

ax.title.set_fontsize(20)
ax.set_xlabel('u-g')
ax.xaxis.label.set_fontsize(20)
ax.set_ylabel('g-r')
ax.yaxis.label.set_fontsize(20)

for cls in 'OBAFGKM':
    b = df['subclass'].map(lambda s: s[0] ==cls)
    x = df[b]['u-g']
    y = df[b]['g-r']
    c = colors[cls]
    ax.scatter(x, y, c = c, s=6, edgecolors='none', label = cls)
legend = ax.legend(scatterpoints=1,markerscale = 6, shadow=True)
frame = legend.get_frame()
frame.set_facecolor('0.90')

Como los colores se superponen, no es facil de ver donde se situa cada clase espectral, por lo que vamos a hacer un mosaico de diagramas, uno para cada clase espectral. Empezaremos escribiendo una función para representar gráficamente cada clase:

In [25]:

def plot_gr_ug(df, cls, axes):
    x = df['u-g']
    y = df['g-r']
    axes.scatter(x, y, c = colors[cls], s=6, edgecolors='none')
    axes.set_xlim(-1, 6)
    axes.set_ylim(-1,2.5)
    axes.set_xlabel('u-g')
    axes.set_ylabel('g-r')
    axes.grid()
    axes.set_title('Secuencia principal: class %s' %cls);

Por otro lado, queremos representar también la línea que señala el diagrama color-color de un cuerpo negro, para ver como se posiciona respecto a ella el locus de cada clase. Necesitaremos unas funciones adicionales, ya vistas en anteriores posts:

In [26]:

# Definición de las constantes del sistema fotométrico del SDSS
lambda_u = 3540 * pq.angstrom
delta_u = 599 * pq.angstrom
lambda_g = 4770 * pq.angstrom
delta_g = 1379 * pq.angstrom
lambda_r = 6222 * pq.angstrom
delta_r = 1382 * pq.angstrom
lambda_i = 7632 * pq.angstrom
delta_i = 1535 * pq.angstrom
lambda_z = 9049 * pq.angstrom
delta_u = 1370 * pq.angstrom

In [27]:

def B(wl,T):
    '''wl es un array de longitudes de onda con unidades de longitud
    T es una temperatura expresada en Kelvin
    el resultado es un array de valores de la radiancia espectral
    con unidades W/(m**2 * nm * sr)
    '''
    I = 2 * pq.constants.h * (pq.c)**2 / wl**5 *  \
        1 / (np.exp((pq.constants.h*pq.c \
        / (wl*pq.constants.k*T)).simplified)-1)
    return I.rescale(pq.watt/(pq.m**2 * pq.nm *pq.sr))

Definamos ahora las funciones que nos permitirán calcular los índices de color en el sistema fotométrico del SDSS para un cuerpo negro a una temperatura dada:

In [28]:

# El flujo recibido lo aproximamos como el producto de la radiancia espectral
# B(lambda, T) por la anchura delta lambda

def get_UG(T):    
    F = B(lambda_g, T) * delta_g/(B(lambda_u, T)* delta_u)
    return 2.5 * np.log10(F)

def get_GR(T):    # el flujo recibido lo aproximamos como el producto de la radiancia espectral el lambda por la anchura delta lambda
    F = B(lambda_r, T) * delta_r/(B(lambda_g, T)* delta_g)
    return 2.5 * np.log10(F)

def get_RI(T):
    F = B(lambda_i, T) * delta_i/(B(lambda_r, T)* delta_r)
    return 2.5 * np.log10(F)

Y finalmente vamos a generar un mosaico de diagramas color-color g-r sobre r-i para cada una de las clases espectrales presentes en la secuencia principal, en los cuales vamos a representar también la línea que correspondería a un cuerpo negro perfecto.

In [29]:

fig = plt.figure(figsize=(10,20))
c ='OBAFGKM'
for i in range(0,7):   
    ax1 = fig.add_subplot(4,2,i+1, axisbg='0.4')
    cls = c[i]
    b = df['subclass'].map(lambda s: s[0] ==cls)
    plot_gr_ug(df[b],cls=cls,axes=ax1)
    T1 = 2000 * pq.kelvin
    T2 = 15000 * pq.kelvin
    GmenosR_1 = get_GR(T1)
    UmenosG_1 = get_UG(T1)
    GmenosR_2 = get_GR(T2)
    UmenosG_2 = get_UG(T2)
    # Trazamos la línea "blackbody"
    ax1.plot([UmenosG_1, UmenosG_2], [GmenosR_1, GmenosR_2], c='k',ls='--')

Se observa como las estrellas de la secuencia principal presentan un rango de valores característicos en su índice de color en la región visible, el "g-r". Según vamos avanzando en la secuencia "OBAFGKM", este índice de color va tomando valores mayores, es decir, las estrellas van siendo menos luminosas en esa zona del espectro.

Sin embargo, su índice de color en la región del ultravioleta, "u-g" presenta una dispersión mucho mayor. Basta observar la forma alargada del locus, siendo además más amplia la escala en el eje de abcisas. Realmente se hace dificil distinguir entre las clases F, G, K y M si atendemos solo a este índice.

Otro aspecto a destacar es que las estrellas de clase O, al ser las más calientes, presentan unos valores de los índices de color muy próximos a lo que sería el modelo de radiación de un cuerpo negro.

Representación color-color r-i vs. g-r de la secuencia principal¶¶

Definamos la función para crear el mosaico

In [30]:

def plot_ri_gr(df, cls, axes):
    x = df['g-r']
    y = df['r-i']
    axes.scatter(x, y, c = colors[cls], s=6, edgecolors='none')
    axes.set_xlim(-1, 2.5)
    axes.set_ylim(-1, 2)
    axes.set_xlabel('g-r')
    axes.set_ylabel('r-i')
    axes.grid()
    axes.set_title('Secuencia principal: class %s' %cls);

In [31]:

fig = plt.figure(figsize=(10,20))
for i in range(0,7):
    c ='OBAFGKM'
    ax2 = fig.add_subplot(4,2,i+1, axisbg='0.4')
    cls = c[i]
    b = df['subclass'].map(lambda s: s[0] ==cls)
    plot_ri_gr(df[b],cls=cls,axes=ax2)
    T1 = 2000 * pq.kelvin
    T2 = 15000 * pq.kelvin
    GmenosR_1 = get_GR(T1)
    RmenosI_1 = get_RI(T1)
    GmenosR_2 = get_GR(T2)
    RmenosI_2 = get_RI(T2)

    ax2.plot([GmenosR_1, GmenosR_2], [RmenosI_1, RmenosI_2], c='k', ls='--')

Me resulta interesante observar en este caso como el centro del locus de cada clase espectral se va desplazando a lo largo de la línea del cuerpo negro, desde la zona de temperaturas mas altas a la de temperaturas mas bajas. De nuevo se aprecia un comportamiento especialmente ajustado al espectro del cuerpo negro en las estrellas de la clase O.

Finalizaré aquí esta tercera entrada dedicada al color de las estrellas y los diagramas color-color. Hasta la próxima entrega...