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:
%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
Preparación del acceso a la base de datos del SDSS¶
# 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:
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:
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):
df_V = SDSS_select(SQL_V)
df_V.shape
Cambiamos el nombre de las columnas y echamos un vistazo
df_V.columns = ['objID','subclass','flags','zwarning','u','g','r','i','z']
df_V.head()
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:
df_V[df_V['zwarning']!=0].shape
Hay una cantidad respetable. Para estar seguros de que los datos no tienen ningún problema, vamos eliminar de nuestro dataframe todos estos casos.
df_V = df_V[df_V['zwarning']==0]
Y quitamos la columna correspondiente, que ya no tiene utilidad
df_V = df_V.drop('zwarning', axis=1)
df_V.shape
Asignaremos como Ãndice la primera columna:
df_V = df_V.set_index('objID')
df_V.head()
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
flg = df_V.loc[1237662640123281773,'flags']
hex(flg)
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:
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:
f = lambda s: (s & mask) == 0
b =df_V['flags'].map(f)
b.value_counts()
Es decir, hay 4660 casos con alguno de estos flags activo. Vamos a eliminarlos:
df_V_clean =df_V[b]
df_V_clean.shape
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"
f = lambda s: s[0]
clases = df_V_clean['subclass'].map(f)
clases.value_counts()
Se ve que falta las clases "O" y hay aparentemente una clase que empieza por "C". Veamos que clase espectral es esta:
b = df_V_clean['subclass'].apply(lambda s : s[0] == 'C')
df_V_clean['subclass'][b].value_counts()
¡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:
df_V_clean = df_V_clean[~b]
df_V_clean.shape
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"
# 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'
# 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
# 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
Finalmente concatenamos los dos dataframes, creando un nuevo dataframe "df" que es con el que vamos a trabajar:
df = pd.concat([df_V_clean, df_O_clean])
df.shape
Ampliemos el dataframe con nuevas columnas con los Ãndices de color:
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:
colors = {'O':'#0000FF', 'B':'#CCCCFF', 'A':'#FFFFFF',
'F':'#FFFFB2', 'G':'#FFFF00', 'K':'#FFA500',
'M':'#FF0000'}
Y generamos el diagrama color-color:
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:
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:
# 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
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:
# 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.
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
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);
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...
Hola Eduardo. Veo que usas diferentes fuentes o bases de datos en diferentes posts de tu blog y quiero preguntarte como consigues la información y el acceso a esas bases de datos?
ResponderEliminar