01. Introduccion a DP03#

DP03_01_Introduccion_a_DP03_ES

Introducción a DP0.3
Autoría: Bob Abel, Douglas Tucker y 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 vista general del contenido de los catálogos de objetos en movimiento de DP0.3.

Habilidades: utilizar el servicio TAP y ADQL para acceder a las tablas de DP0.3.

Productos de Datos LSST: TAP dp03_catalogs.

Paquetes: lsst.rsp.get_tap_service

Créditos: el conjunto de datos DP0.3 fue generado por miembros de los equipos de Pipelines del Sistema Solar y Puesta a Punto de Rubin --Rubin Solar System Pipelines and Commissioning teams-- con ayuda del grupo de la Colaboración Científica del Sistema Solar de LSST --LSST Solar System Science Collaboration--, en particular: Pedro Bernardinelli, Jake Kurlander, Joachim Moeyens, Samuel Cornwall, Ari Heinze, Steph Merritt, Lynne Jones, Siegfried Eggl, Meg Schwamb y Mario Jurić.

Mario Jurić brindó asistencia esencial en la etapa inicial de la creación de este notebook. Quienes construyeron este notebook además reconocen la ayuda de Leanne Guy, Pedro Bernardinelli, Sarah Greenstreet, Megan Schwamb, Brian Rogers, Niall McElroy y Jake Vanderplass, entre otras personas.

Soporte: se pueden encontrar recursos y documentación relacionados con DP0.3 en dp0-3.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. Introducción¶

Este notebook muestra cómo acceder al conjunto de datos simulado de Vista Previa de Datos 0.3 (DP0.3 --Data Preview 0.3) en la Plataforma Científica de Rubin (RSP - Rubin Science Platform).

Para la simulación DP0.3, sólo se simularon objetos en movimiento y sólo fueron creados los cátalogos (no hay imágenes). La simulación DP0.3 es enteramente independiente y separada de la simulación DP0.2.

DP0.3 es un catálogo híbrido que contiene objetos del Sistema Solar tanto reales como simulados (asteroides, objetos próximos a la Tierra, Troyanos, objetos transneptunianos, cometas e incluso una nave espacial simulada... pero ningún planeta principal ni la Luna). Ver la documentación de DP0.3 para más información acerca de cómo fue creado el catálogo híbrido.

En las Operaciones de Rubin --Rubin Operations--, estas tablas estarían cambiando constantemente, siendo actualizadas todos los días con los resultados de las observaciones de la noche anterior. Sin embargo, para DP0.3, se han simulado tanto un catálogo estático de 1 año como un catálogo estático de 10 años.

Aviso: para reiterar, en DP0.3, hay tablas disponibles tanto para catálogos de 1 año como para catálogos de 10 años. Sin embargo, durante el resto de este notebook, a menos que se indique lo contrario, sólo consideraremos las tablas del catálogo de 10 años. Se recomienda explorar las tablas de los catálogos de 1 año por cuenta propia.

1.1. Importación de paquetes¶

Importar paquetes generales de python y las utilidades del servicio TAP de Rubin.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colormaps
import pandas as pd
from lsst.rsp import get_tap_service

2. Crear el cliente del servicio SSO TAP de Rubin¶

El conjunto de datos DP0.3 está disponible mediante el servicio de protocolo de acceso a datos tabulados (TAP --Table Access Protocol) y puede ser accedido mediante sentencias ADQL (Astronomical Data Query Language --Lenguaje de consulta de datos astronómicos).

TAP 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). ADQL es similar a SQL (Structured Query Language --lenguaje de consulta estructurado). La documentación de ADQL incluye más información sobre sintaxis y palabras clave.

Obtener una instancia del servicio SSO TAP, y corroborar que éste exista.

Aviso: el servicio TAP de DP0.3 se llama ssotap (mientras que tap se utiliza para DP0.2).

In [2]:

service = get_tap_service("ssotap")
assert service is not None

2.1. ¿Qué esquemas de `ssotap` están disponibles?¶

Crear una consulta ADQL para seleccionar todos (*) los esquemas disponibles, y usar el servicio TAP para ejecutar la consulta.

Usar to_table para convertir a una tabla astropy, la cual visualizará los resultados de manera amigable.

In [3]:

query = "SELECT * FROM tap_schema.schemas"
results = service.search(query).to_table()
results

Out[3]:

Table length=4

description	schema_index	schema_name	utype
str512	int32	str64	str512
Data Preview 0.3 (ten-year version): Contains the catalog products of a Solar System Science Collaboration simulation of the results of SSO analysis of the wide-fast-deep data from the full LSST ten-year dataset.	0	dp03_catalogs_10yr
Data Preview 0.3 (one-year version): Contains the catalog products of a Solar System Science Collaboration simulation of the results of SSO analysis of the wide-fast-deep data, based on analyzing only the first year of LSST observations.	1	dp03_catalogs_1yr
A TAP-standard-mandated schema to describe tablesets in a TAP 1.1 service	100000	tap_schema
UWS Metadata	120000	uws

2.2 ¿Qué tablas de DP0.3 están disponibles?¶

Hay cuatro tablas para la simulación de 1 año y la de 10 años: MPCORB, SSObject, SSSource y DiaSource.

Se podrán encontrar las descripciones de las tablas y sus esquemas, más información y consejos sobre acceso y consulta de las tablas de DP0.3 (incluyendo qué columnas están actualmente vacías), en la documentación de las definiciones de los productos de datos de DP0.3.

Este tutorial explorará las tablas individualmente, y la combinación de tablas será mostrada en un tutorial futuro (ver también este consejo sobre combinaciones de tablas en DP0.3).

Seleccionar todas las tablas disponibles y mostrar información acerca de ellas.

In [4]:

query = "SELECT * FROM tap_schema.tables " \
        "WHERE tap_schema.tables.schema_name LIKE 'dp03_catalogs%' " \
        "ORDER BY table_index ASC"
results = service.search(query).to_table()
results

Out[4]:

Table length=8

description	schema_name	table_index	table_name	table_type	utype
str512	str512	int32	str64	str8	str512
LSST-computed per-object quantities. 1:1 relationship with MPCORB. Recomputed daily, upon MPCORB ingestion.	dp03_catalogs_1yr	1	dp03_catalogs_1yr.SSObject	table
LSST-computed per-object quantities. 1:1 relationship with MPCORB. Recomputed daily, upon MPCORB ingestion.	dp03_catalogs_10yr	1	dp03_catalogs_10yr.SSObject	table
The orbit catalog produced by the Minor Planet Center. Ingested daily. O(10M) rows by survey end. The columns are described at https://minorplanetcenter.net//iau/info/MPOrbitFormat.html	dp03_catalogs_1yr	2	dp03_catalogs_1yr.MPCORB	table
The orbit catalog produced by the Minor Planet Center. Ingested daily. O(10M) rows by survey end. The columns are described at https://minorplanetcenter.net//iau/info/MPOrbitFormat.html	dp03_catalogs_10yr	2	dp03_catalogs_10yr.MPCORB	table
Table to store 'difference image sources'; - sources detected at SNR >=5 on difference images.	dp03_catalogs_1yr	3	dp03_catalogs_1yr.DiaSource	table
Table to store 'difference image sources' - sources detected at SNR >=5 on difference images.	dp03_catalogs_10yr	3	dp03_catalogs_10yr.DiaSource	table
LSST-computed per-source quantities. 1:1 relationship with DIASource. Recomputed daily, upon MPCORB ingestion.	dp03_catalogs_1yr	4	dp03_catalogs_1yr.SSSource	table
LSST-computed per-source quantities. 1:1 relationship with DiaSource. Recomputed daily, upon MPCORB ingestion.	dp03_catalogs_10yr	4	dp03_catalogs_10yr.SSSource	table

3. La tabla `MPCORB`¶

Durante las Operaciones de Rubin, el Procesamiento del Sistema Solar ocurrirá durante el día, después de una noche de observación. Asociará las detecciones de objetos en movimiento en imágenes diferenciales y reportará los descubrimientos al Centro de Planetas Menores (MPC --Minor Planet Center), como así también computará las propiedades derivadas (magnitudes, ajustes de curvas de fase, coordenadas en varios sistemas).

El MPC calculará los parámetros orbitales y esos resultados se enviarán de vuelta a Rubin, y serán almacenados y puestos a disposición a la comunidad usuaria como la tabla MPCORB (las otras propiedades derivadas son almacenadas en las otras tres tablas exploradas abajo). La tabla MPCORB de DP0.3 es una simulación de cómo será este producto de datos después de 10 años de LSST.

Aviso: el MPC contiene todos los objetos en movimiento reportados en el Sistema Solar, y no se limita a aquellos detectados por LSST. De esta manera, la tabla MPCORB tendrá más filas que la tabla SSObject.

Aviso: para DP0.3, el MPC no realizó ajustes y la tabla MPCORB contiene los elementos orbitales usados en la simulación (la tabla MPCORB es una tabla con valores de verdad).

Para más información acerca de los planes de Rubin para el Procesamiento del Sistema Solar, ver la Sección 3.2.2 del Documento de Definición de Productos de Datos (Data Products Definitions Document, DPDD). Es importante tener en cuenta que todavía existen diferencias entre la Tabla 4 del DPDD, que contiene el esquema anticipado para las tablas de objetos en movimiento, y los esquemas de tablas de DP0.3.

3.1. Tamaño¶

Usar el servicio TAP para contar el número de filas en la tabla MPCORB.

In [5]:

results = service.search("SELECT COUNT(*) FROM dp03_catalogs_10yr.MPCORB")
results.to_table()

Out[5]:

Table length=1

COUNT1
int64
14462388

Hay 14.5 millones de filas en la tabla MPCORB.

3.2. Columnas¶

Usar el servicio TAP para consultar sobre información de las columnas de MPCORB.

Imprimir los resultados como una tabla de pandas.

In [6]:

results = service.search("SELECT column_name, datatype, description, "
                         "       unit from TAP_SCHEMA.columns "
                         "WHERE table_name = 'dp03_catalogs_10yr.MPCORB'")
results.to_table().to_pandas()

Out[6]:

	column_name	datatype	description	unit
0	arc	float	MPCORB: Arc (days), for single-opposition objects	d
1	arcEnd	char	MPCORB: Year of last observation (for multi-op...
2	arcStart	char	MPCORB: Year of first observation (for multi-o...
3	computer	char	MPCORB: Computer name
4	e	double	MPCORB: Orbital eccentricity
5	epoch	double	MPCORB: Epoch (in MJD, .0 TT)	d
6	flags	int	MPCORB: 4-hexdigit flags. See https://minorpla...
7	fullDesignation	char	MPCORB: Readable designation
8	incl	double	MPCORB: Inclination to the ecliptic, J2000.0 (...	deg
9	lastIncludedObservation	float	MPCORB: Date of last observation included in o...	d
10	mpcDesignation	char	MPCORB: Number or provisional designation (in ...
11	mpcG	float	MPCORB: Slope parameter, G
12	mpcH	float	MPCORB: Absolute magnitude, H	mag
13	mpcNumber	int	MPC number (if the asteroid has been numbered;...
14	n	double	MPCORB: Mean daily motion (degrees per day)	deg/d
15	nobs	int	MPCORB: Number of observations
16	node	double	MPCORB: Longitude of the ascending node, J2000...	deg
17	nopp	int	MPCORB: Number of oppositions
18	peri	double	MPCORB: Argument of perihelion, J2000.0 (degrees)	deg
19	pertsLong	char	MPCORB: Precise indicator of perturbers (blank...
20	pertsShort	char	MPCORB: Coarse indicator of perturbers (blank ...
21	q	double	MPCORB: Perihelion distance (AU)	AU
22	reference	char	MPCORB: Reference
23	rms	float	MPCORB: r.m.s residual (")	arcsec
24	ssObjectId	long	LSST unique identifier (if observed by LSST)
25	tperi	double	MPCORB: MJD of pericentric passage	d
26	uncertaintyParameter	char	MPCORB: Uncertainty parameter, U

En algunos casos, las descripciones de las columnas se cortan en la tabla anterior.

Ejecutar lo siguiente para ver, por ejemplo, la descripción completa para la columna mpcDesignation.

In [7]:

results['description'][10]

Out[7]:

'MPCORB: Number or provisional designation (in packed form)'

3.3. Obtener un subconjunto al azar¶

Para obtener un subconjunto de filas al azar, hacer uso del hecho de que ssObjectId es un entero sin signo de tipo long unsigned integer de 64-bits asignado aleatoriamente. Ya que ADQL interpreta al tipo long unsigned integer de 64-bits como un entero con signo de 63-bits, estos varían desde un valor entero negativo muy grande hasta un valor entero positivo muy grande. Esto será arreglado en el futuro para que todos los identificadores sean números positivos.

Aviso: por el uso de ssObjectId, la siguiente consulta devuelve un subconjunto de filas al azar de MPCORB que están asociadas con una fila en la tabla SSObject. En otras palabras, esto limita la consulta a sólo obtener objetos en movimiento desde la tabla MPCORB que han sido detectados por LSST.

Primero, averiguar el rango completo de los valores de ssObjectId usando las funciones MIN y MAX de ADQL.

In [8]:

results = service.search("SELECT min(ssObjectId), max(ssObjectId)"
                         "FROM dp03_catalogs_10yr.MPCORB")
results.to_table()

Out[8]:

Table length=1

min1	max2
int64	int64
-9223370383071521539	9223370875126069107

Definir un rango de búsqueda para ssObjectId que devuelva no más del 1% de todos los objetos en MPCORB. Esto lo haremos estimando un nuevo mínimo ssObjectId que esté 1% por debajo del máximo ssObjectId para todo el rango de los valores de ssObjectId.

Aviso: ya que el rango de ssObjectId (-9223370383071521539 --> +9223370875126069107) es mucho más grande que el número de filas en la tabla MPCORB (14600302), no esperamos obtener con este método exactamente el 1% de las filas de MPCORB, pero deberíamos obtener aproximadamente el 1%, siempre y cuando los valores de ssObjectId estén distribuidos de manera razonablemente uniforme sobre su extenso rango.

In [9]:

min_val = int(results[0].get('min1'))
max_val = int(results[0].get('max2'))
print('Rango completo: ', min_val, max_val)

min_val = int(max_val - 0.01*(max_val-min_val))
print('Rango de búsqueda: ', min_val, max_val)

Rango completo:  -9223370383071521539 9223370875126069107
Rango de búsqueda:  9038903462544093184 9223370875126069107

Ejecutar la busqueda y obtener todas (*) las columnas de la tabla MPCORB.

Guardar los resultados en df como un marco de datos de pandas (dataframe).

In [10]:

results = service.search("SELECT * FROM dp03_catalogs_10yr.MPCORB "
                         "WHERE ssObjectId BETWEEN 9038903462544093184 "
                         "AND 9223370875126069107")
df = results.to_table().to_pandas()

Mostrar el dataframe resultante. Notar que es truncado automáticamente tanto en las columnas como en las filas (busca los "..." a mitad de camino hacia abajo y a mitad de camino a través de la tabla).

In [11]:

df

Out[11]:

	arc	arcEnd	arcStart	computer	e	epoch	flags	fullDesignation	incl	lastIncludedObservation	...	nopp	peri	pertsLong	pertsShort	q	reference	rms	ssObjectId	tperi	uncertaintyParameter
0	NaN				0.143340	60065.000000	0	2011 1062 T-3	6.230810	NaN	...	0	62.518599			1.975493		NaN	9154661686854389525	59969.017393
1	NaN				0.112557	60065.000000	0	2011 1093 T-2	6.322560	NaN	...	0	72.646946			2.061410		NaN	9059993893438746957	60310.434546
2	NaN				0.049899	60065.000000	0	2011 1184 T-2	4.306140	NaN	...	0	279.168196			2.426167		NaN	9170955388068648464	60191.630055
3	NaN				0.105047	60065.000000	0	2011 1186 T-3	16.188500	NaN	...	0	357.091672			2.757674		NaN	9077066921741880727	59093.291192
4	NaN				0.176216	60065.000000	0	2011 1464 T-2	4.495110	NaN	...	0	81.310304			2.121439		NaN	9053987333942294876	59808.773603
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
144467	NaN				0.919897	60255.126890	0	2011 LPCD750	106.184487	NaN	...	0	209.420137			6.642560		NaN	9133287149408159006	60251.935331
144468	NaN				0.527282	61444.848873	0	2011 LPCD823	11.480147	NaN	...	0	92.507010			16.379375		NaN	9140884185118476131	61444.842271
144469	NaN				0.536121	62656.471653	0	2011 LPCD889	117.420822	NaN	...	0	12.259819			15.988392		NaN	9189745973987900787	62656.437497
144470	NaN				0.602921	63093.436806	0	2011 LPCD940	44.898648	NaN	...	0	81.527019			15.572242		NaN	9181407122537090723	63092.423677
144471	NaN				0.655530	62223.336444	0	2011 LPCD941	177.067536	NaN	...	0	276.729975			12.396072		NaN	9211044985939642883	62223.322988

144472 rows × 27 columns

Como se puede ver, se obtuvieron 144,472 filas, las cuales -- como se esperaba -- son casi (pero no exactamente) el 1% de las 14,462,388 filas en la tabla MPCORB.

Aviso: hay numerosas columnas que actualmente contienen valores NaN (not a number -- no es un número). Para los datos simulados de DP0.3 estas columnas podrían ser reemplazadas en el futuro cercano y para las divulgaciones de datos reales no habrá ninguna columna con todos sus valores en NaN. Si se desea, se pueden quitar todas las columnas con NaN haciendo, e.g. df.dropna(axis=1, how='all', inplace=True). De todas formas, la mejor práctica es entender las columnas y obtener sólo lo que se vaya a utilizar.

Opcional: usar la función info del dataframe de pandas para aprender más acerca de los valores en las columnas obtenidas. Descomentar la siguiente celda de código (eliminar el símbolo #) y ejecutar la celda para ver la información sobre df.

In [12]:

# df.info()

Opcional: usar la función describe de pandas para mostrar estadísticas de las columnas numéricas en el dataframe (conteo, media, desviación estándar, etc.).

In [13]:

# df.describe()

3.4. Graficar histogramas de columnas seleccionadas¶

Wikipedia provee una guía adecuada de elementos orbitales para principiantes.

Para una referencia rápida, abajo se muestran distribuciones para cinco elementos orbitales clave y la magnitud absoluta $H$ (ver la Sección 4.4 para una descripción de $H$).

In [14]:

fig, ax = plt.subplots(2, 3, figsize=(10, 6), sharey=False)
ax[0, 0].hist(df['e'], bins=100, log=True)
ax[0, 0].set_xlabel('Excentricidad')
ax[0, 0].set_ylabel('log(Número)')
ax[0, 1].hist(df['incl'], bins=100, log=True)
ax[0, 1].set_xlabel('Inclinación (deg)')
ax[0, 1].set_ylabel('log(Número)')
ax[0, 2].hist(df['mpcH'], bins=100, log=True)
ax[0, 2].set_xlabel('Magnitud Absoluta, H (mag)')
ax[0, 2].set_ylabel('log(Número)')
ax[1, 0].hist(df['node'], bins=50)
ax[1, 0].set_xlabel('Longitud del Nodo Ascendente (deg)')
ax[1, 0].set_ylabel('Número')
ax[1, 0].set_ylim(0,3500)
ax[1, 1].hist(df['peri'], bins=50)
ax[1, 1].set_xlabel('Argumento del Perihelio (deg)')
ax[1, 1].set_ylabel('Número')
ax[1, 1].set_ylim(0,3500)
ax[1, 2].hist(df['q'], bins=100, log=True)
ax[1, 2].set_xlabel('Distancia al Perihelio (au)')
ax[1, 2].set_ylabel('log(Número)')
fig.suptitle('Histogramas de Elementos Orbitales Clave')
fig.tight_layout()
plt.show()

4. La tabla `SSObject`¶

Durante las Operaciones de Rubin, el Procesamiento de Datos Inmediatos (Prompt Processing) se realizará durante la noche, detectando fuentes en imágenes diferenciales (a las cuales llamaremos DIASources --fuentes resultantes del Análisis de Imágenes Diferenciales o fuentes DIA--, ver Sección 6, tabla DiaSource) y asociándolas a transitorios y variables estáticos en el cielo (llamados DIAObjects de la tabla DIAObject, no incluidos en DP0.3).

El Procesamiento del Sistema Solar que ocurre durante el día, luego de una noche de observación, vincula las detecciones de objetos en movimiento DiaSource con objetos SSObject (Objeto del Sistema Solar). Mientras que la tabla MPCORB contiene los elementos orbitales para esos objetos en movimiento, los objetos de SSObject contienen las propiedades medidas por Rubin como ajustes de curvas de fase y magnitudes absolutas.

Aviso: no se inyectaron artefactos ni fuentes de imágenes diferenciales espurias en los catálogos de DP0.3.

Aviso: aunque hay columnas para ellas, no se simularon datos para las bandas u o Y en DP0.3. A menos que se indique lo contrario, ignoraremos los datos de las bandas u o Y para el resto de este notebook.

4.1. Tamaño¶

Usar la función de conteo de ADQL para devolver el tamaño.

In [15]:

results = service.search("SELECT COUNT(*) from dp03_catalogs_10yr.SSObject")
results.to_table().to_pandas()

Out[15]:

	COUNT1
0	4443479

El conjunto de datos DP0.3 contiene 4.4 millones de objetos del sistema solar detectados por Rubin.

Esto es menos que los 14.5 millones de objetos en el catálogo MPCORB. Se deja como ejercicio en la Sección 7 determinar las características de esos objetos de la tabla MPCORB que faltan en la tabla SSObject.

4.2. Columnas¶

Opcional: mostrar una lista de nombres de columnas, tipos de datos, descripciones y unidades.

In [16]:

# results = service.search("SELECT column_name, datatype, description, "
#                          "unit from TAP_SCHEMA.columns "
#                          "WHERE table_name = 'dp03_catalogs_10yr.SSObject'")
# results.to_table().to_pandas()

4.3. Obtener un subconjunto al azar¶

Usar básicamente la misma consulta que se utilizo para la tabla MPCORB, arriba.

In [17]:

results = service.search("SELECT * FROM dp03_catalogs_10yr.ssObject "
                         "WHERE ssObjectId BETWEEN 9038903462544093184 "
                         "AND 9223370875126069107")
df = results.to_table().to_pandas()

In [18]:

df

Out[18]:

	arc	discoverySubmissionDate	firstObservationDate	flags	g_Chi2	g_G12	g_G12Err	g_H	g_H_gG12_Cov	g_HErr	...	y_H_yG12_Cov	y_HErr	y_Ndata	z_Chi2	z_G12	z_G12Err	z_H	z_H_zG12_Cov	z_HErr	z_Ndata
0	2225.923096	61634.00000	61640.42378	102	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	0	NaN	NaN	NaN	NaN	NaN	NaN	0
1	3148.112305	60308.00000	60311.20320	66	1.821222	0.538449	0.029092	16.201378	0.000087	0.003903	...	NaN	NaN	0	1.004799	0.570081	2.714864e-02	15.410084	6.308406e-05	4.074336e-03	68
2	3314.104736	60544.00000	60546.12723	66	1.481194	0.447420	0.102737	19.085361	0.001247	0.015581	...	NaN	NaN	0	0.496874	0.493693	1.690698e-01	18.293407	4.123830e-03	3.045361e-02	38
3	3436.664307	60362.00000	60369.30842	66	2.093236	0.617962	0.031888	17.350933	0.000167	0.006134	...	NaN	NaN	0	1.211593	0.586116	4.590447e-02	16.562313	2.647262e-04	7.582015e-03	79
4	3325.815674	60476.00000	60476.15480	66	2.019966	0.590387	0.057489	18.055870	0.000590	0.010950	...	NaN	NaN	0	0.913779	0.598814	6.471661e-02	17.270544	5.449418e-04	1.152790e-02	49
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
44606	2897.801514	60574.00000	60580.30096	98	0.261178	0.202494	0.468495	21.202513	0.042931	0.101005	...	NaN	NaN	0	NaN	NaN	NaN	NaN	NaN	NaN	0
44607	500.044525	62742.18036	62735.18036	2150	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	0	NaN	NaN	NaN	NaN	NaN	NaN	0
44608	3341.052979	60219.00000	60226.23987	102	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	0	NaN	NaN	NaN	NaN	NaN	NaN	0
44609	2105.805664	61488.00000	61493.23812	98	1.452305	0.541928	0.451856	20.889555	0.032904	0.080253	...	NaN	NaN	0	NaN	NaN	NaN	NaN	NaN	NaN	0
44610	1859.874756	61727.00000	61731.34561	70	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	0	0.859407	-767251.750000	9.258375e+11	36.474541	-1.229126e+18	1.327583e+06	6

44611 rows × 55 columns

Hay 44611 filas.

Opcional: usar las funciones info o describe sobre el dataframe de pandas para aprender más sobre los resultados obtenidos.

In [19]:

# df.info()

In [20]:

# df.describe()

4.4. Graficar un diagrama color-color¶

En el dataframe mostrado arriba, pareciera que para varios objetos de SSObject los ajustes de las curvas de fase para obtener magnitudes absolutas no fueron exitosos, ya que las columnas <banda>_H (donde <banda> es una banda u, g, r, i, z o y) tienen valores NaN. (Recordar en particular que no se simularon datos para las bandas u o y en DP0.3).

Antes de calcular y graficar los colores, quitar todas las filas para las cuales los ajustes de las curvas de fase no fueron exitosos para las bandas g, r, i, y/o z.

In [21]:

df.dropna(subset=['g_H', 'r_H', 'i_H', 'z_H'], inplace=True)
df.reset_index(inplace=True)
print('Número de filas después de eliminar filas: %d' % len(df))

Número de filas después de eliminar filas: 30853

Para objetos del Sistema Solar, las magnitudes absolutas se definen como la magnitud de un objeto a 1 ua del Sol y a 1 ua del observador, y en un ángulo de fase (el ángulo Sol-objeto-Tierra) de 0 grados.

Las magnitudes absolutas se obtienen corrigiendo por distancia, ajustando una función a la relación entre la magnitud absoluta y la fase y evaluando la función en una fase de 0 grados.

El proceso de ajuste de las curvas de fase será abordado en otro tutorial.

Calcular los colores en los filtros de Rubin para los objetos SSObject que tienen magnitudes absolutas.

In [22]:

df['gr'] = df['g_H'] - df['r_H']
df['ri'] = df['r_H'] - df['i_H']
df['iz'] = df['i_H'] - df['z_H']

Graficar diagramas color-color como histogramas bi-dimensionales (mapas de calor -- heatmaps).

In [23]:

fig, ax = plt.subplots(1, 2, figsize=(8, 3))

h = ax[0].hist2d(df['gr'], df['ri'],
                 bins=(np.arange(-0.5, 1.5, 0.01),
                       np.arange(-0.5, 0.75, 0.01)),
                 norm='log')
ax[0].set_xlabel('g - r')
ax[0].set_ylabel('r - i')
ax[0].grid()
plt.colorbar(h[3])

h = ax[1].hist2d(df['ri'], df['iz'],
                 bins=(np.arange(-0.5, 0.75, 0.01),
                       np.arange(-0.75, 0.5, 0.01)),
                 norm='log')
ax[1].set_xlabel('r - i')
ax[1].set_ylabel('i - z')
ax[1].grid()
plt.colorbar(h[3])

fig.suptitle('Diagramas Color-Color para el catálogo SSObject')
fig.tight_layout()
plt.show()

Hay dos colores usados en la simulación - pero este no es el caso para los objetos reales del Sistema Solar. Estos gráficos se verán muy diferentes en el futuro, cuando estén hechos con datos reales de Rubin.

5. La tabla `SSSource`¶

Como se describe arriba, el Procesamiento del Sistema Solar vincula las fuentes de DiaSource (detecciones en las imágenes diferenciales individuales) con objetos móviles en los objetos de SSObject.

La tabla SSSource (Fuentes del Sistema Solar) contiene las coordenadas 2-d (cielo) y las distancias y velocidades 3-d para cada objeto de SSObject al momento de cada observación de LSST de ese objeto SSObject. Las tablas SSSourcey DiaSource son 1:1, ya que cada una contiene datos por observación, mientras que SSObject contiene datos por objeto.

5.1. Tamaño¶

Puede tomar hasta un minuto obtener el tamaño del catálogo SSSource.

In [24]:

results = service.search("SELECT COUNT(*) from dp03_catalogs_10yr.SSSource")
results.to_table().to_pandas()

Out[24]:

	COUNT1
0	653005444

¡Esta tabla contiene más de 653 millones de fuentes!

5.2. Columnas¶

Opcional: imprimir la información de las columnas para la tabla SSSource.

In [25]:

# results = service.search("SELECT column_name, datatype, description, "
#                          "unit from TAP_SCHEMA.columns "
#                          "WHERE table_name = 'dp03_catalogs_10yr.SSSource'")
# results.to_table().to_pandas()

5.3. Obtener datos de un objeto de `SSObject`¶

Es posible obtener los datos de SSSource para un conjunto de objetos de SSObject. Por ejemplo, para obtener todas las fuentes de SSSource para los objetos de SSObject obtenidos en la Sección 4, usar una consulta como:

SELECT * FROM dp03_catalogs_10yr.SSSource
WHERE ssObjectId BETWEEN 9038903462544093184 AND 9223370875126069107

Sin embargo, la mejor manera de mostrar los datos en la tabla SSSource es mirar sólo un objeto de SSObject, y el que tiene un ssObjectId = 6793512588511170680 es una opción divertida.

Obtener las coordenadas X e Y heliocéntricas (centradas en el sol) y topocéntricas (centradas en el observatorio).

In [26]:

results = service.search("SELECT heliocentricX, heliocentricY, "
                         "topocentricX, topocentricY, ssObjectId "
                         "FROM dp03_catalogs_10yr.SSSource "
                         "WHERE ssObjectId = 6793512588511170680")
df_xy = results.to_table().to_pandas()
print('Retrieved ', len(df_xy), ' rows.')

Retrieved  319  rows.

Opcional: mostrar la tabla en su totalidad o utilizar las funciones info o describe.

In [27]:

# df_xy

In [28]:

# df_xy.info()

In [29]:

# df_xy.describe()

5.3. Graficar las posiciones de un objeto de `SSObject`¶

Graficar las posiciones del objeto seleccionado de la tabla SSObject en el momento de cada observación de LSST utilizando las coordenadas X e Y heliocéntricas (centradas en el Sol; estrella naranja) y las coordenadas topocéntricas (centradas en el observatorio; círculo azul). Esto puede considerarse una proyección de la órbita en el plano del Sistema Solar.

Notar cómo los puntos no están espaciados regularmente. Esto es debido a que hay un punto por observación realizada por el LSST al objeto y algunos años éste recibe más o menos observaciones.

Observar cómo los puntos aparecen en una elipse alrededor del Sol con coordenadas heliocéntricas (izquierda). Esto se debe a que el objeto seleccionado está en el cinturón principal de asteroides y lo suficientemente cerca del Sol como para completar al menos una órbita durante los 10 años de estudio del LSST. Si el objeto seleccionado hubiese estado en el Sistema Solar exterior, o si este tutorial estuviera utilizando el conjunto de datos de 1 año, el gráfico a continuación mostraría un arco en lugar de una elipse.

El gráfico de las coordenadas topocéntricas (derecha) no se ve elíptico porque el movimiento de la Tierra respecto al objeto durante los 10 años del LSST está reflejado en los datos. Para las coordenadas topocéntricas, la rotación de la Tierra también contribuye, pero es un efecto mucho menor en la escala de estos gráficos.

In [30]:

fig, ax = plt.subplots(1, 2, figsize=(8, 4))
ax[0].grid()
ax[0].plot(df_xy['heliocentricX'], df_xy['heliocentricY'],
           'o', ms=4, mew=0, color='black')
ax[0].plot(0, 0, '*', ms=15, color='darkorange')
ax[0].set_xlabel('X heliocéntrico (au)')
ax[0].set_ylabel('Y heliocéntrico (au)')
ax[0].set_title('Coordenadas Heliocéntricas')
ax[1].grid()
ax[1].plot(df_xy['topocentricX'], df_xy['topocentricY'],
           'o', ms=4, mew=0, color='black')
ax[1].plot(0, 0, 'o', ms=15, color='dodgerblue')
ax[1].set_xlabel('X topocéntrico (au)')
ax[1].set_ylabel('Y topocéntrico (au)')
ax[1].set_title('Coordenadas Topocéntricas')
fig.suptitle('Trayectoria XY de 6793512588511170680')
fig.tight_layout()
plt.show()

6. El catálogo `DiaSource`¶

Por último, pero definitivamente no menos importante, la tabla DiaSource - la cual en realidad es la primera en ser generada por el pipeline del Procesamiento de Datos Inmediatos. Esta tabla va a contener mediciones para todas las fuentes detectadas con una relación señal/ruido de al menos 5 en una imagen diferencial, incluyendo fuentes en movimiento y estáticas en el cielo (static-sky).

Sin embargo, para la simulación DP0.3 la tabla DiaSource contiene sólo objetos en movimiento, y no objetos transitorios con posición fija en el cielo (ni defectos del detector).

Para DP0.3 la tabla simulada DiaSource contiene sólo un subconjunto de las columnas que tendrá la tabla DiaSource real; ver Tabla 1 del Documento de Definición de Productos de Datos de Rubin. Más aún, la tabla DiaSource contiene unas pocas columnas con valores de verdad extras, tales como raTrue, decTrue, magTrue.

Para DP0.3 no se simularon variaciones fotométricas debidas a, e.g., la rotación de cuerpos no esféricos, colisiones o eventos de desgasificación. Toda la evolución en magnitudes aparentes con el tiempo son debidas a la curva de fase, la cual es explorada en otro tutorial.

6.1. Tamaño¶

Opcional: obtener el tamaño de la tabla DiaSource (es 1:1 con la tabla SSSource, por lo que tienen el mismo tamaño).

In [31]:

# results = service.search("SELECT COUNT(*) from dp03_catalogs_10yr.DiaSource")
# results.to_table().to_pandas()

6.2. Columnas¶

Opcional: imprimir la información de las columnas.

In [32]:

# results = service.search("SELECT column_name, datatype, description, "
#                          "unit from TAP_SCHEMA.columns "
#                          "WHERE table_name = 'dp03_catalogs_10yr.DiaSource'")
# results.to_table().to_pandas()

6.3. Obtener datos para un objeto de `SSObject`¶

Parecido a lo que hicimos con los datos de SSSource en la Sección 5, es posible obtener los datos de DiaSource para un conjunto de objetos SSObject. Por ejemplo, para obtener todas las fuentes en DiaSource para los objetos de SSObject que se obtuvieron en la Sección 4, usar una consulta como:

SELECT * FROM dp03_catalogs_10yr.DiaSource
WHERE ssObjectId BETWEEN 9038903462544093184 AND 9223370875126069107

Sin embargo, como vimos para los datos de las fuentes de SSSource en la Sección 5, la mejor forma de mostrar los datos en la tabla DiaSource es mirar sólo un objeto de SSObject. Aquí, como en la Sección 5, usaremos el que tiene un ssObjectId = 6793512588511170680.

Obtener la ascensión recta, declinación y el TAI (Tiempo Atómico Internacional) correctos en el punto medio de exposición (MJD) de la observación.

In [33]:

results = service.search("SELECT ra, dec, midpointMjdTai "
                         "FROM dp03_catalogs_10yr.DiaSource "
                         "WHERE ssObjectId = 6793512588511170680")
df = results.to_table().to_pandas()

df

Out[33]:

	ra	dec	midpointMjdTai
0	348.416159	-2.408505	62096.01616
1	351.997357	1.212909	62043.12454
2	53.506941	19.244016	63430.44136
3	182.188871	-1.424286	62709.97076
4	137.800325	7.798673	62560.10471
...	...	...	...
314	312.379754	-12.443726	63148.08742
315	313.660126	-12.585707	63160.04304
316	147.164625	4.901712	62510.30546
317	220.150693	-12.830205	61555.00272
318	17.293268	9.615766	60830.43848

319 rows × 3 columns

6.4. Graficar datos sobre transitorios para un objeto de `SSObject`¶

Aquí, graficaremos las coordenadas ecuatoriales para todas las observaciones, coloreadas por la MJD de la observación.

Este gráfico muestra la posición del objeto en el cielo, visto desde la Tierra. Recordar de los gráficos de la Sección 5 que este objeto está aproximadamente a sólo unas 2 ua del Sol, por lo que se encuentra en el borde interno del Cinturón Principal de Asteroides. La paralaje de la órbita terrestre contribuye a cómo se ve este gráfico.

In [34]:

cmap = colormaps['viridis']
fig, ax = plt.subplots(1, 1, figsize=(12, 3))
im = ax.scatter(df['ra'], df['dec'], c=df['midpointMjdTai'], cmap=cmap)
ax.set_xlabel('ra (deg)')
ax.set_ylabel('dec (deg)')
fig.suptitle("Trayectoria de un objeto SSObject en coordenadas ecuatoriales a lo largo de 10 años")
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.02, 0.7])
fig.colorbar(im, cax=cbar_ax, label='TAI')
ax.grid(True)

El uso de la magnitud y el filtro de los datos obtenidos de la tabla DiaSource es dejado como ejercicio en la Sección 7.

7. Ejercicios¶

1. Calcular el semieje mayor para un subconjunto de objetos.

El elemento orbital del semieje mayor ($a$) no está pre-computado en la tabla MPCORB debido a que se puede obtener a partir de la elipticidad de la órbita ($e$) y de la distancia al perihelio ($q$) con $a = q/(1-e)$.

Copiar la consulta de la Sección 3.3 y modificarla para obtener sólo $e$ y $q$ para un subconjunto de objetos. Luego agregar una columna a la tabla resultante para el semieje mayor y graficar la excentricidad versus el semieje mayor.

Pistas:

Restringir la consulta a objetos con excentricidades entre 0 y 1.
df['a'] = df['q'] / (1.0 - df['e'])

2. ¿Cómo cambia la magnitud de un objeto de la tabla SSObject con el tiempo?

En la Sección 6.4, no se obtuvieron ni la magnitud ni el filtro para todas las detecciones del objeto de SSObject en la tabla DiaSource; modificar la consulta para obtener mag y band, y graficar la magnitud como función del tiempo para el filtro de tu elección (recordar: no habrá magnitudes en las bandas u o y). Notar que la única causa simulada de los cambios fotométricos en los objetos de DP0.3 es la distancia a la Tierra y la curva de fase. La corrección por distancia y las aplicaciones de curvas de fase, se abordan en otro tutorial.