¿Cómo trazar los datos de astrometría de Gaia en imágenes TESS utilizando Python?

Larga historia corta: quiero trazar datos de astrometría de Gaia a imágenes TESS en Python. ¿Como es posible? Vea a continuación la versión elaborada.


Tengo imágenes TESS de 64×64 píxeles de una estrella con ID de Gaia 4687500098271761792 . La página 8 de la Guía del Observatorio TESS dice que 1 píxel es ~ 21 arcsec. Utilizando el Archivo de Gaia , busco esta estrella (debajo de las funciones principales , haga clic en Buscar ) y envío una consulta para ver las estrellas dentro de 1000 arcsec, aproximadamente el radio que necesitamos. El nombre que uso para la búsqueda es Gaia DR2 4687500098271761792 , como se muestra a continuación:

introduzca la descripción de la imagen aquí

Enviar consulta y obtengo una lista de 500 estrellas con coordenadas RA y DEC . Seleccionando CSV y Download results , obtengo la lista de estrellas alrededor de 4687500098271761792 . Este archivo resultante también se puede encontrar aquí . Esta es la entrada de Gaia que queremos usar.

Desde el TESS, tenemos 4687500098271761792_med.fits , un archivo de imagen. Lo graficamos utilizando:

 from astropy.io import fits from astropy.wcs import WCS import matplotlib.pyplot as plt hdul = fits.open("4687500098271761792_med.fits")[0] wcs = WCS(hdul.header) fig = plt.figure(figsize=(12,12)) fig.add_subplot(111, projection=wcs) plt.imshow(hdul.data) 

y obtener una buena foto:

introduzca la descripción de la imagen aquí

y un montón de advertencias, la mayoría de las cuales se explicaron amablemente aquí (advertencias en la Q, explicación en los comentarios) .

Tenga en cuenta que es bueno que estemos usando la proyección WCS . Para comprobarlo, solo hdul.data los datos en hdul.data sin preocuparnos por la proyección:

 plt.imshow(hdul.data) 

El resultado:

introduzca la descripción de la imagen aquí

Casi lo mismo que antes, pero ahora las tags de los ejes son solo números de píxeles, no RA y DEC , como sería preferible. Los valores de DEC y RA en el primer gráfico son aproximadamente -72 ° y 16 ° respectivamente, lo que es bueno, dado que el catálogo de Gaia nos dio estrellas en la proximidad de 4687500098271761792 con aproximadamente estas coordenadas. Así que la proyección parece estar razonablemente bien.

Ahora intentemos trazar las estrellas de Gaia sobre la ttwig imshow() . Leemos el archivo CSV que descargamos anteriormente y extraemos los valores RA y DEC de los objetos:

 import pandas as pd df=pd.read_csv("4687500098271761792_within_1000arcsec.csv") ralist=df['ra'].tolist() declist=df['dec'].tolist() 

Parcela para comprobar:

 plt.scatter(ralist,declist,marker='+') 

introduzca la descripción de la imagen aquí

La forma no es un círculo como se esperaba. Esto podría ser un indicador de problemas futuros.

Intentemos transformar estos valores de RA y DEC en WCS , y trazarlos de esa manera:

 for index, each in enumerate(ralist): ra, dec = wcs.all_world2pix([each], [declist[index]], 1) plt.scatter(ra, dec, marker='+', c='k') 

El resultado es:

introduzca la descripción de la imagen aquí

La función all_world2pix es desde aquí . El parámetro 1 solo establece donde está el origen. La descripción de all_world2pix dice:

Aquí, el origen es la coordenada en la esquina superior izquierda de la imagen. En los estándares FITS y Fortran, esto es 1. En los estándares Numpy y C, este es 0.

Sin embargo, la forma de la distribución de puntos que obtenemos no es nada prometedora. Vamos a juntar los datos de TESS y Gaia:

 hdul = fits.open("4687500098271761792_med.fits")[0] wcs = WCS(hdul.header) fig = plt.figure(figsize=(12,12)) fig.add_subplot(111, projection=wcs) plt.imshow(hdul.data) for index, each in enumerate(ralist): ra, dec = wcs.all_world2pix([each], [declist[index]], 1) plt.scatter(ra, dec, marker='+', c='k') 

Obtenemos:

introduzca la descripción de la imagen aquí

Lo que no está ni cerca de lo que sería ideal. Espero tener una imagen subyacente de imshow() con muchos marcadores, y los marcadores deben estar donde están las estrellas en la imagen de TESS. El cuaderno de Jupyter en el que trabajé está disponible aquí .

¿Qué pasos me faltan o qué hago mal?


Nuevos desarrollos

En respuesta a otra pregunta , Keflavich amablemente sugirió usar un argumento de transform para trazar en coordinados mundiales . Lo probé con algunos puntos de ejemplo (la cruz doblada en la gráfica a continuación). También trazaron los datos de Gaia en la ttwig sin procesarlos, terminaron concentrados en un espacio muy estrecho. Aplicado a ellos el método de transform , obtuvo un resultado aparentemente muy similar al anterior. El código (y también aquí ):

 import pandas as pd df=pd.read_csv("4687500098271761792_within_1000arcsec.csv") ralist=df['ra'].tolist() declist=df['dec'].tolist() from astropy.io import fits from astropy.wcs import WCS import matplotlib.pyplot as plt hdul = fits.open("4687500098271761792_med.fits")[0] wcs = WCS(hdul.header) fig = plt.figure(figsize=(12,12)) fig.add_subplot(111, projection=wcs) plt.imshow(hdul.data) ax = fig.gca() ax.scatter([16], [-72], transform=ax.get_transform('world')) ax.scatter([16], [-72.2], transform=ax.get_transform('world')) ax.scatter([16], [-72.4], transform=ax.get_transform('world')) ax.scatter([16], [-72.6], transform=ax.get_transform('world')) ax.scatter([16], [-72.8], transform=ax.get_transform('world')) ax.scatter([16], [-73], transform=ax.get_transform('world')) ax.scatter([15], [-72.5], transform=ax.get_transform('world')) ax.scatter([15.4], [-72.5], transform=ax.get_transform('world')) ax.scatter([15.8], [-72.5], transform=ax.get_transform('world')) ax.scatter([16.2], [-72.5], transform=ax.get_transform('world')) ax.scatter([16.6], [-72.5], transform=ax.get_transform('world')) ax.scatter([17], [-72.5], transform=ax.get_transform('world')) for index, each in enumerate(ralist): ax.scatter([each], [declist[index]], transform=ax.get_transform('world'),c='k',marker='+') for index, each in enumerate(ralist): ax.scatter([each], [declist[index]],c='b',marker='+') 

y la ttwig resultante:

introduzca la descripción de la imagen aquí

Se espera esta cruz, ya que TESS no está alineada con las líneas de latitud y longitud constantes (es decir, los arms de la cruz no tienen que estar paralelos a los lados de la imagen TESS, trazados con imshow() ). Ahora tratemos de trazar líneas constantes de RA y DEC (o, por ejemplo, líneas de latitud y longitud constantes) para comprender mejor por qué los puntos de datos de Gaia están fuera de lugar. Expande el código de arriba por unas pocas líneas:

 ax.coords.grid(True, color='green', ls='solid') overlay = ax.get_coords_overlay('icrs') overlay.grid(color='red', ls='dotted') 

El resultado es alentador:

introduzca la descripción de la imagen aquí

(Ver cuaderno aquí ).

Primero tengo que decir, gran pregunta! Muy detallado y reproducible. Repasé tu pregunta e intenté rehacer el ejercicio a partir de tu repository de git y descargando el catálogo del archivo GAIA.

EDITAR

Programáticamente, su código está bien (vea ANTIGUO PARTE a continuación para un enfoque ligeramente diferente). El problema con los puntos faltantes es que solo se obtienen 500 puntos de datos al descargar el archivo csv desde el archivo GAIA. Por lo tanto, parece que todos los puntos de la consulta están en una forma extraña. Sin embargo, si restringe el radio de la búsqueda a un valor más pequeño, puede ver que hay puntos que se encuentran dentro de la imagen TESS:

introduzca la descripción de la imagen aquí

por favor, compare con la versión que se muestra a continuación en la PARTE ANTIGUA. El código es el mismo que el de abajo, solo el archivo csv descargado es para un radio más pequeño. Por lo tanto, parece que acaba de descargar una parte de todos los datos disponibles del archivo GAIA al exportar a csv. La forma de sortear esto es hacer la búsqueda como lo hiciste. Luego, en la página de resultados, haga clic en Show query in ADQL form en la parte inferior y en la consulta que se muestra en el cambio de formato SQL:

 Select Top 500 

a

 Select 

Al comienzo de la consulta.

PARTE ANTIGUA (el código está bien y funcionando, pero mi conclusión es incorrecta):

Para el trazado, utilicé aplpy (uso matplotlib en segundo plano) y terminé con el siguiente código:

 from astropy.io import fits from astropy.wcs import WCS import aplpy import matplotlib.pyplot as plt import pandas as pd from astropy.coordinates import SkyCoord import astropy.units as u from astropy.io import fits fits_file = fits.open("4687500098271761792_med.fits") central_coordinate = SkyCoord(fits_file[0].header["CRVAL1"], fits_file[0].header["CRVAL2"], unit="deg") figure = plt.figure(figsize=(10, 10)) fig = aplpy.FITSFigure("4687500098271761792_med.fits", figure=figure) cmap = "gist_heat" stretch = "log" fig.show_colorscale(cmap=cmap, stretch=stretch) fig.show_colorbar() df = pd.read_csv("4687500098271761792_within_1000arcsec.csv") # the epoch found in the dataset is J2015.5 df['coord'] = SkyCoord(df["ra"], df["dec"], unit="deg", frame="icrs", equinox="J2015.5") coords = df["coord"].tolist() coords_degrees = [[coord.ra.degree, coord.dec.value] for coord in df["coord"]] ra_values = [coord[0] for coord in coords_degrees] dec_values = [coord[1] for coord in coords_degrees] width = (40*u.arcmin).to(u.degree).value height = (40*u.arcmin).to(u.degree).value fig.recenter(x=central_coordinate.ra.degree, y=central_coordinate.dec.degree, width=width, height=height) fig.show_markers(central_coordinate.ra.degree,central_coordinate.dec.degree, marker="o", c="white", s=15, lw=1) fig.show_markers(ra_values, dec_values, marker="o", c="blue", s=15, lw=1) fig.show_circles(central_coordinate.ra.degree,central_coordinate.dec.degree, radius=(1000*u.arcsec).to(u.degree).value, edgecolor="black") fig.save("GAIA_TESS_test.png") 

Sin embargo, esto se traduce en una ttwig similar a la suya:

introduzca la descripción de la imagen aquí

Para comprobar mi sospecha de que las coordenadas del archivo GAIA se muestran correctamente, dibujo un círculo de 1000 arcsec desde el centro de la imagen TESS. Como puede ver, se alinea aproximadamente con la forma circular del lado exterior (visto desde el centro de la imagen) de la nube de puntos de datos de las posiciones GAIA. Simplemente creo que estos son todos los puntos en el archivo GAIA DR2 que se encuentran dentro de la región que buscó. La nube de datos incluso parece tener un límite cuadrado en el interior, que podría provenir de algo como un campo de visión cuadrado.