首页 > 编程语言> > 如何使用Python将Gaia天体测量数据绘制成TESS图像？

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

2019-10-05 05:56:07 作者：互联网

长话短说：我想将盖亚天体测量数据绘制成Python中的TESS图像.这怎么可能？
请参阅下面的详细说明.

我有一个64×64像素TESS星的图像与Gaia ID 4687500098271761792.TESS天文台指南的第8页说1像素是~21弧秒.使用Gaia Archive,我搜索这颗恒星(在顶部特征下面,点击搜索.)并提交查询以查看1000弧秒内的星星,大致是我们需要的半径.我用于搜索的名称是Gaia DR2 4687500098271761792,如下所示：

提交查询,我得到一个包含RA和DEC坐标的500颗星的列表.选择CSV和下载结果,我得到了大约4687500098271761792的星星列表.这个结果文件也可以找到here.这是我们想要使用的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)

得到一个很好的照片：

还有一堆警告,其中大部分都是here (warnings in the Q, explanation in the comments).

请注意,我们使用WCS投影确实很好.为了检查,我们只是在hdul.data中绘制数据而不关心投影：

plt.imshow(hdul.data)

结果：

与以前几乎相同,但现在轴的标签只是像素数,而不是RA and 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来自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笔记本电脑是here.

我错过了哪些步骤,或者我做错了什么？

进一步发展

在response到另一个question,keflavich友好地建议使用plotting in world coordintes的变换参数.尝试使用一些示例点(下图中的弯曲十字).在没有处理它的情节下绘制了盖亚数据,他们最终集中在一个非常狭窄的空间.应用于它们的变换方法,得到了比以前看似非常相似的结果.代码(&也是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线(或者说,恒定的纬度和经度线)更好地理解为什么来自盖亚的数据点是错误的.将以上代码展开几行：

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

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

结果令人鼓舞：

(见笔记本here.)

解决方法:

首先我要说,好问题！非常详细和可重复.我回答了你的问题并尝试从你的git repo开始重做练习,并从GAIA档案中下载目录.

编辑

以编程方式您的代码很好(请参阅下面的OLD PART以获得略有不同的方法).丢失点的问题是,从GAIA存档下载csv文件时,只能获得500个数据点.因此,看起来查询中的所有点都被塞进了奇怪的形状.但是,如果将搜索半径限制为较小的值,则可以看到TESS图像中存在点：

请与下面OLD PART中显示的版本进行比较.代码与下面的代码相同,只有下载的csv文件的半径更小.因此,在导出到csv时,您似乎只是从GAIA存档中下载了所有可用数据的一部分.绕过这种情况的方法是像你一样进行搜索.然后,在结果页面上单击底部的以ADQL形式显示查询,在查询中以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弧秒的圆圈.正如您所看到的,它大致与GAIA位置数据点云的外部圆形(从图像中心看)对齐.我只是认为这些都是GAIA DR2档案中属于您搜索范围内的所有点.数据云甚至似乎在内部有一个方形边界,这可能来自某个方形视场.

标签：python,python-3-x,astronomy,astropy,fits
来源： https://codeday.me/bug/20191005/1855201.html