Introducción a Jupyter Notebooks para Data Preview 0.2
Autoría: Melissa Graham
Última verificación de ejecución: 2023-11-29
Versión de las Pipelines Científicas de LSST: Weekly 2023_47
Tamaño del contenedor (Container size): medium
Nivel de aprendizaje: principiante
Descripción: una introducción al uso de los Jupyter Notebooks y los paquetes de Python de Rubin para acceder a los productos de datos de LSST (imágenes y catálogos).
Habilidades: ejecutar código Python en un Jupyter Notebook. Usar el servicio TAP para obtener datos del catálogo Object. Utilizar Butler para obtener y visualizar una imagen deepCoadd.
Productos de Datos LSST: TAP dp02_dc2_catalogs. Tabla Object. Imagen deepCoadd con Butler.
Paquetes: lsst.rsp.get_tap_service, lsst.rsp.retrieve_query, lsst.daf.butler, lsst.afw.display, lsst.geom, pandas, matplotlib
Créditos: desarrollado originalmente por Melissa Graham y el equipo científico de la comunidad de Rubin (Rubin Community Science Team) en el contexto de la Vista Previa de Datos de Rubin DP0.1.
Soporte: se pueden encontrar recursos y documentación relacionada a DP0 en dp0-2.lsst.io. Las preguntas son bienvenidas como nuevos temas en la categoría Support - Data Preview 0 del foro de la comunidad de Rubin. El equipo de Rubin responderá a todas las preguntas publicadas allí.
1.0. Introducción¶
Este Jupyter Notebook da una introducción a cómo funcionan los notebooks. Muestra cómo ejecutar celdas de código y de texto markdown, cómo importar paquetes de Python y aprender acerca de sus módulos, y provee enlaces a documentacion adicional.
Este Notebook también muestra la funcionalidad básica de la Plataforma Científica de Rubin (RSP - Rubin Science Platform) instalada en el Centro Interino de Datos (IDF - Interim Data Facility; Google Cloud) incluyendo cómo usar el servicio TAP para consultar y obtener datos de catálogos; matplotlib para graficar datos de catálogos; el paquete Butler de LSST para consultar y obtener datos de imágenes; y el paquete afwDisplay de LSST para visualizar imágenes.
Este notebook utiliza el conjunto de datos Vista Previa de Datos 0.2 (DP0.2 - Data Preview 0.2). Este conjunto de datos utiliza un subconjunto de imágenes simuladas del Desafío de Datos 2 de DESC (DC2 - DESC's Data Challenge 2), que han sido reprocesadas por el Observatorio Rubin usando la versión 23 de las Pipelines Científicas de LSST (LSST Science Pipelines). Se puede encontrar más información sobre los datos simulados en el artículo científico de DC2 de DESC y en la documentación de la divulgación de datos de DP0.2.
1.1. Cómo usar un Jupyter Notebook¶
Los Jupyter Notebooks contienen una combinación de código, salida, visualizaciones y texto narrativo. La fuente de documentación más completa acerca de los Jupyter Notebooks es https://jupyter-notebook.readthedocs.io, pero hay numerosos y muy buenos tutoriales y demostraciones para principiantes disponibles. Generalmente una búsqueda en la web de una pregunta, como por ejemplo "cómo hacer una tabla markdown en jupyter notebook", llevará a varios buenos ejemplos. A menudo las respuestas se hallarán en StackOverflow.
Un Jupyter Notebook es una serie de celdas (cells). Hay tres tipos de celdas: código, markdown y raw (en bruto). Este texto se generó a partir de una celda markdown. Arriba en la barra de menú vas a encontrar un menú desplegable para establecer el tipo de celda.
Advertencia: todas las celdas de código en un notebook deben ser ejecutadas en el orden en que aparecen.
Hacer clic en la siguiente celda de código: con el cursor en la celda, presionar simultaneamente "shift" (tecla Mayús) y "enter" (tecla Entrar) para ejecutar el código de la celda.
# Esta es una celda de código. Presionar shift-enter para ejecutar.
# El # convierte estas líneas en comentarios, no código. No se ejecutan.
print('Hello, world!')
Hello, world!
Hacer doble clic sobre ESTAS PALABRAS EN ESTA CELDA MARKDOWN para ver el código fuente markdown.
1.1.1. Guías prácticas para Jupyter Notebooks¶
Cómo ejecutar rápidamente todas las celdas: ir a la barra de menú de más arriba y seleccionar "Kernel", luego "Restart Kernel and Run All Cells" (Reiniciar Kernel y Ejecutar Todas las Celdas).
Cómo detener de emergencia un notebook: si una celda de código está demorando mucho tiempo en ejecutarse (e.g., si se inició por accidente un proceso para obtener los datos de un catálogo entero) mátalo yendo a "Kernel" en el menú de más arriba y seleccionando "Restart Kernel and Clear All Outputs" (Reiniciar el Kernel y Limpiar Todas las Salidas). Puede igualmente tomar algunas decenas de segundos, pero detendrá el proceso y reiniciará el kernel.
El kernel es el motor de cómputo para el notebook (la RSP usa un kernel de python3)
y se puede pensar como un compilador en vivo.
Reiniciar el kernel y limpiar todas las salidas significa que todas las variables o funciones definidas son eliminadas de la memoria,
y todas las celdas de código vuelven al estado "sin-ejecutar".
Cómo ver una tabla de contenidos para este notebook: hacer clic en el ícono de la lista con viñetas en la barra del menú vertical a la izquierda, y aparecerá a la izquierda una tabla de contenidos (ToC - Table of Contents) generada automáticamente. Hacer clic en el ícono de la carpeta de archivos que está en lo más alto de la barra de menú vertical a la izquierda, para volver a la vista de directorio.
Cómo saber qué versión de las Pipelines Científicas de LSST se está ejecutando: mirar a lo largo de la barra inferior del navegador, y buscar la versión de las Pipelines Científicas de LSST que se seleccionó como "imagen". Es probablemente "Recommended (Weekly aaaa_ss)" -donde "aaaa" se refiere al año y "ss" a la semana- y debería corresponderse con la versión verificada que se lista en la cabecera del notebook. Alternativamente, descomentar de la siguiente celda de código las dos líneas y ejecutarla.
# ! echo $IMAGE_DESCRIPTION
# ! eups list -s | grep lsst_distrib
1.2. Importación de paquetes¶
La mayoría de los Jupyter Notebooks comienzan cargando / importando todos los paquetes que necesitarán en la primera celda de código.
No es necesario conocer en profundidad esos paquetes para completar este tutorial, pero aquí hay un poco de información básica y algunos links para seguir aprendiendo.
numpy: un paquete fundamental para cómputo científico con arreglos (arrays) en Python. Su documentación completa está disponible en numpy.org, e incluye guías rápidas para principiantes. (El paquete numpy no se utiliza en este notebook, pero se importa como ejemplo porque es un paquete de uso muy frecuente).
matplotlib: este paquete es una biblioteca muy completa para crear visualizaciones estáticas, animadas e interactivas en Python. Toda su documentación se encuentra en matplotlib.org. La galería de matplotlib es un gran punto de partida y tiene enlaces a ejemplos.
pandas: un paquete que permite trabajar de manera eficiente con datos tabulares en dataframes (marcos de datos). Se puede aprender más en la documentación de Pandas.
astropy: un paquete de Python con herramientas de astronomía muy útiles. Se puede aprender más en la documentación de astropy.
lsst: estos paquetes son todos de las Pipelines Científicas de LSST.
El paquete lsst.rsp habilita el acceso a imágenes y catálogos a través del servicio TAP (ver Sección 2);
el paquete lsst.daf.butler habilita el acceso a imágenes y catálogos a través de Butler (ver Sección 3);
el paquete lsst.geom posee funciones auxiliares para metadatos de imágenes y el paquete lsst.afw.display habilita la visualización de las imágenes (ver Sección 3).
Importar los paquetes utilizados en este notebook ejecutando la siguiente celda.
import numpy
import matplotlib
import matplotlib.pyplot as plt
import pandas
from lsst.rsp import get_tap_service, retrieve_query
import lsst.daf.butler as dafButler
import lsst.geom
import lsst.afw.display as afwDisplay
1.2.1. Obtener información sobre los paquetes de Python importados¶
Mostrar las versiones de numpy y matplotlib.
print('numpy version: ', numpy.__version__)
print('matplotlib version: ', matplotlib.__version__)
numpy version: 1.24.4 matplotlib version: 3.8.2
Se puede ver una lista emergente de cualquier módulo de un paquete escribiendo el nombre del paquete, luego un punto y finalmente presionando la tecla Tab. Usar las flechas hacia arriba y abajo para desplazarse a través de la lista emergente. Esto funciona tanto si la línea de código está comentada como si no lo está. En la celda de abajo, numpy. está comentado porque no es una sentencia de código que se pueda ejecutar y, si no estuviera el #, esta celda fallaría en la ejecución (pruébalo -- elimina el #, presiona shift-enter y mira cómo falla).
# numpy.
Utilizar la funcion "help" para ver la documentación de ayuda para el paquete. Borrar el símbolo # para descomentar cualquiera de las líneas y ejecutar la siguiente celda. La documentación de ayuda puede ser demasiado larga. Re-comentar la línea reintroduciendo el #, luego re-ejecutar la celda y la salida desaparecerá.
# help(numpy)
# help(matplotlib)
# help(numpy.abs)
# help(matplotlib.pyplot)
2.0. Datos de Catálogos¶
2.1. Servicio de Protocolo de Acceso a Datos Tabulados (TAP)¶
El protocolo de acceso a datos tabulados TAP (Table Access Protocol) provee un acceso estandarizado a los datos de catálogos para exploración, búsqueda y acceso. La documentación completa de TAP es provista por la Alianza Internacional de Observatorios Virtuales (IVOA --International Virtual Observatory Alliance).
El servicio TAP utiliza un lenguaje de consultas similar a SQL (Structured Query Language --lenguaje de consulta estructurado) denominado ADQL (Astronomical Data Query Language --lenguaje de consulta de datos astronómicos). La documentación para ADQL incluye más información sobre sintaxis y palabras clave.
Aviso: no todas las funcionalidades de ADQL están soportadas por RSP para la Vista Previa de Datos 0 (DP0).
Iniciar el servicio TAP.
service = get_tap_service("tap")
2.2. Explorando tablas y columnas de catálogos con TAP¶
Este ejemplo utiliza el catálogo Object de objetos de DP0.2, que contiene fuentes detectadas en imágenes coagregadas (coadded) (también llamadas imágenes apiladas --stacked--, combinadas --combined--, o deepCoad).
Los contenidos del catálogo también se pueden explorar con el navegador de esquemas de DP0.2 .
Los resultados de una búsqueda usando el servicio TAP se visualizan de mejor manera usando una de estas dos funciones:
.to_table(): convierte los resultados en una tabla de astropy.
.to_table().to_pandas(): convierte a una tabla de astropy y luego a un dataframe de Pandas.
Advertencia: no usar el método .to_table().show_in_notebook(). Esto puede provocar problemas en el entorno de Jupyterlab de RSP que podría bloquear tu notebook indefinidamente.
Los tres ejercicios opcionales de abajo enseñan diferentes maneras de explorar usando el servicio TAP. Muestran cómo usar el servicio TAP con sentencias ADQL para explorar qué catálogos existen y qué columnas contienen. Cada celda utiliza un método diferente para visualizar los resultados de la búsqueda con TAP. Eliminar todos los #, ejecutar cada celda y observar que se crean muchas salidas -- agregar el # nuevamente a cada línea, re-ejecutar la celda y la salida desaparecerá.
2.2.1. Ejercicio 1¶
Obtener y mostrar una lista de todos los nombres y descripciones de las tablas que están disponibles a través del servidor TAP.
# my_adql_query = "SELECT description, table_name FROM TAP_SCHEMA.tables"
# results = service.search(my_adql_query)
# results_table = results.to_table().to_pandas()
# results_table
2.2.2. Ejercicio 2¶
Obtener y mostrar una lista de los nombres de los campos (nombres de columna) en el esquema TAP del catálogo de Objetos de DP0.2. Notar que los resultados se pueden nombrar de cualquier otra manera; aquí, se usa alternativamente 'res'.
# my_adql_query = "SELECT * from TAP_SCHEMA.columns "+\
# "WHERE table_name = 'dp02_dc2_catalogs.Object'"
# res = service.search(my_adql_query)
# print(res.fieldnames)
2.2.3. Ejercicio 3¶
Obtener los nombres, tipos de datos, descripciones y unidades para todas las columnas en el catálogo Object. Mostrar el número de columnas.
# my_adql_query = "SELECT column_name, datatype, description, unit "+\
# "FROM TAP_SCHEMA.columns "+\
# "WHERE table_name = 'dp02_dc2_catalogs.Object'"
# results = service.search(my_adql_query)
# results_table = results.to_table().to_pandas()
# print('Number of columns available in the Object catalog: ', len(results_table))
Mostrar todos los 991 nombres de columnas y su información. ¡Son tantas salidas! Comentar la línea en la celda y re-ejecutar la celda para hacer que todas las salidas desaparezcan.
# results_table
Mostrar sólo nombres y descripciones para columnas que contengan la cadena (texto) "cModelFlux". Probar otras cadenas como "coord", "extendedness", "deblend", o "detect".
# my_string = 'cModelFlux'
# for col,des in zip(results_table['column_name'],results_table['description']):
# if col.find(my_string) > -1:
# print('%-40s %-200s' % (col,des))
2.3. Obteniendo Datos con TAP¶
Algunos consejos sobre cómo hacer consultas eficientes en los catálogos de DP0.2.
Restricciones en RA, Dec producen consultas más rápidas: los Servicios de Consulta de LSST (Qserv - Query Services) proveen acceso a la base de datos de catálogos del LSST. Los catálogos se pueden consultar utilizando el lenguaje estándar de consultas SQL con unas pocas limitaciones. Qserv almacena los datos de catálogos fragmentados por coordenadas (RA, Dec). Las sentencias de consultas ADQL que incluyen restricciones por coordenadas no requieren una búsqueda en todo el catálogo, y son generalmente más rápidas (y pueden ser mucho más rápidas) que las sentencias de consulta ADQL que sólo incluyen restricciones para otras columnas.
Obtener una pequeña muestra de filas:
como se mostró en la Sección 2.3.2, usar maxrec=10 o SELECT TOP 10 cuando se explora el conjunto de datos para obtener como resultado unas pocas filas con las que jugar (esto también puede acortar los tiempos de consulta para consultas exploratorias sin las sentencias WHERE).
Reestricción recomendada en detect_isPrimary:
cuando corresponda, se recomienda incluir detect_isPrimary = True en las consultas para los catálogos Object, Source y ForcedSource.
Este parámetro toma el valor True si una fuente
no tiene hijos,
si está en la zona interior de una parcela coagregada (coadd patch),
si está en la zona interior de una región coagregada (coadd tract),
y si no se detecta en un pseudo-filtro.
Al incluir esta restricción se eliminan duplicados (i.e., no se incluirán al mismo tiempo a padres e hijos obtenidos en la separación --deblend).
2.3.1. Convirtiendo flujos a magnitudes¶
Los catálogos de objetos y fuentes almacenan sólo flujos. Hay cientos de columnas relacionadas a flujos y almacenarlas también como magnitudes sería redundante y un desperdicio de espacio.
Todas las unidades de flujo son nanojanskys ($nJy$). La página de Wikipedia sobre Magnitudes AB proporciona un recurso conciso para quienes no estén familiarizados con el uso de magnitudes AB y flujos jansky. Para convertir $nJy$ a magnitudes AB usar: $m_{AB} = -2.5log( f_{nJy}) + 31.4$.
Como se muestra en la Sección 2.3.2, para agregar columnas de magnitudes después de obtener columnas de flujos, se puede hacer lo siguiente:
results_table['r_calibMag'] = -2.50 * numpy.log10(results_table['r_calibFlux']) + 31.4
results_table['r_cModelMag'] = -2.50 * numpy.log10(results_table['r_cModelFlux']) + 31.4
Como se muestra en la Sección 2.3.3, para obtener columnas de flujos como magnitudes en una consulta ADQL, se puede hacer lo siguiente:
scisql_nanojanskyToAbMag(g_calibFlux) as g_calibMag,
y las columnas de errores de magnitudes se pueden obtener con:
scisql_nanojanskyToAbMagSigma(g_calibFlux, g_calibFluxErr) as g_calibMagErr.
2.3.2. Diez objetos de cualquier tipo¶
Para mostrar rápidamente cómo obtener datos del catálogo Object, usar una búsqueda de cono y solicitar que se devuelvan sólo 10 registros. La Figura 2 del artículo de DC2 de DESC muestra que la región del cielo cubierta por la simulación de DC2 contiene las coordenadas RA,Dec = 62,-37.
Este ejemplo usa maxrec=10 en la función service.search(), pero se pueden lograr los mismos resultados reemplazando SELECT por SELECT TOP 10 en la consulta ADQL.
Atención: el catálogo Object contiene cientos de millones de filas. Las búsquedas que no especifican una región y/o un número máximo de registros pueden tomar mucho tiempo y devolver demasiadas filas para mostrar en un notebook.
Obtener coordenadas y magnitudes g,r,i para diez objetos dentro de un radio de 0.5 grados de 62,-37.
use_center_coords = "62, -37"
Abajo, SELECT TOP 10 se utiliza en la sentencia de consulta para limitar a 10 objetos los datos que se devuelven.
Una alternativa es usar la palabra clave maxrec en la sentencia de búsqueda: service.search(my_adql_query, maxrec=10).
Sin embargo, usar maxrec puede devolver un DALOverflowWarning que alerta que se han devuelto resultados parciales (aún cuando se deseaban resultados parciales).
my_adql_query = "SELECT TOP 10 "+ \
"coord_ra, coord_dec, detect_isPrimary, " + \
"r_calibFlux, r_cModelFlux, r_extendedness " + \
"FROM dp02_dc2_catalogs.Object " + \
"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), " + \
"CIRCLE('ICRS', " + use_center_coords + ", 0.01)) = 1 "
results = service.search(my_adql_query)
results_table = results.to_table()
results_table['r_calibMag'] = -2.50 * numpy.log10(results_table['r_calibFlux']) + 31.4
results_table['r_cModelMag'] = -2.50 * numpy.log10(results_table['r_cModelFlux']) + 31.4
results_table
| coord_ra | coord_dec | detect_isPrimary | r_calibFlux | r_cModelFlux | r_extendedness | r_calibMag | r_cModelMag |
|---|---|---|---|---|---|---|---|
| deg | deg | nJy | nJy | ||||
| float64 | float64 | bool | float64 | float64 | float64 | float64 | float64 |
| 62.0095695 | -37.0030527 | False | 115.5597619 | 107.20676 | 1.0 | 26.242983404124434 | 26.32444457256579 |
| 61.999653 | -37.0037443 | False | 142.1429817 | 76.2996347 | 0.0 | 26.018186446761295 | 26.69369385328298 |
| 62.0024481 | -37.0066927 | False | -- | -- | -- | -- | -- |
| 61.9954064 | -37.0080442 | False | 1062.1604372 | 1092.7958693 | 1.0 | 23.834524697476997 | 23.803652388163506 |
| 61.9977835 | -37.0087985 | False | 261.1418936 | 197.5926925 | -- | 25.357808627606243 | 25.660572802022923 |
| 61.9961705 | -37.0056238 | False | 117.6636971 | 48.4746206 | 1.0 | 26.223393774519916 | 27.18621395340402 |
| 61.9977822 | -37.0095762 | False | 94.7492789 | 42.5903896 | -- | 26.45856021648528 | 27.326720967885816 |
| 61.9956804 | -37.0035825 | False | 46.7946246 | 32.0731837 | -- | 27.224510081023972 | 27.63464482073427 |
| 61.9958404 | -37.0015949 | False | 21.184015 | 39.0455005 | -- | 28.08497931181146 | 27.421072514620995 |
| 61.9962257 | -37.0006286 | False | 152.1184444 | 87.7436118 | 1.0 | 25.944545311095364 | 26.541961231628306 |
2.3.3 Diez mil objetos puntuales¶
Además de la búsqueda de cono, fijar restricciones de consulta para que detect_isPrimary tenga valor True (esto no devolverá fuentes "hijas" producto de la separación --deblending), que el flujo calibrado sea mayor que 360 nJy (aproximadamente magnitud 25) y que los parámetros de extensión sean 0 (fuentes puntuales).
Obtener las magnitudes en las bandas g, r e i para 10000 objetos que tienen alta probabilidad de ser estrellas.
results = service.search("SELECT TOP 10000 "
"coord_ra, coord_dec, "
"scisql_nanojanskyToAbMag(g_calibFlux) as g_calibMag, "
"scisql_nanojanskyToAbMag(r_calibFlux) as r_calibMag, "
"scisql_nanojanskyToAbMag(i_calibFlux) as i_calibMag, "
"scisql_nanojanskyToAbMagSigma(g_calibFlux, g_calibFluxErr) as g_calibMagErr "
"FROM dp02_dc2_catalogs.Object "
"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), "
"CIRCLE('ICRS', "+use_center_coords+", 1.0)) = 1 "
"AND detect_isPrimary = 1 "
"AND g_calibFlux > 360 "
"AND r_calibFlux > 360 "
"AND i_calibFlux > 360 "
"AND g_extendedness = 0 "
"AND r_extendedness = 0 "
"AND i_extendedness = 0")
results_table = results.to_table()
print(len(results_table))
10000
La muestra de la tabla se truncará automáticamente.
results_table
| coord_ra | coord_dec | g_calibMag | r_calibMag | i_calibMag | g_calibMagErr |
|---|---|---|---|---|---|
| deg | deg | ||||
| float64 | float64 | float64 | float64 | float64 | float64 |
| 61.6570015 | -37.2804304 | 22.031272582886093 | 20.841869048065924 | 19.633224165064846 | 0.006548812053841391 |
| 61.5654545 | -37.2963864 | 22.568669591960138 | 21.94191211446436 | 21.719320973459546 | 0.010549103249775003 |
| 61.5639412 | -37.2881681 | 21.0869429923351 | 20.71128125317726 | 20.569107250519423 | 0.0029110825474043095 |
| 61.6580862 | -37.2905435 | 23.35878140951924 | 22.14970701036416 | 21.063171986383644 | 0.021889802393035156 |
| 61.5719555 | -37.2879985 | 24.99380135837644 | 24.035047121790495 | 22.914051359802265 | 0.09619117973374319 |
| 61.6925881 | -37.2897205 | 24.01602993552452 | 21.964467256686607 | 21.386827715941486 | 0.04074438061067066 |
| 61.5770609 | -37.3281371 | 18.083387757958448 | 17.74096309651904 | 17.616044503268515 | 0.00029055416952623517 |
| 61.6620278 | -37.3342887 | 20.06331506003999 | 19.03255728561365 | 18.59664835524638 | 0.0012304012370138574 |
| 61.6613321 | -37.3334399 | 19.21919376721202 | 18.973003454214414 | 18.900900425071242 | 0.0006474951150186445 |
| ... | ... | ... | ... | ... | ... |
| 62.8928622 | -37.3063833 | 20.781802163607647 | 20.39203435662265 | 20.246455248021046 | 0.0021238197541374104 |
| 62.7489924 | -37.2732423 | 22.12803987442131 | 21.00940296650269 | 20.534428444366867 | 0.0064227266307295364 |
| 62.8035395 | -37.3383226 | 20.640204463047148 | 20.146494823846375 | 19.962718467569054 | 0.0017968678550262553 |
| 62.8585036 | -37.3399787 | 19.80909406850452 | 19.601135027482712 | 19.547192908791843 | 0.0009610062487019912 |
| 62.8325929 | -37.3069738 | 18.560153891879253 | 17.711876921428697 | 17.392439121029557 | 0.0004006822165554178 |
| 62.5229987 | -37.5871432 | 21.530962950197896 | 20.177988086812682 | 18.44079036859312 | 0.004129621099510703 |
| 62.515065 | -37.5538895 | 21.64928125141401 | 21.295857019853116 | 21.181519169469237 | 0.004389355944178713 |
| 62.6426043 | -37.5515414 | 23.6765541695956 | 22.89866911057431 | 22.555546085717182 | 0.029090963314695923 |
| 62.6446212 | -37.5506902 | 18.818481568148286 | 17.673830648251833 | 17.163218434572794 | 0.0004918488295918679 |
| 62.6905634 | -37.5568301 | 23.791633100616124 | 22.522001918588344 | 20.802947135673712 | 0.03425834109730646 |
Guardar los resultados en un dataframe de pandas para tener acceso fácil al contenido. Estos datos son utilizados para crear un diagrama color-magnitud en la Sección 2.4.
data = results_table.to_pandas()
Para quienes no estén familiarizados con Pandas, aquí hay algunas líneas de código opcionales que muestran cómo imprimir los nombres de las columnas, la información de la columna 'ra' o los valores de la columna 'ra'. Descomentar (eliminar #) y ejecutar la celda para ver la salida de la demostración.
# data.columns
# data['coord_ra']
# data['coord_ra'].values
2.4 Hacer un diagrama color-magnitud¶
Usar la función plot para graficar del paquete matplotlib.pyplot (que fue importado como plt). Los parámetros de la función plot se describen en su totalidad en este sitio web de matplotlib, pero en resumen son: valores de x, valores de y, forma del símbolo ('o' es círculo), tamaño de marcador (marker size, ms) y transparencia del marcador (alpha).
plt.plot(data['r_calibMag'].values - data['i_calibMag'].values,
data['g_calibMag'].values, 'o', ms=2, alpha=0.2)
plt.xlabel('mag_r - mag_i', fontsize=16)
plt.ylabel('mag_g', fontsize=16)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.xlim([-0.5, 2.0])
plt.ylim([25.5, 16.5])
plt.show()
Este gráfico genera muchos interrogantes, como por ejemplo "¿Por qué los colores están concentrados discretamente?" y "¿Son todas esas realmente estrellas?". Las respuestas exceden el ámbito de este notebook y son dejadas como temas potenciales de análisis científico que pueden abordarse con el conjunto de datos DC2.
2.5 Opcional: graficar magnitud versus error de magnitud¶
Para ilustrar tanto las magnitudes como los errores de las magnitudes obtenidos mediante la consulta TAP de más arriba, aquí hay una opción para graficar el error de la magnitud como función de la magnitud.
# plt.plot(data['g_calibMag'].values, data['g_calibMagErr'].values, 'o', ms=2, alpha=0.2)
# plt.show()
3.0. Datos de Imágenes¶
Los dos tipos de imágenes más comunes con los que van a interactuar quienes tengan acceso a DP0 son calexps y deepCoadds.
calexp: una única imagen en un único filtro.
deepCoadd: una combinación de imágenes individuales apiladas en profundidad o coagregadas.
Las Pipelines Científicas de LSST (Science Pipelines) procesan y almacenan imágenes en regiones (tracts) y parcelas (patches).
región (tract): una porción del cielo dentro de la teselación del cielo completo (mapa del cielo) de LSST (LSST all-sky tessellation); dividido en parcelas.
parcela (patch): una subregión cuadrilátera de una región, de un tamaño que puede almacenarse fácilmente en la memoria de una computadora de escritorio.
Para obtener y visualizar una imagen en las cordenadas deseadas, se debe especificar el tipo de imagen, la región (tract) y la parcela (patch).
3.1. Crear una instancia de Butler¶
Butler --que en inglés significa mayordomo-- (documentación) es un paquete de software de las Pipelines Científicas de LSST que sirve para obtener datos de LSST sin necesidad de conocer su ubicación o formato. Además Butler también puede ser utilizado para explorar y descubrir qué datos existen. Otros tutoriales muestran la funcionalidad completa de Butler.
Crear una instancia de Butler utilizando la siguiente configuración y colección de DP0.2. Devolverá información sobre las credenciales encontradas.
butler = dafButler.Butler('dp02', collections='2.2i/runs/DP0.2')
3.2. Identificar y obtener una imagen deepCoadd¶
Hay un cúmulo de galaxias DC2 que se ve genial en RA = 03h42m59.0s, Dec = -32d16m09s (en grados, 55.745834, -32.269167).
Usar lsst.geom para definir un SpherePoint --punto en la esfera-- para las coordenadas del cúmulo (documentación de lsst.geom).
my_ra_deg = 55.745834
my_dec_deg = -32.269167
my_spherePoint = lsst.geom.SpherePoint(my_ra_deg*lsst.geom.degrees,
my_dec_deg*lsst.geom.degrees)
print(my_spherePoint)
(55.7458340000, -32.2691670000)
Obtener el mapa del cielo de DC2 desde Butler y usarlo para identificar la región y parcela para las coordenadas del cúmulo (documentación de skymap).
skymap = butler.get('skyMap')
tract = skymap.findTract(my_spherePoint)
patch = tract.findPatch(my_spherePoint)
my_tract = tract.tract_id
my_patch = patch.getSequentialIndex()
print('my_tract: ', my_tract)
print('my_patch: ', my_patch)
my_tract: 4431 my_patch: 17
Utilizar Butler para obtener la imagen deepCoadd en la banda i para la región y parcela.
dataId = {'band': 'i', 'tract': my_tract, 'patch': my_patch}
my_deepCoadd = butler.get('deepCoadd', dataId=dataId)
3.3. Visualizar la imagen con afwDisplay¶
Los datos de imágenes obtenidos con Butler se pueden visualizar de diferentes maneras. Una opción sencilla es utilizar el paquete afwDisplay de las Pipelines Científicas de LSST. Hay cierta documentación de afwDisplay disponible y otros tutoriales de DP0 entran en más detalle acerca de todas las opciones de visualización (e.g., superponer datos de máscaras para mostrar píxeles malos).
Establecer el backend de afwDisplay para matplotlib.
afwDisplay.setDefaultBackend('matplotlib')
Usar afwDisplay para mostrar los datos obtenidos de la imagen.
La siguiente celda de código crea una figura matplotlib.pyplot; define una instancia de lsst.afw.display.Display como afw_display; establece la escala para las intensidades o tonos de los píxeles (pixel shading); visualiza la imagen usando mtv; y activa las etiquetas para los ejes x, y (coordenadas del píxel).
fig = plt.figure(figsize=(10, 8))
afw_display = afwDisplay.Display(1)
afw_display.scale('asinh', 'zscale')
afw_display.mtv(my_deepCoadd.image)
plt.gca().axis('on')
(11899.5, 16099.5, 7899.5, 12099.5)
Para conocer más sobre el paquete afwDisplay y sus funcionalidades, usar la función de ayuda help.
# help(afw_display.scale)
# help(afw_display.mtv)