Как построить данные астрономии 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 не заботясь о проекции:
plt.imshow(hdul.data)
Результат:
Почти так же, как и раньше, но теперь метки осей - это просто числа пикселей, а не RA и DEC, как было бы предпочтительнее. DEC а также RA значения на первом графике составляют около -72° и 16° соответственно, что хорошо, учитывая, что каталог Gaia дал нам звезды в непосредственной близости от 4687500098271761792 примерно с этими координатами. Таким образом, прогноз, кажется, достаточно хорошо.
Теперь давайте попробуем нанести звезды Гайи на 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, в котором я работал, доступен здесь.
Какие шаги я пропускаю или что я делаю не так?
Дальнейшие разработки
В ответ на другой вопрос keflavich любезно предложил использовать 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, которые находятся в пределах области, которую вы искали. Облако данных, кажется, даже имеет квадратную границу внутри, которая может исходить из чего-то квадратного поля зрения.
Действительно хороший пример. Просто отметим, что вы также можете интегрировать запрос в архив 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,5 минут). Полный запрос вернет 103574 строки.
