Lire un catalogue au format texte et générer des graphiques
Table of Contents
Notions de statistiques pour les traitements des données de catalogues astronomiques.
Point de départ
La documentation en ligne de la bibliothèque Astropy (https://www.astropy.org/) comprend des didacticiels, accessibles à l'adresse http://learn.astropy.org/tutorials.html
Pour explorer ces ressources, j'ai choisi un sujet qui correspond à mes sujets de curiosité du moment : l'exploitation des catalogues d'étoiles.
Les notes qui suivent correspondent à la lecture, l'annotation, l'expérimentation depuis la ressource : http://learn.astropy.org/rst-tutorials/plot-catalog.html?highlight=filtertutorials
Remarque : l'url suivante semble aussi très intéressante :
http://learn.astropy.org/rst-tutorials/conesearch.html?highlight=filtertutorials
En-tête de la ressource
Auteurs
Adrian Price-Whelan, Kelle Cruz, Stephanie T. Douglas
Apprentissages
- Lire un fichier texte avec Astropy
- Convertir les coordonnées en utilisant astropy.coordinates (heures et degrés).
- Réaliser une projection sphérique avec matplotlib.
Sommaire
Ce didacticiel montre comment utiliser astropy.io.ascii pour lire des données contenues dans un fichier texte, astropy.coordinates et astropy.units pour convertir RA (en sexagésimal) en degrés décimaux, et matplotlib pour réaliser un diagramme couleur-magnitude sur une projection de Mollweide.
Début : initialisation des bibliothèques
import numpy as np import matplotlib.pyplot as plt
Une chose n'est pas claire : où est le fichier de données à lire ? J'en ai bricolé un en capturant celui qui s'affiche, mais il est très réduit.
from astropy.io import ascii tbl = ascii.read("simple_table.csv") print(tbl)
name ra dec ------ ---------- --------- BLG100 17:51:00.0 -29:59:48 BLG101 17:53:40.2 -29:49:52 BLG102 17:56:20.2 -29:30:51 BLG103 17:56:20.2 -30:06:22
Les titres de colonnes sont automatiquement affichés au-dessus des colonnes de données.
Maintenant, voyons comment lister la colonne "ra" :
print(tbl["ra"])
ra ---------- 17:51:00.0 17:53:40.2 17:56:20.2 17:56:20.2
L'acension droite, ra, est donnée en coordonnées sexagésimales. Pour effectuer la conversion en degrés et décimales, on va montrer l'utilisation de deux bibliothèques :
import astropy.coordinates as coord import astropy.units as u first_row = tbl[0] ra = coord.Angle(first_row["ra"], unit=u.hour) print(ra.degree)
267.75
267.75
Remarque : les fichiers de données nécessaires sont ici : https://github.com/astropy/astropy-tutorials/tree/master/tutorials/notebooks/plot-catalog
Exploitation d'un catalogue CSV
Rappel : le format CSV est un format texte représentant des données tabulaires (tableau de données) sous forme de valeurs séparées par des virgules (ou des ; ou des tabulations…).
Chargement et exploration
Chargement
tbl = ascii.read("Young-Objects-Compilation.csv") print(tbl.colnames)
['col1', 'col2', 'col3', 'col4', 'col5', 'col6', 'col7', 'col8', 'col9', 'col10', 'col11', 'col12', 'col13', 'col14', 'col15', 'col16', 'col17', 'col18', 'col19', 'col20', 'col21', 'col22', 'col23', 'col24', 'col25', 'col26', 'col27', 'col28', 'col29', 'col30', 'col31', 'col32', 'col33', 'col34']
Commentaires
Les titres de colonnes 'col1', 'col2', 'col3'… ne sont pas les "vrais", mais une substitution qui signifie que la fonction a été incapable de trouver les réponses.
D'ailleurs, si l'on liste…
print(tbl[0])
col1 col2 col3 col4 col5 col6 col7 col8 col9 col10 col11 col12 col13 col14 col15 col16 col17 col18 col19 col20 col21 col22 col23 col24 col25 col26 col27 col28 col29 col30 col31 col32 col33 col34 ---- ---- ---- ---- ---------------- ---- ---- ---- ---- ----- --------------- ----- ----- ----- ----- ----- ----- ----- ------- ----- ----- ----- ---------- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -- -- -- -- 2MASS Photometry -- -- -- -- -- WISE Photometry -- -- -- -- -- -- -- Spectra -- -- -- Astrometry -- -- -- -- -- -- -- -- -- -- --
En fait, le fichier commence par des méta-données qui constituent un en-tête. En supposant que cet en-tête n'occupe qu'une seule ligne, on peut paramétrer la lecture de façon à commencer à la ligne suivante (la numéro 1).
tbl = ascii.read("Young-Objects-Compilation.csv", header_start=1) print(tbl.colnames)
['Name', 'Designation', 'RA', 'Dec', 'Jmag', 'J_unc', 'Hmag', 'H_unc', 'Kmag', 'K_unc', 'W1', 'W1_unc', 'W2', 'W2_unc', 'W3', 'W3_unc', 'W4', 'W4_unc', 'Spectral Type', 'Spectra (FITS)', 'Opt Spec Refs', 'NIR Spec Refs', 'pm_ra (mas)', 'pm_ra_unc', 'pm_dec (mas)', 'pm_dec_unc', 'pi (mas)', 'pi_unc', 'radial velocity (km/s)', 'rv_unc', 'Astrometry Refs', 'Discovery Refs', 'Group/Age', 'Note']
C'est mieux ainsi, mais maintenant que nous avons l'en-tête, il faut aussi charger les données, qui commencent à la ligne suivante.
tbl = ascii.read("Young-Objects-Compilation.csv", header_start=1, data_start=2)
Examen des données
Certaines colonnes ne contiennent pas les données attendues. Par exemple, certaines valeurs de la colonne RA sont manquantes et sont signalées par "-" quand listées.
print(tbl['RA'])
RA
---------
1.01201
6.92489
8.23267
9.42942
11.33929
--
--
--
21.19163
21.5275
...
300.20171
--
303.46467
321.71
--
--
332.05679
333.43715
342.47273
--
350.72079
Length = 64 rows
On appelle cela une colonne masquée (Masked column) parce que certaines valeurs manquantes sont masquées au moment de l'affichage.
On ne peut laisser de vides dans les colonnes si l'on veut exploiter les données, aussi va-t-on leur substituer la valeur conventionnelle nan.
tbl['RA'].filled(np.nan) print(tbl['RA'])
RA
---------
1.01201
6.92489
8.23267
9.42942
11.33929
--
--
--
21.19163
21.5275
...
300.20171
--
303.46467
321.71
--
--
332.05679
333.43715
342.47273
--
350.72079
Length = 64 rows
Ce qui ne ressemble pas aux valeurs attendues.
Sélection des données et représentation graphique simplifiée
Projet : réaliser une représentation graphique magnitude / couleur.
Mise en relation (scatter plot) J-K indice de couleur sur l'axe de X et magnitude J sur l'axe Y. On utilise un "truc" pour retourner l'axe Y : plt.ylim(reversed(plt.ylim())) . Appelé sans paramètre, plt.ylim() renvoie un tuple avec les axes bounds ?? (exemple (0, 10). Si l'on appelle la fonction avec des paramètres, permet de définir les limites des valeurs sur l'axe. Il est alors demandé de mettre les les limites dans l'ordre inverse.
Cela convient bien pour des graphiques réalisés rapidement, mais on a moins de contrôles sur le rendu de la figure.
data = ascii.read("Young-Objects-Compilation.csv", header_start=1, data_start=2) plt.scatter(data["Jmag"] - data["Kmag"], data["Jmag"]) # plot J-K vs. J plt.ylim(reversed(plt.ylim())) # flip the y-axis plt.xlabel("$J-K_s$", fontsize=20) plt.ylabel("$J$", fontsize=20) plt.show()
Projections
Il s'agit ici de réaliser une représentation de la position de certains types d'étoiles sur le ciel… afin de visualiser si leur position est aléatoire ou non.
L'approche est davantage orientée programmation objet :
- Création d'un objet figure.
- Ajout d'un "subplot" à la figure.
- Définition du type de projection utilisé via le paramètre projection
- Dans l'exemple suivant, on utilise la projection de Mollweide.
- Dans celle-ci la convention longitude/RA croît vers la gauche.
L'axe ax, attend les coordonnées angulaires. Mais :
- les valeurs doivent être en radians et l'azimut doit être compris entre -180° et 180°.
- Ceci n'est pas directement paramétrable, aussi faut-il convertir les valeurs.
Astropy contient une classe pour "s'occuper de" des valeurs angulaires astropy.coordinates.Angle
Il est possible de convertir les valeurs de la colonne RA en radians et "wrap" les limites en utilisant cette classe.
ra = coord.Angle(data['RA'].filled(np.nan)*u.degree) ra = ra.wrap_at(180*u.degree) dec = coord.Angle(data['Dec'].filled(np.nan)*u.degree)
Il est alors possible de paramétrer la projection graphique.
fig = plt.figure(figsize=(8,6)) ax = fig.add_subplot(111, projection="mollweide") ax.scatter(ra.radian, dec.radian) plt.show()
Par défaut, matplotlib ajoute des graduations en degrés. Mais il est possible de remplacer par la graduation en heures et d'ajouter une grille de coordonnées.
fig = plt.figure(figsize=(8,6)) ax = fig.add_subplot(111, projection="mollweide") ax.scatter(ra.radian, dec.radian) ax.set_xticklabels(['14h','16h','18h','20h','22h','0h','2h','4h','6h','8h','10h']) ax.grid(True) plt.show()
Sauvergarde de l'image dans un fichier PDF.
fig.savefig("map.pdf")