本記事では、GAIAのDR3 の、python を用いた基本的な扱い方について記述します。

今回は、比較の例としてHerschelのアーカイブデータも用います。データは
http://herschel.esac.esa.int/Science_Archive.shtml
から
Taurus
と検索して、一番上に出てくる level3 のfits
hspireplw_30pxmp_0436_p2515_1476896810029.fits
を使用します。

まずは必要なものを import します。

from astropy.io import fits
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import aplpy
from astropy import units as u
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
from astroquery.gaia import Gaia
import pandas as pd
pd.set_option('display.max_columns', 100)

とりあえず、Herschel のデータを見ます。

hdu = fits.open("hspireplw_30pxmp_0436_p2515_1476896810029.fits")[1] ### 普通は[0]
w = WCS(hdu)

fig = plt.figure(figsize=(8, 8))
f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(cmap="jet", aspect="equal")
plt.show()

Gaia データのダウンロード

Gaiaのデータは膨大ですので、全部ダウンロードすると大変なことになります。基本的には必要な部分のみをダウンロードして使います。
ウェブのフォームからもダウンロード可能なのですが、一度に取得できる量に制限があるので、python による方法をお勧めします。

まずは以下の変数を書き換えます。

Gaia.MAIN_GAIA_TABLE = "gaiadr3.gaia_source"  # DR3
# Gaia.MAIN_GAIA_TABLE = "gaiaedr3.gaia_source"  # EDR3
Gaia.ROW_LIMIT = -1 # -1を指定すると上限青天井

今回は例として、適当にこのあたりの星を取得します。

coords = "4:19:52.0 +27:12:30.0"
c = SkyCoord(coords, frame="fk5", unit=(u.hourangle, u.deg))

fig = plt.figure(figsize=(8, 8))
f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(cmap="Greys", aspect="equal")
f.show_markers(c.fk5.ra.deg, c.fk5.dec.deg)
f.recenter(c.fk5.ra.deg, c.fk5.dec.deg, width=0.5, height=0.5)
plt.show()

半径 0.5 度の範囲を取得します。

radius = 0.5*u.deg
j = Gaia.cone_search_async(c, radius)
r = j.get_results()

少し時間がかかります。

結果を見ます。

r.pprint

#<bound method Table.pprint of <Table length=8131>
#        dist         ...
#                     ...
#      float64        ...
#-------------------- ...
# 0.01910917798307798 ...
#0.021585256432401608 ...
# 0.02757292230391897 ...
#0.029302247109999367 ...
#0.031447477553026355 ...
#0.032318879498195584 ...
# 0.03500904509716016 ...
# 0.03884074378865377 ...
#0.039060427046933945 ...
#                 ... ...
#  0.4995100823717482 ...
# 0.49954443834279916 ...
#  0.4996255278820678 ...
#  0.4996465140393764 ...
# 0.49977119463277786 ...
#  0.4997756139360711 ...
# 0.49977785016602705 ...
# 0.49991212461028306 ...
# 0.49994741414110117 ...
#  0.4999627788926277 ...>

もしくは

r[0]

などで見られます。

このままでは少し扱いづらいので、pandas 形式に変換します。

df = r.to_pandas()

これでpandas 形式になりました。
csv で保存したいときは

df.to_csv("stars_test.csv")

とします。

使い方の例

例えば、"radial_velocity" のデータが入っているものだけ欲しい場合は、

df_vr = df[df["radial_velocity"]==df["radial_velocity"]]

などで絞ります。(※ バージョンによってキーの名前が微妙に異なっている可能性があります。)

"radial_velocity"が 10 km/s 以上のものだけ拾いたい場合は

df[df["radial_velocity"]>10.0]

などとします。

星までの距離を"parallax"から計算します。

distance_all = 1.0/(np.abs(df["parallax"])/1000.0)

中身を見ます。

plt.hist(distance_all, bins=np.linspace(0, 3000.0, 100))
plt.xlabel("pc")
plt.ylabel("Number of stars")

plt.hist(df[df["radial_velocity"]==df["radial_velocity"]]["radial_velocity"], bins=50)
plt.xlabel("km/s")
plt.ylabel("Number of stars")

マップの上に plot します。

fig = plt.figure(figsize=(12, 12))
f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(cmap="Greys", aspect="equal")
f.show_markers(df["ra"], df["dec"], c="r", s=5)
f.recenter(c.fk5.ra.deg, c.fk5.dec.deg, width=1.0, height=1.0)
plt.show()

綺麗に減光されていることがわかります。

どの距離から減光されて居なくなっているのかを調べます。

dist_list = np.linspace(0, 500, 25)
print(dist_list)
#[  0.          20.83333333  41.66666667  62.5         83.33333333
# 104.16666667 125.         145.83333333 166.66666667 187.5
# 208.33333333 229.16666667 250.         270.83333333 291.66666667
# 312.5        333.33333333 354.16666667 375.         395.83333333
# 416.66666667 437.5        458.33333333 479.16666667 500.        ]

num_stars_list = [len(distance_all[distance_all>dist]) for dist in dist_list]

plt.plot(dist_list, num_stars_list)
plt.xlabel("pc")
plt.ylabel("Number of stars")

思ったよりわかりませんでした。

Herschel の強度も使ってみます。

ra_ch_list, dec_ch_list = w.wcs_world2pix(df["ra"], df["dec"], 0)
ra_ch_list = [int(round(ra_ch)) for ra_ch in ra_ch_list]
dec_ch_list = [int(round(dec_ch)) for dec_ch in dec_ch_list]

int_array = np.array([hdu.data[dec_ch, ra_ch] for  ra_ch, dec_ch in zip(ra_ch_list, dec_ch_list)])

例えば、Herschel の強度が 20 MJy/sr 以上のものだけ plot します。

mask = int_array>=20.0

fig = plt.figure(figsize=(12, 12))
f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(cmap="Greys", aspect="equal")
f.show_markers(df["ra"][mask], df["dec"][mask], c="r", s=5)
f.recenter(c.fk5.ra.deg, c.fk5.dec.deg, width=1.0, height=1.0)
plt.show()

num_stars_list_over20 = [len(distance_all[(distance_all>dist) & (int_array>=20.0)]) for dist in dist_list]

plt.plot(dist_list, num_stars_list_over20)
plt.xlabel("pc")
plt.ylabel("Number of stars")

距離 150 pc 付近に大きな段差があるような気もしますし、ないような気もします。

といった感じで、pandas 形式にしてしまえば、あとは自由な解析が可能です。

以上です。

リンク
目次