-->

Etiquetas

Construcción de diagramas de Hertzsprung-Russell en Python

El objeto de este post es construir un diagrama de Hertsprung-Russel y posicionar en el mismo diferentes clases de estrellas

Importaciones y referencias

In [1]:
%matplotlib inline

from __future__ import division

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# 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
Out[1]:
SoftwareVersion
Python2.7.9 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython2.3.1
OSLinux 3.13.0 45 generic x86_64 with debian jessie sid
numpy1.9.1
matplotlib1.4.2
pandas0.15.2
Thu Feb 05 13:11:48 2015 CET

Una interesante actividad educativa acerca de los diagramas de Hertsprung-Russel es The Hertzsprung-Russell Diagram en el SLOAN DIGITAL SKY SURVEY III (SDSS III).

Algo de teoría

Los hechos relevantes para la elaboración de un diagrama de Hertsprung-Russel (diagrama H-R en breve) son los siguientes:

  • Por un lado existe una relación entre la magnitud absoluta de una estrella y su luminosidad, la cual ha sido descrita en este post: El color de los objetos celestes: parte I

  • También existe una relación entre la clase espectral de la estrella (O B A F G K M), su temperatura superficial y su índice de color (por ej. B-V). Esta relación entre índice de color y temperatura se puede consultar también en el mismo post anterior. En cuanto a la relación entre clases espectrales y temperaturas superficiales, puede consultarse en mi post El color de los objetos celestes: parte II

  • El objeto de los diagramas H-R es mostrar la relación existente entre ambos grupos de magnitudes. El diagrama H-R "clásico" muestra la clase espectral de las estrellas (o su correspondiente índice de color B-V) en el eje horizontal y las magnitudes absolutas en el eje vertical.

  • La posición de una estrella en el diagrama H-R nos dice algo acerca de sus carácterísticas. Así por ejemplo si una estrella presenta una clase espectral que indica una temperatura superficial "fria", por ejemplo clases K o M, pero su magnitud absoluta / luminosidad es elevada, esto apunta a que el astro posee una gran superficie, con lo que la elevada superficie radiativa "compensa" la baja intensidad de su radiación debido a su baja temperatura. Es decir, nos hallaríamos ante una estrella "gigante"o "supergigante". Como los índices de color correspondientes a temperaturas frias tienden al rojo, hablaríamos en este caso por ejemplo de gigantes rojas.

El catálogo HIPPARCOS

La mayor dificultad para construir un diagrama H-R es determinar la luminosidad o magnitud absoluta de las estrellas, ya que ello supone conocer sus distancias. En este sentido el catálogo de objetos celestes del SDSS que he venido utilizando hasta ahora en artículos anteriores no es el más adecuado ya que su principal propósito es la fotometría de objetos extragalácticos, disponiendo de poca información de estrellas de nuestra galaxia, y desde luego careciendo de información precisa de la distancia que nos separa de ellas.

Por esta razón, el catálogo de estrellas en el que vamos a basarnos es HIPPARCOS, ya que el objetivo priomordial de esta misión fue precisamente realizar mediciones precisas de cerca de 120000 estrellas. Además, la información en dicho catálogo incluye la clase espectral de la gran mayoría de estos objetos, por lo que es la fuente más adecuada para realizar un diagrama H-R.

Una descripción de los campos en el catálogo puede descargarse en formato pdf en la página siguiente: Hipparcos Overview. En particular, serán relevantes para este artículo los campos siguientes:

  • H1 identificador en el catálogo (HIP)
  • H5 magnitud en el filtro V del sistema Johnson (Vmag)
  • H11 paralaje trigonométrico en milisegundos de arco (Plx)
  • H37 índice de color (B-V)
  • H76 Spectral type (SpType)

Las magnitudes absolutas se calcularán a partir de los valores de V. Para ello se puede partir de la relación entre la magnitud absoluta \(M\), la magnitud aparente \(m\) y la paralaje \(p\) en arcosegundos (ver por ejemplo el artículo Absolute magnitude):

\[ M = m + 5 (1 + \log_{10} p )\]

Cuando la paralaje se expresa en milisegundos de arco (Plx), y tomando como magnitud aparente el valor \(V\) del catálogo, obtenemos:

\[ M = Vmag + 5 (1 + \log_{10} \frac{Plx}{1000}) \]

Que, simplificando queda:

\[ M = Vmag + 5 * log_{10} \frac{Plx}{100} \]

Para descargar el catálogo se debe acceder a la página VizieR Search Page del catálogo I/239 del Centro de Datos Astronómicos de Estrasburgo (CDS). Seleccionar solo los campos indicados anteriormente: HIP, Vmag, Plx, B-V, SpType. En el panel izquierdo, en "Preferences" seleccionar "unlimited" y "; -Separated Values". Pulsar el botón "Submit" y descargar el fichero, que yo he llamado "I_239_selection.tsv" y solo ocupa 3.8 MB. Abriendolo con un editor de textos se comprueba que se trata de un fichero de texto con unas líneas de cabecera explicativas de su contenido, seguidas de 118322 entradas con los cinco campos seleccionados separados por el carácter ;. Al final del fichero hay una línea en blanco que debemos omitir al leer el fichero.

A partir de este momento supondremos pues descargado el fichero "I_239_selection.tsv" en el directorio "../Hipparcos" A continuación se va a leer el fichero creando un dataframe de Pandas:

In [2]:
filename = '../Hipparcos/I_239_selection.tsv'
df = pd.read_table(filename, skiprows=44, sep=';', header=None, index_col=0,
                   names = ['HIP', 'Vmag', 'Plx', 'B-V', 'SpType'],
                   skipfooter=1, engine='python')
In [3]:
df.head()
Out[3]:
Vmag Plx B-V SpType
HIP
1 9.10 3.54 0.482 F5
2 9.27 21.90 0.999 K3V
3 6.61 2.81 -0.019 B9
4 8.06 7.75 0.370 F0V
5 8.55 2.87 0.902 G8III
In [4]:
df.tail()
Out[4]:
Vmag Plx B-V SpType
HIP
118318 6.99 1.92 1.595 K2
118319 8.23 10.63 0.639 G2V
118320 7.59 5.00 0.999 K0
118321 9.20 19.22 0.698 G5V
118322 4.49 8.71 -0.075 B9IV

Preparación del dataframe con los datos

Vamos ahora a trabajar sobre este dataframe. En primer lugar, se analizará si hay entradas en el catálogo sin datos válidos (por ejemplo, con valores NaN)

In [5]:
df.describe()
Out[5]:
Vmag Plx B-V SpType
count 118218 118218 118218 115184
unique 1127 5617 2437 4124
top 8.69 K0
freq 504 263 1281 8570

El valor de "count" indica que de las 118322 entradas disponibles en el catálogo, faltan valores en cada una de las tres columnas, que corresponderán a valores no existentes (NaN). Vamos a identificar y eliminar esas entradas.

Previamente vamos a cambiar cualquier campo con espacios a valores NaN, ya que estos valores no los considera Pandas como NaN, pero no obstante deben ser eliminados ya que con ellos no se podrá realizar cálculos (el método a continuación está tomado de la siguiente respuesta en Stackoverflow

In [6]:
df_clean = df.applymap(lambda x: np.nan if isinstance(x, basestring) and x.isspace() else x)
In [7]:
df_clean.describe()
Out[7]:
Vmag Plx B-V SpType
count 118217 117955 116937 115184
unique 1126 5616 2436 4124
top 8.69 2.93 1.000 K0
freq 504 183 317 8570

Se puede comprobar que los "count" son menores que los anteriores, lo que significa que habia valores, sobre todo en la columna Plx, con espacios, los cuales no fueron considerados como NaN. Ahora están marcados todos como NaN y finalmente, se eliminarán todas las filas que contienen algún valor NaN

In [8]:
df_clean= df_clean.dropna()
df_clean.describe()
Out[8]:
Vmag Plx B-V SpType
count 114472 114472 114472 114472
unique 1072 5361 2426 4070
top 8.69 2.93 1.000 K0
freq 502 182 308 8537
In [9]:
df_clean.shape
Out[9]:
(114472, 4)

Se comprueba que el contador de valores no nulos ahora coincide con el número de líneas del dataframe. A partir de ahora se trabajará con el dataframe df_clean.

A continuación se va a calcular, mediante la fórmula anteriormente indicada, la magnitud absoluta en el filtro V, a la que llamaremos 'M_V'. Se añadirá una columna al dataframe, pero previamente se van a cambiar las columnas del dataframe con las que se va a operar, convirtiendolas de strings a valores numéricos.

In [10]:
df_clean['Vmag'] = df_clean['Vmag'].astype(np.float)
df_clean['Plx'] = df_clean['Plx'].astype(np.float)
df_clean['B-V'] = df_clean['B-V'].astype(np.float)
df_clean['M_V'] = df_clean['Vmag'] + 5 * np.log10(df_clean['Plx']/100.)
df_clean.head()
Out[10]:
Vmag Plx B-V SpType M_V
HIP
1 9.10 3.54 0.482 F5 1.845016
2 9.27 21.90 0.999 K3V 5.972221
3 6.61 2.81 -0.019 B9 -1.146468
4 8.06 7.75 0.370 F0V 2.506509
5 8.55 2.87 0.902 G8III 0.839409

A continuación se van a estudiar los tipos espectrales presentes en el dataframe. Realmente la única información que se necesita de momento de la columna SpType son los dos primeros caracteres. Además solo interesan aquellos casos en los que el primer carácter es alfabético y el segundo numérico, como G8 o F5. Una vez eliminadas las anomalías, se creará una nueva columna con solo los dos primeros caracteres:

In [11]:
# Se eliminan las filas que no cumplen la condición alfa + num
f = lambda s: (len(s) >= 2)  and (s[0].isalpha()) and (s[1].isdigit())
i  = df_clean['SpType'].apply(f)
df_clean = df_clean[i]

# Se crea una columna nueva con los dos primeros caracteres de 'SpType'
f = lambda s: s[0:2]
df_clean['SpType2'] = df_clean['SpType'].apply(f)
In [12]:
df_clean.shape
Out[12]:
(111562, 6)
In [13]:
df_clean.head()
Out[13]:
Vmag Plx B-V SpType M_V SpType2
HIP
1 9.10 3.54 0.482 F5 1.845016 F5
2 9.27 21.90 0.999 K3V 5.972221 K3
3 6.61 2.81 -0.019 B9 -1.146468 B9
4 8.06 7.75 0.370 F0V 2.506509 F0
5 8.55 2.87 0.902 G8III 0.839409 G8

Veamos cuales son las clases espectrales presentes

In [14]:
f = lambda s: s[0]
clases = df_clean['SpType'].map(f)
clases.value_counts()
Out[14]:
K    31578
F    25201
G    22213
A    17651
B    10281
M     4212
O      256
C       82
N       48
R       23
S       17
dtype: int64

Se van a eliminar las líneas con las clases especiales C, N, R y S para quedarnos solo con las de la secuencia OBAFGKM

In [15]:
f = lambda s: s[0] in 'OBAFGKM'
df_clean = df_clean[df_clean['SpType'].map(f)]

Comprobación:

In [16]:
f = lambda s: s[0]
clases = df_clean['SpType'].map(f)
clases.value_counts()
Out[16]:
K    31578
F    25201
G    22213
A    17651
B    10281
M     4212
O      256
dtype: int64

El objetivo que se persigue es realizar un gráfico donde en el eje horizontal figuren las clases espectrales en el sistema de Morgan–Keenan (MKK): OBAFGKM, seguidas de un dígito 0-9. Pero además se pretende que las clases espectrales aparezcan ordenadas precisamente en este orden (B5 antes de A0...). Por ello se van a sustituir las letras por dígitos en secuencia creciente en la columna SpType2 con el fin de poder ordenar por esa columna.

In [17]:
orden = {'O':'0', 'B':'1', 'A':'2', 'F':'3', 'G':'4', 'K':'5', 'M':'6'}
f = lambda s: orden[s[0]]+s[1]
df_clean['SpType2'] = df_clean['SpType2'].apply(f)
df_clean.head()
Out[17]:
Vmag Plx B-V SpType M_V SpType2
HIP
1 9.10 3.54 0.482 F5 1.845016 35
2 9.27 21.90 0.999 K3V 5.972221 53
3 6.61 2.81 -0.019 B9 -1.146468 19
4 8.06 7.75 0.370 F0V 2.506509 30
5 8.55 2.87 0.902 G8III 0.839409 48

Representación gráfica

In [18]:
fig, ax = plt.subplots(figsize=(8,10))

ax.set_xlim(0, 70)
ax.set_ylim(15, -10)
ax.grid()
ax.set_title(u'Diagrama H-R \n (Catálogo Hipparcos)')

ax.title.set_fontsize(20)
ax.set_xlabel('Clase espectral')
ax.xaxis.label.set_fontsize(20)
ax.set_ylabel('Magnitud absoluta')
ax.yaxis.label.set_fontsize(20)

ax.scatter(df_clean['SpType2'].astype(np.int), df_clean['M_V'],
           s=50, edgecolors='none', alpha=0.015, c='k')
ax.set_xticks(range(5,75,10))
ax.set_xticklabels(['O', 'B', 'A', 'F', 'G', 'K', 'M'])
ax.tick_params(axis='both', labelsize=14)

El diagrama anterior no es muy satisfactorio, en el sentido de que, al corresponder el eje horizontal a un conjunto discreto de valores (las diferentes clases espectrales), se forman bandas verticales por acumulación de observaciones en esos valores. Sería preferible utilizar una escala continua de valores decimales. Por este motivo se va a repetir el diagrama utilizando en esta ocasión el índice de color B-V en lugar de las clases espectrales. Además, ello permitirá comprobar que ambas representaciones presentan la misma forma confirmando la relación entre clases espectrales e índices de color.

In [19]:
fig, ax = plt.subplots(figsize=(8,10))

ax.set_xlim(-0.5, 2.5)
ax.set_ylim(15, -10)
ax.grid()
ax.set_title(u'Diagrama H-R \n (Catálogo Hipparcos)')

ax.title.set_fontsize(20)
ax.set_xlabel(u'Índice de color B-V')
ax.xaxis.label.set_fontsize(20)
ax.set_ylabel('Magnitud absoluta')
ax.yaxis.label.set_fontsize(20)

ax.scatter(df_clean['B-V'], df_clean['M_V'],
#           s=50, edgecolors='none', alpha=0.015, c='k')
           s=1, edgecolors='none', c='k')

ax.tick_params(axis='both', labelsize=14)

Clases de luminosidad

Además del tipo espectral OBAFGKM, muchas estrellas en el catálogo Hipparcos contienen, en la codificación del tipo de espectro (SpType), la clase de luminosidad en números romanos: I a VII, que nos indican el tipo de estrella de que se trata. Puede leerse una explicación de este sistema de clasificación en el artículo de wikipedia: Stellar classification.

A continuación se va a analizar la presencia de estas clases de luminosidad en el catálogo:

In [20]:
f = lambda s: 'VII' in s
b = df_clean['SpType'].map(f)
print "Clase VII: enanas blancas, hay %d estrellas" %sum(b)

f = lambda s: ('VI' in s) and ('VII' not in s)
b = df_clean['SpType'].map(f)
print "Clase VI: subenanas, hay %d estrellas" %sum(b)

f = lambda s: ('V' in s) and ('VI' not in s) and ('IV' not in s)
b = df_clean['SpType'].map(f)
print "Clase V: secuencia principal, hay %d estrellas" %sum(b)

f = lambda s: 'IV' in s
b = df_clean['SpType'].map(f)
print "Clase IV: subgigantes, hay %d estrellas" %sum(b)

f = lambda s: 'III' in s
b = df_clean['SpType'].map(f)
print "Clase III: gigantes, hay %d estrellas" %sum(b)

f = lambda s: ('II' in s) and ('III' not in s) and ('VII' not in s)
b = df_clean['SpType'].map(f)
print "Clase II: gigantes brillantes, hay %d estrellas" %sum(b)

f = lambda s: ('I' in s) and ('II' not in s) and ('V' not in s)
b = df_clean['SpType'].map(f)
print "Clase I: supergigantes, hay %d estrellas" %sum(b)
Clase VII: enanas blancas, hay 1 estrellas
Clase VI: subenanas, hay 16 estrellas
Clase V: secuencia principal, hay 24683 estrellas
Clase IV: subgigantes, hay 7955 estrellas
Clase III: gigantes, hay 22519 estrellas
Clase II: gigantes brillantes, hay 1239 estrellas
Clase I: supergigantes, hay 937 estrellas

Como se puede ver a continuación, aproximadamente la mitad de las estrellas no tienen codificada una clase de luminosidad en el campo SpType:

In [21]:
f = lambda s: ('I' not in s) and ('V' not in s)
b = df_clean['SpType'].map(f)
print sum(b)
55605

No obstante, existe un importante número de entradas en el catálogo que sí contienen información de la clase estelar, por lo que a continuación se va a mostrar la posición que ocupa cada una de estas clases de luminosidad en el diagrama H-R

In [22]:
def plot_lum_class(b,c, label):
    ''' b: Series booleana para hacer la selección
        c: Color
        label: etiqueta para la leyenda
    '''
    x = df_clean['B-V'][b]
    y = df_clean['M_V'][b]
    ax.scatter(x, y, c = c, s=6, edgecolors='none', label = label)

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

ax.set_xlim(-0.5, 2.5)
ax.set_ylim(15, -15)
ax.grid()
ax.set_title(u'Diagrama H-R \n (Catálogo Hipparcos)')

ax.title.set_fontsize(20)
ax.set_xlabel(u'Índice de color B-V')
ax.xaxis.label.set_fontsize(20)
ax.set_ylabel('Magnitud absoluta')
ax.yaxis.label.set_fontsize(20)

f = lambda s: 'VII' in s
b = df_clean['SpType'].map(f)
plot_lum_class(b,'white', 'VII: enanas blancas')

f = lambda s: ('VI' in s) and ('VII' not in s)
b = df_clean['SpType'].map(f)
plot_lum_class(b,'blue', 'VI: subenanas')

f = lambda s: ('V' in s) and ('VI' not in s) and ('IV' not in s)
b = df_clean['SpType'].map(f)
plot_lum_class(b,'black', 'V: secuencia principal')

f = lambda s: 'IV' in s
b = df_clean['SpType'].map(f)
plot_lum_class(b,'grey', 'IV: subgigantes')

f = lambda s: 'III' in s
b = df_clean['SpType'].map(f)
plot_lum_class(b,'green', 'III: gigantes')

f = lambda s: ('II' in s) and ('III' not in s) and ('VII' not in s)
b = df_clean['SpType'].map(f)
plot_lum_class(b,'orange', 'II: gigantes brillantes')

f = lambda s: ('I' in s) and ('II' not in s) and ('V' not in s)
b = df_clean['SpType'].map(f)
plot_lum_class(b,'yellow', 'I: supergigantes')

ax.tick_params(axis='both', labelsize=14)
legend = ax.legend(scatterpoints=1,markerscale = 6, shadow=True)
frame = legend.get_frame()
frame.set_facecolor('0.90')

Aquí concluye este artículo. No dudeis en enviarme vuestros comentarios. ¡Hasta la próxima!

3 comentarios:

  1. Excelente trabajo! Muchísimas gracias por compartirlo!

    Tengo interés en hacer algo parecido con 2MASS, pero no me deja bajar el catálogo completo como vos lo hiciste con Hipparcos a través de VizieR. Tenés alguna idea al respecto?

    Saludos,
    Germán.

    ResponderEliminar