如何使用 Python 将 Gaia 天体测量数据绘制成 TESS 图像？
How to plot Gaia astrometry data to TESS images using Python?
长话短说：我想将 Gaia 天体测量数据绘制到 Python 中的 TESS 图像中。这怎么可能？
详见下文。
我有 64x64 像素 TESS imagery of a star with Gaia ID 4687500098271761792. Page 8 of the TESS Observatory Guide says 1 pixel is ~21 arcsec. Using the Gaia Archive，我搜索这颗星（在 top features 下方，点击 Search。）然后提交查询以查看 1000 弧秒内的星星，大致是我们需要的半径。我用来搜索的名字是Gaia DR2 4687500098271761792，如下图：
提交查询，我得到一个包含 500 颗星的列表，坐标为 RA 和 DEC。 Select CSV 和 Download results，我得到了 4687500098271761792 周围的星星列表。这个生成的文件也可以找到 here. This is the input from 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)
拍一张漂亮的照片：
和一堆警告，其中大部分都得到了很好的解释 。
请注意，我们使用 WCS 投影确实很好。为了检查，让我们只绘制 hdul.data 中的数据而不关心投影：
plt.imshow(hdul.data)
结果：
几乎和以前一样，但现在轴的标签只是像素数，而不是 RA and DEC，这将是更可取的。第一张图中的 DEC 和 RA 值分别在 -72° 和 16° 左右，这很好，因为盖亚目录给了我们 4687500098271761792[ 附近的恒星=134=] 大概有这些坐标。所以投影似乎还不错。
现在让我们尝试在 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 来自 here。 1 参数只是设置原点的位置。 all_world2pix 的描述说：
Here, origin is the coordinate in the upper left corner of the image.
In FITS and Fortran standards, this is 1. In Numpy and C standards
this is 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 Notebook 可用 here.
我错过了什么步骤，或者我做错了什么？
进一步发展
在 to another , keflavich kindly suggested to use a transform argument for plotting in world coordintes. Tried it with some example points (the bended cross on the plot below). Also plotted the Gaia data on the plot without processing it, they ended up being concentrated at a very narrow space. Applied to them the transform method, got a seemingly very similar result than before. The code ( & also here):
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() 绘制).现在让我们尝试绘制恒定的赤经和赤纬线（或者说，恒定的纬度和经度线）以更好地理解为什么来自盖亚的数据点放错了位置。将上面的代码扩展几行：
ax.coords.grid(True, color='green', ls='solid')
overlay = ax.get_coords_overlay('icrs')
overlay.grid(color='red', ls='dotted')
结果令人鼓舞：
（参见笔记本 here。）
首先我要说，好问题！非常详细和可重现。我审阅了你的问题并尝试从你的 git 存储库开始重做练习并从 GAIA 存档下载目录。
编辑
您的代码在编程上没有问题（请参阅下面的旧部分，了解稍微不同的方法）。缺失点的问题是从 GAIA 存档下载 csv 文件时只能获得 500 个数据点。因此，看起来查询中的所有点都被塞进了一个奇怪的形状。但是，如果您将搜索半径限制为较小的值，您会看到有些点位于 TESS 图像内：
请与以下旧部分中显示的版本进行比较。代码与下面相同，只是下载的 csv 文件用于较小的半径。因此，在导出到 csv 时，您似乎只是从 GAIA 存档中下载了所有可用数据的一部分。避免这种情况的方法是像您一样进行搜索。然后，在结果页面上单击底部的 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 档案中的坐标是否正确显示的怀疑，我从 TESS 图像的中心画了一个 1000 arcsec 的圆。如您所见，它与 GAIA 位置数据点云的外侧（从图像中心看）的圆形形状大致对齐。我只是认为这些都是 GAIA DR2 存档中属于您搜索的区域的所有点。数据云的内部甚至似乎有一个方形的边界，这可能来自方形视野之类的东西。
很好的例子。顺便提一下，您还可以使用 astropy
中包含的 astroquery.gaia 模块将查询集成到 Gaia 档案中
https://astroquery.readthedocs.io/en/latest/gaia/gaia.html
通过这种方式，您将能够运行 Gaia 存档中的相同查询 UI 并以更简单的方式更改为不同的来源
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 添加了订单，这将消除这种奇怪的 non-circular 行为。该索引对于不强制初始测试的完整输出非常有用。
此外，我已将Simbad的赤经坐标输出声明为小时。
最后，我使用了异步查询，它对执行时间和响应中的最大行数限制较少。
您也可以将查询更改为
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 分钟）。完整查询将 return 103574 行。
长话短说：我想将 Gaia 天体测量数据绘制到 Python 中的 TESS 图像中。这怎么可能？ 详见下文。
我有 64x64 像素 TESS imagery of a star with Gaia ID 4687500098271761792. Page 8 of the TESS Observatory Guide says 1 pixel is ~21 arcsec. Using the Gaia Archive，我搜索这颗星（在 top features 下方，点击 Search。）然后提交查询以查看 1000 弧秒内的星星，大致是我们需要的半径。我用来搜索的名字是Gaia DR2 4687500098271761792，如下图：
提交查询，我得到一个包含 500 颗星的列表，坐标为 RA 和 DEC。 Select CSV 和 Download results，我得到了 4687500098271761792 周围的星星列表。这个生成的文件也可以找到 here. This is the input from 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)
拍一张漂亮的照片：
和一堆警告，其中大部分都得到了很好的解释
请注意，我们使用 WCS 投影确实很好。为了检查，让我们只绘制 hdul.data 中的数据而不关心投影：
plt.imshow(hdul.data)
结果：
几乎和以前一样，但现在轴的标签只是像素数，而不是 RA and DEC，这将是更可取的。第一张图中的 DEC 和 RA 值分别在 -72° 和 16° 左右，这很好，因为盖亚目录给了我们 4687500098271761792[ 附近的恒星=134=] 大概有这些坐标。所以投影似乎还不错。
现在让我们尝试在 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 来自 here。 1 参数只是设置原点的位置。 all_world2pix 的描述说：
Here, origin is the coordinate in the upper left corner of the image. In FITS and Fortran standards, this is 1. In Numpy and C standards this is 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 Notebook 可用 here.
我错过了什么步骤，或者我做错了什么？
进一步发展
在transform argument for plotting in world coordintes. Tried it with some example points (the bended cross on the plot below). Also plotted the Gaia data on the plot without processing it, they ended up being concentrated at a very narrow space. Applied to them the transform method, got a seemingly very similar result than before. The code ( & also here):
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() 绘制).现在让我们尝试绘制恒定的赤经和赤纬线（或者说，恒定的纬度和经度线）以更好地理解为什么来自盖亚的数据点放错了位置。将上面的代码扩展几行：
ax.coords.grid(True, color='green', ls='solid')
overlay = ax.get_coords_overlay('icrs')
overlay.grid(color='red', ls='dotted')
结果令人鼓舞：
（参见笔记本 here。）
首先我要说，好问题！非常详细和可重现。我审阅了你的问题并尝试从你的 git 存储库开始重做练习并从 GAIA 存档下载目录。
编辑
您的代码在编程上没有问题（请参阅下面的旧部分，了解稍微不同的方法）。缺失点的问题是从 GAIA 存档下载 csv 文件时只能获得 500 个数据点。因此，看起来查询中的所有点都被塞进了一个奇怪的形状。但是，如果您将搜索半径限制为较小的值，您会看到有些点位于 TESS 图像内：
请与以下旧部分中显示的版本进行比较。代码与下面相同，只是下载的 csv 文件用于较小的半径。因此，在导出到 csv 时，您似乎只是从 GAIA 存档中下载了所有可用数据的一部分。避免这种情况的方法是像您一样进行搜索。然后，在结果页面上单击底部的 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 档案中的坐标是否正确显示的怀疑，我从 TESS 图像的中心画了一个 1000 arcsec 的圆。如您所见，它与 GAIA 位置数据点云的外侧（从图像中心看）的圆形形状大致对齐。我只是认为这些都是 GAIA DR2 存档中属于您搜索的区域的所有点。数据云的内部甚至似乎有一个方形的边界，这可能来自方形视野之类的东西。
很好的例子。顺便提一下，您还可以使用 astropy
中包含的 astroquery.gaia 模块将查询集成到 Gaia 档案中https://astroquery.readthedocs.io/en/latest/gaia/gaia.html
通过这种方式，您将能够运行 Gaia 存档中的相同查询 UI 并以更简单的方式更改为不同的来源
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 添加了订单，这将消除这种奇怪的 non-circular 行为。该索引对于不强制初始测试的完整输出非常有用。
此外，我已将Simbad的赤经坐标输出声明为小时。
最后，我使用了异步查询，它对执行时间和响应中的最大行数限制较少。
您也可以将查询更改为
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 分钟）。完整查询将 return 103574 行。