Как построить данные астрономии Gaia для изображений TESS, используя Python?

Короче говоря: я хочу представить данные астрономии Гайи изображениям TESS в Python. Как это возможно? Смотрите ниже для подробной версии.

У меня есть изображение TESS размером 64x64 пикселя со звездой Gaia ID 4687500098271761792. На странице 8 Руководства обсерватории TESS говорится, что 1 пиксель составляет ~ 21 угловую секунду. Используя архив Gaia, я ищу эту звезду (под верхними объектами, нажимаю "Поиск") и отправляю запрос, чтобы увидеть звезды в пределах 1000 угловых секунд, примерно того радиуса, который нам нужен. Имя, которое я использую для поиска, - Gaia DR2 4687500098271761792, как показано ниже:

Отправьте запрос, и я получу список из 500 звезд с координатами RA и DEC. Выберите CSV и Download results, я получаю список звезд вокруг 4687500098271761792. Этот результирующий файл также можно найти здесь. Это вклад от Gaia, который мы хотим использовать.

Из TESS у нас есть 4687500098271761792_med.fits, файл изображения. Мы строим это, используя:

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)

и получить хорошую картинку:

и кучу предупреждений, большая часть которых была любезно объяснена здесь (предупреждения в Q, объяснение в комментариях).

Обратите внимание, что действительно хорошо, что мы используем проекцию WCS. Чтобы проверить, давайте просто hdul.data данные в hdul.data не заботясь о проекции:

plt.imshow(hdul.data)

Результат:

Почти так же, как и раньше, но теперь метки осей - это просто числа пикселей, а не RA и DEC, как было бы предпочтительнее. Значения DEC и RA на первом графике составляют около -72 ° и 16 ° соответственно, что хорошо, учитывая, что каталог Gaia дал нам звезды в непосредственной близости от 4687500098271761792 примерно с этими координатами. Таким образом, прогноз, кажется, достаточно хорошо.

Теперь давайте попробуем нанести звезды Gaia над imshow(). Мы читаем в CSV файле, который мы скачали ранее, и извлекаем из него значения RA и DEC:

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

Участок для проверки:

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

Форма не круг, как ожидалось. Это может быть индикатором будущих неприятностей.

Попробуем преобразовать эти значения RA и DEC в WCS и построить их таким образом:

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

Результат:

Функция all_world2pix отсюда. Параметр 1 просто устанавливает, где находится источник. Описание all_world2pix гласит:

Здесь начало координат - это координата в верхнем левом углу изображения. В стандартах FITS и Fortran это значение равно 1. В стандартах Numpy и C это значение равно 0.

Тем не менее, форма распределения точек, которую мы получаем, не является многообещающей. Позвольте собрать воедино данные TESS и 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')

Мы получаем:

что не близко к тому, что было бы идеально. Я ожидаю, что на нем будет imshow() с множеством маркеров, и маркеры должны быть там, где звезды находятся на изображении TESS. Блокнот Jupyter, в котором я работал, доступен здесь.

Какие шаги я пропускаю или что я делаю не так?

Дальнейшие разработки

В ответ на другой вопрос Кефлавич любезно предложил использовать аргумент transform для построения координат в мировых координатах. Пробовал это с некоторыми примерами точек (изогнутый крест на графике ниже). Также нанесены данные Gaia на график без их обработки, в итоге они оказались сосредоточены в очень узком пространстве. Применив к ним метод transform, получил, казалось бы, очень похожий результат, чем раньше. Код (и также здесь):

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='+')

и полученный сюжет:

Этот изгибный крест ожидается, так как TESS не выровнен с линиями постоянной широты и долготы (т.е. плечи креста не должны быть параллельны сторонам изображения TESS, нанесенного с помощью imshow()). Теперь давайте попробуем построить постоянные линии RA и DEC (или, скажем, постоянные линии широты и долготы), чтобы лучше понять, почему точки данных из Gaia смещены. Разверните код выше несколькими строками:

ax.coords.grid(True, color='green', ls='solid')

overlay = ax.get_coords_overlay('icrs')
overlay.grid(color='red', ls='dotted')

Результат обнадеживает:

(Смотри блокнот здесь.)

Ответ 1

Сначала я должен сказать, отличный вопрос! Очень подробный и воспроизводимый. Я прошел ваш вопрос и попытался повторить упражнение, начиная с вашего git-репо и загружая каталог из архива GAIA.

РЕДАКТИРОВАТЬ

Программно ваш код в порядке (см. СТАРУЮ ЧАСТЬ ниже для немного другого подхода). Проблема с отсутствующими точками заключается в том, что при загрузке файла CSV из архива GAIA можно получить только 500 точек данных. Поэтому все выглядит так, как будто все точки из запроса забиты в странную форму. Однако, если вы ограничите радиус поиска меньшим значением, вы увидите, что в изображении TESS есть точки:

пожалуйста, сравните с версией, показанной ниже в СТАРОЙ ЧАСТИ. Код такой же, как и ниже, только загруженный CSV файл предназначен для меньшего радиуса. Поэтому кажется, что вы только что загрузили часть всех доступных данных из архива GAIA при экспорте в CSV. Способ обойти это - выполнить поиск, как вы. Затем на странице результатов нажмите " Show query in ADQL form внизу и в запросе, отображаемом в изменении формата SQL:

Select Top 500

Select

в начале запроса.

СТАРАЯ ЧАСТЬ (код в порядке и работает, но мой вывод неверен):

Для построения графика я использовал aplpy - использует matplotlib в фоновом режиме - и получил следующий код:

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")

Однако это приводит к сюжету, похожему на ваш:

Чтобы проверить мое подозрение, что координаты из архива GAIA отображаются правильно, я рисую круг в 1000 угловых секунд от центра изображения TESS. Как вы можете видеть, он примерно совпадает с круглой формой внешней (видимой из центра изображения) стороны облака точек данных позиций GAIA. Я просто думаю, что это все точки в архиве GAIA DR2, которые находятся в пределах области, которую вы искали. Облако данных даже, кажется, имеет квадратную границу внутри, которая может исходить из чего-то квадратного поля зрения.

Ответ 2

Действительно хороший пример. Просто отметим, что вы также можете интегрировать запрос в архив Gaia с помощью модуля astroquery.gaia, включенного в astropy.

https://astroquery.readthedocs.io/en/latest/gaia/gaia.html

Таким образом, вы сможете выполнять те же запросы, которые находятся внутри пользовательского интерфейса архива Gaia, и проще переходить на другие источники.

from astroquery.simbad import Simbad
import astropy.units as u
from astropy.coordinates import SkyCoord
from astroquery.gaia import Gaia

result_table = Simbad.query_object("Gaia DR2 4687500098271761792")
raValue = result_table['RA']
decValue = result_table['DEC']

coord = SkyCoord(ra=raValue, dec=decValue, unit=(u.hour, u.degree), frame='icrs')

query = """SELECT TOP 1000 * FROM gaiadr2.gaia_source 
           WHERE CONTAINS(POINT('ICRS',gaiadr2.gaia_source.ra,gaiadr2.gaia_source.dec), 
           CIRCLE('ICRS',{ra},{dec},0.2777777777777778))=1 ORDER BY random_index""".format(ra=str(coord.ra.deg[0]),dec=str(coord.dec.deg[0]))


job = Gaia.launch_job_async(query)  
r = job.get_results()

ralist = r['ra'].tolist()
declist = r['dec'].tolist()

import matplotlib.pyplot as plt
plt.scatter(ralist,declist,marker='+')
plt.show()

Обратите внимание, что я добавил порядок с помощью random_index, чтобы устранить это странное некруглое поведение. Этот индекс довольно полезен, чтобы не форсировать полный вывод для начальных тестов.

Кроме того, я объявил вывод координат в виде часов.

Наконец, я использовал асинхронный запрос, который имеет меньше ограничений по времени выполнения и максимум строк в ответе.

Вы также можете изменить запрос на

query = """SELECT * FROM gaiadr2.gaia_source 
               WHERE CONTAINS(POINT('ICRS',gaiadr2.gaia_source.ra,gaiadr2.gaia_source.dec), 
               CIRCLE('ICRS',{ra},{dec},0.2777777777777778))=1""".format(ra=str(coord.ra.deg[0]),dec=str(coord.dec.deg[0]))

(снятие ограничения до 1000 строк) (в этом случае использование случайного индекса необязательно) для получения полного ответа от сервера.

Конечно, этот запрос занимает некоторое время (около 1 минуты). Полный запрос вернет 103574 строки.