Как построить данные астрономии Gaia для изображений TESS, используя Python?
Короче говоря: я хочу представить данные астрономии Гайи изображениям TESS в Python. Как это возможно? Смотрите ниже для подробной версии.
У меня есть изображение TESS размером 64x64 пикселя со звездой Gaia ID 4687500098271761792. На странице 8 Руководства обсерватории TESS говорится, что 1 пиксель составляет ~ 21 угловую секунду. Используя архив Gaia, я ищу эту звезду (под верхними объектами, нажимаю "Поиск") и отправляю запрос, чтобы увидеть звезды в пределах 1000 угловых секунд, примерно того радиуса, который нам нужен. Имя, которое я использую для поиска, - Gaia DR2 4687500098271761792, как показано ниже:
![enter image description here]()
Отправьте запрос, и я получу список из 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)
и получить хорошую картинку:
![enter image description here]()
и кучу предупреждений, большая часть которых была любезно объяснена здесь (предупреждения в Q, объяснение в комментариях).
Обратите внимание, что действительно хорошо, что мы используем проекцию WCS. Чтобы проверить, давайте просто hdul.data данные в hdul.data не заботясь о проекции:
plt.imshow(hdul.data)
Результат:
![enter image description here]()
Почти так же, как и раньше, но теперь метки осей - это просто числа пикселей, а не 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='+')
![enter image description here]()
Форма не круг, как ожидалось. Это может быть индикатором будущих неприятностей.
Попробуем преобразовать эти значения 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')
Результат:
![enter image description here]()
Функция 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')
Мы получаем:
![enter image description here]()
что не близко к тому, что было бы идеально. Я ожидаю, что на нем будет 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='+')
и полученный сюжет:
![enter image description here]()
Этот изгибный крест ожидается, так как 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')
Результат обнадеживает:
![enter image description here]()
(Смотри блокнот здесь.)
Ответы
Ответ 1
Сначала я должен сказать, отличный вопрос! Очень подробный и воспроизводимый. Я прошел ваш вопрос и попытался повторить упражнение, начиная с вашего git-репо и загружая каталог из архива GAIA.
РЕДАКТИРОВАТЬ
Программно ваш код в порядке (см. СТАРУЮ ЧАСТЬ ниже для немного другого подхода). Проблема с отсутствующими точками заключается в том, что при загрузке файла CSV из архива GAIA можно получить только 500 точек данных. Поэтому все выглядит так, как будто все точки из запроса забиты в странную форму. Однако, если вы ограничите радиус поиска меньшим значением, вы увидите, что в изображении TESS есть точки:
![enter image description here]()
пожалуйста, сравните с версией, показанной ниже в СТАРОЙ ЧАСТИ. Код такой же, как и ниже, только загруженный 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")
Однако это приводит к сюжету, похожему на ваш:
![enter image description here]()
Чтобы проверить мое подозрение, что координаты из архива 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()
![First 1000 rows, ordered by random index]()
Обратите внимание, что я добавил порядок с помощью 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 строки.
![All sources. 103574 rows]()
