如何将Gaia天体测量数据绘制成TESS图像

Question

长话短说:我想用Python把Gaia天体测量数据绘制成TESS图像。这怎么可能呢？详细的版本见下文。

我有64x64像素的苔丝图像的一颗星与盖亚ID 4687500098271761792。“苔丝天文台指南”第8页说，1像素是21弧秒。使用盖亚档案馆，我搜索这个星型(在顶部功能下面，单击search )。并提交一个查询，在1000弧秒内看到恒星，大致相当于我们需要的半径。用于搜索的名称是Gaia DR2 4687500098271761792，如下所示：

​

​

提交查询，我得到一个具有RA和DEC坐标的500颗恒星的列表。选择CSV和Download results，我会得到大约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。

然而，我们得到的点分布的形状一点也不乐观。让我们把苔丝和盖亚的数据放在一起：

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图像上的星星位置。我工作过的木星笔记本是可用的这里。

我错过了什么步骤，或者我做错了什么？

进一步发展

在响应 to另一个问题中，凯夫利奇善意地建议为世界坐标标绘使用transform参数。尝试了一些例子点(在下面的情节上弯曲的十字)。还绘制了盖亚的数据，但没有处理它，他们最终集中在一个非常狭窄的空间。将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='+')

以及由此产生的阴谋：

​

​

这个弯曲的十字是预期的，因为苔丝没有与恒定的纬度和经线对齐(即十字架的手臂不必与苔丝图像的两侧平行，用imshow()绘制)。现在，让我们尝试绘制常量RA & DEC线(或者说，恒定纬度和经度线)，以更好地理解为什么来自Gaia的数据点被错误放置。将上面的代码展开几行：

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

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

其结果是令人鼓舞的：

​

​

(见笔记本电脑这里。)

Answer 1

首先，我要说的是，好问题！非常详细和可复制。我仔细研究了你的问题，试着重新做练习，从你的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")

然而，这导致了一个类似于你的阴谋：

​

​

为了检查我对盖亚档案中的坐标是否正确显示的怀疑，我从苔丝图像的中心画了一个1000弧秒的圆圈。正如您所看到的，它大致与GAIA位置的数据点云的外部(从图像中心看到)的圆形形状相一致。我只是认为这些都是GAIA DR2存档中属于您搜索区域的所有点。数据云的内部似乎有一个方形的边界，这可能来自于一个方形的视野。

Answer 2

很好的例子。请注意，您还可以使用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的命令，这将消除这种奇怪的非循环行为。此索引对于不强制初始测试的全部输出非常有用。

另外，我已经宣布了从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分钟)。完整查询将返回103574行。

​

​