天文データ解析入門 その23 (GAIA DR3 の扱い方)

本記事では、[GAIAのDR3](https://www.cosmos.esa.int/web/gaia/dr3) の、python を用いた基本的な扱い方について記述します。

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

まずは必要なものを import します。
```python
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 のデータを見ます。
```python
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()
```
![ダウンロード (17).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/3f49e771-007f-f8ae-7a4e-bbf654d03733.png)

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

まずは以下の変数を書き換えます。
```python
Gaia.MAIN_GAIA_TABLE = "gaiadr3.gaia_source"  # DR3
# Gaia.MAIN_GAIA_TABLE = "gaiaedr3.gaia_source"  # EDR3
Gaia.ROW_LIMIT = -1 # -1を指定すると上限青天井
```

今回は例として、適当にこのあたりの星を取得します。
```python
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()
```
![ダウンロード (21).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/17a91015-be61-185e-617b-714d6c0e6552.png)


半径 0.5 度の範囲を取得します。
```python
radius = 0.5*u.deg
j = Gaia.cone_search_async(c, radius)
r = j.get_results()
```

少し時間がかかります。

結果を見ます。
```python
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 ...>
```

もしくは
```python
r[0]
```
などで見られます。

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

```python
df = r.to_pandas()
```
これでpandas 形式になりました。
csv で保存したいときは
```python
df.to_csv("stars_test.csv")
```
とします。

# 使い方の例

例えば、"radial_velocity" のデータが入っているものだけ欲しい場合は、
```python
df_vr = df[df["radial_velocity"]==df["radial_velocity"]]
```
などで絞ります。(※ バージョンによってキーの名前が微妙に異なっている可能性があります。)

"radial_velocity"が 10 km/s 以上のものだけ拾いたい場合は
```python
df[df["radial_velocity"]>10.0]
```
などとします。

星までの距離を"parallax"から計算します。
```python
distance_all = 1.0/(np.abs(df["parallax"])/1000.0)
```
中身を見ます。
```python
plt.hist(distance_all, bins=np.linspace(0, 3000.0, 100))
plt.xlabel("pc")
plt.ylabel("Number of stars")
```
![ダウンロード (19).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/4168ad3f-4ec4-5cf4-8b3a-6563811a3c45.png)

```python
plt.hist(df[df["radial_velocity"]==df["radial_velocity"]]["radial_velocity"], bins=50)
plt.xlabel("km/s")
plt.ylabel("Number of stars")
```
![ダウンロード (20).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/2934c0ea-77bf-bcd9-3092-8f310fdeb8c0.png)

マップの上に plot します。
```python
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()
```
![ダウンロード (22).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/5c38dd2d-f08c-ab8b-f232-84d5642fee6c.png)

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

どの距離から減光されて居なくなっているのかを調べます。
```python
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")
```
![ダウンロード (23).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/00e1afe2-947c-3baa-2455-dd1b3fa438e3.png)

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

Herschel の強度も使ってみます。
```python
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 します。

```python
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()
```
![ダウンロード (24).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/27aaefac-d28f-5818-bf32-d52624adf25e.png)

```python
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")
```
![ダウンロード (26).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/af98cfef-8824-a491-8469-83837eb1385f.png)

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

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

以上です。

リンク
[目次](https://qiita.com/Shinji_Fujita/items/ecb869e5b6a9fa468483)