如何使用Python将Gaia占星数据绘制到TESS图像上? [英] How to plot Gaia astrometry data to TESS images using Python?
问题描述
长话短说:我想用Python将Gaia天文测量数据绘制到TESS图像中.这怎么可能? 详见下文.
Long story short: I want to plot Gaia astrometry data to TESS imagery in Python. How is it possible? See below for elaborated version.
我有一个64x64像素的 TESS 一颗具有 Gaia ID 4687500098271761792 的恒星图像em>. TESS天文台指南的第8页说1像素约为21弧秒.使用 Gaia存档，我搜索了这颗星(在主要功能下方 >，点击 Search .)并提交查询，以查看1000弧秒以内的恒星，大约是我们所需的半径.我用于搜索的名称是Gaia DR2 4687500098271761792，如下所示:
I have 64x64 pixel 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, I search for this star (below top features, click Search.) and submit a query to see the stars within 1000 arcsec, roughly the radius we need. The name I use for the search is Gaia DR2 4687500098271761792, as shown below:
提交查询，我得到500个具有RA和DEC坐标的星星的列表.选择CSV和Download results，我会得到 4687500098271761792 周围的星星列表.也可以在此处找到该文件.这是我们要使用的 Gaia 的输入.
Submit Query, and I get a list of 500 stars with RA and DEC coordinates. Select CSV and Download results, I get the list of stars around 4687500098271761792. This resulting file also can be found here. This is the input from Gaia we want to use.
在TESS中，我们有 4687500098271761792_med.fits ，图像文件.我们使用以下方式绘制它:
From the TESS, we have 4687500098271761792_med.fits, an image file. We plot it using:
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)
并得到一张漂亮的照片:
and get a nice pic:
and a bunch warnings, most of which was kindly explained here (warnings in the Q, explanation in the comments).
请注意，我们使用 WCS 投影确实很好.为了进行检查，我们只需要在hdul.data中绘制数据，而不必关心投影:
Notice that it is indeed good that we are using the WCS projection. To check, let's just plot the data in hdul.data without caring about the projection:
plt.imshow(hdul.data)
结果:
与以前几乎相同，但是现在轴的标签只是像素编号，而不是 RA和DEC ，这将是更好的选择.第一个图中的DEC和RA值分别约为-72°和16°，这是很好的，因为Gaia目录为我们提供了 4687500098271761792 附近的恒星坐标.因此，预测似乎还可以.
Almost the same as before, but now the labels of the axes are just pixel numbers, not RA and DEC, as would be preferable. The DEC and RA values in the first plot are around -72° and 16° respectively, which is good, given that the Gaia catalogue gave us stars in the proximity of 4687500098271761792 with roughly these coordinates. So the projection seems to be reasonably ok.
现在让我们尝试在imshow()图上方绘制盖亚星.我们读入我们先前下载的CSV文件，并从中提取对象的RA和DEC值:
Now let's try to plot the Gaia stars above the imshow() plot. We read in the CSV file we downloaded earlier and extract the RA and DEC values of the objects from it:
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='+')
形状不是预期的圆.这可能预示着将来的麻烦.
Shape is not a circle as expected. This might be an indicator of future troubles.
让我们尝试将这些RA和DEC值转换为WCS，并以这种方式绘制它们:
Lets try to transform these RA and DEC values to WCS, and plot them that way:
for index, each in enumerate(ralist):
ra, dec = wcs.all_world2pix([each], [declist[index]], 1)
plt.scatter(ra, dec, marker='+', c='k')
结果是:
The function all_world2pix is from here. The 1 parameter just sets where the origin is. The description of all_world2pix says:
在这里，原点是图像左上角的坐标. 在FITS和Fortran标准中，为1.在Numpy和C标准中. 这是0.
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数据放在一起:
Nevertheless, the shape of the point distribution we get is not promising at all. Let's put together the TESS and Gaia data:
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笔记本可在此处使用.
which is not anywhere near what would be ideal. I expect to have an underlying imshow() pic with many markers on it, and the markers should be where stars are on the TESS image. The Jupyter Notebook I worked in is available here.
我错过了哪些步骤，或者我做错了什么?
进一步的发展
在响应中，对另一个 keflavich 建议将transform参数用于此处):
In response to another question, 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()绘制)平行)，因此可以预见到该弯曲十字.现在让我们尝试绘制常数RA& DEC线(或者说恒定的纬度和经度线)以更好地理解为什么Gaia的数据点放错了位置.将上面的代码扩展几行:
This bending cross is expected, since TESS is not aligned with constant latitude and longitude lines (ie the arms of the cross do not have to be parallel with the sides of the TESS image, plotted with imshow() ). Now lets try to plot constant RA & DEC lines (or to say, constant latitude and longitude lines) to better understand why the datapoints from Gaia are misplaced. Expand the code above by a few lines:
ax.coords.grid(True, color='green', ls='solid')
overlay = ax.get_coords_overlay('icrs')
overlay.grid(color='red', ls='dotted')
结果令人鼓舞:
(请在笔记本此处.)
推荐答案
首先，我要说的是一个很好的问题！非常详细且可重现.我经历了您的问题，并尝试从您的git repo开始重做练习，并从GAIA存档中下载目录.
First I have to say, great question! Very detailed and reproducible. I went through your question and tried to redo the exercise starting from your git repo and downloading the catalogue from the GAIA archive.
通过编程，您的代码可以正常运行(有关其他方法，请参见下面的 OLD PART ).缺少点的问题是，从GAIA档案库下载csv文件时，一个人只能获得500个数据点.因此，看起来查询中的所有点都塞满了怪异的形状.但是，如果将搜索半径限制为较小的值，则可以看到TESS图像中有一些点:
Programmatically your code is fine (see OLD PART below for a slightly different approach). The problem with the missing points is that one only gets 500 data points when downloading the csv file from the GAIA archive. Therefore it looks as if all points from the query are crammed into a weird shape. However if you restrict the radius of the search to a smaller value you can see that there are points that lie within the TESS image:
请与以下旧版本中显示的版本进行比较.该代码与以下代码相同，只是下载的csv文件半径较小.因此，当您导出到csv时，您似乎只是从GAIA存档中下载了所有可用数据的一部分.避免这种情况的方法是像您一样进行搜索.然后，在结果页面上，单击底部的Show query in ADQL form，然后在查询中以SQL格式更改显示:
please compare to the version shown below in the OLD PART. The code is the same as below only the downloaded csv file is for a smaller radius. Therefore it seems that you just downloaded a part of all available data from the GAIA archive when exporting to csv. The way to circumvent this is to do the search as you did. Then, on the result page click on Show query in ADQL form on the bottom and in the query you get displayed in SQL format change:
Select Top 500
到
Select
在查询的开头.
对于绘图，我使用了aplpy-在后台使用了matplotlib-并得到了以下代码:
For plotting I used aplpy - uses matplotlib in the background - and ended up with the following code:
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")
但是，这会导致与您的情节相似的情节:
However this results in a plot similar to yours:
要检查我是否怀疑GAIA档案中的坐标正确显示，我从TESS图像的中心绘制了一个1000 arcsec的圆.如您所见，它与GAIA位置的数据点云的外侧(从图像中心看)的圆形大致对齐.我只是认为这些都是GAIA DR2归档文件中属于您搜索区域的所有要点.数据云甚至似乎在内部具有方形的边界，它可能来自方形视场.
To check my suspicion that the coordinates from the GAIA archive are correctly displayed I draw a circle of 1000 arcsec from the center of the TESS image. As you can see it aligns roughly with the circular shape of the outer (seen from the center of the image) side of the data point cloud of the GAIA positions. I simply think that these are all points in the GAIA DR2 archive that fall within the region you searched. The data cloud even seems to have a squarish boundary on the inside, which might come from something as a square field of view.
这篇关于如何使用Python将Gaia占星数据绘制到TESS图像上?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
