赤外線天文衛星「あかり」などの望遠鏡によって撮影された写真をいじってみたり

最近[astropy](http://docs.astropy.org)と[astroquery](https://astroquery.readthedocs.io)の色んな機能を試しています。そして先日@nishimuraatsushiさんのこの記事を見つけたのです。 https://qiita.com/nishimuraatsushi/items/5bf9c869813133352c56

これはastroqueryを使って[SkyView](https://skyview.gsfc.nasa.gov/current/cgi/query.pl)というサイトから望遠鏡によって撮られた天文写真を取得してpythonへ持ち込んで用いる方法です。

この記事から色々勉強しました。面白いので、コードを参考にして自分も色々試しました。

この記事ではSkyViewの天文写真をいじってみたりします。



# SkyViewの天文写真

SkyViewには色んな[掃天観測](https://ja.wikipedia.org/wiki/掃天観測)計画のデータは格納されています。その中には2006から2011年まで日本によって打ち上げられた[赤外線天文衛星「あかり」](https://ja.wikipedia.org/wiki/あかり_(人工衛星))も含まれています。

astroqueryモジュールのskyviewを使ったらすぐSkyViewの写真を取得できます。list_surveysという関数を使ったらすべての計画の名前が得られます。

```python
from astroquery.skyview import SkyView
SkyView.list_surveys()
```

```javascript
(略)
'IR:AKARI': ['AKARI N60', 'AKARI WIDE-S', 'AKARI WIDE-L', 'AKARI N160'],
(略)
```

あかりには遠赤外線の4つの波長帯のフィルターがあります。

- AKARI N60: 50–80 µm
- AKARI WIDE-S: 60–110 µm
- AKARI WIDE-L: 110–180 µm
- AKARI N160: 140–180 µm

まずはあかりの撮った[干潟星雲](https://ja.wikipedia.org/wiki/干潟星雲)(通称M8)という[輝線星雲](https://ja.wikipedia.org/wiki/輝線星雲)の写真を見てみます。[get_images](https://astroquery.readthedocs.io/en/latest/api/astroquery.skyview.SkyViewClass.html#astroquery.skyview.SkyViewClass.get_images)関数を使ったらすぐ写真がダウンロードされます。

```python
from astropy import units as u

position = 'M8' # 欲しいところの位置か天体の名前
survey = ['AKARI N160','AKARI WIDE-L','AKARI WIDE-S','AKARI N60']
radius = 5*u.deg # 範囲
pixels = 800 # 欲しい画像の大きさ
hdu = SkyView.get_images(position=position,
                         survey=survey,
                         radius=radius,
                         pixels=pixels)
print(hdu)
print(hdu[3][0].data.shape)
```

```javascript
[[<astropy.io.fits.hdu.image.PrimaryHDU object at 0x12895efd0>], [<astropy.io.fits.hdu.image.PrimaryHDU object at 0x1286d1c18>], [<astropy.io.fits.hdu.image.PrimaryHDU object at 0x1286bb4e0>], [<astropy.io.fits.hdu.image.PrimaryHDU object at 0x1286afd68>]]
(800, 800)
```

positionは取りたい部分の位置を入れてもいいが、天体の名前を入れたら自動的にその天体の位置はsimbadから取られるので、直接位置を入れなくてもよくて便利です。

surveyにはリストを入れて一度でたくさんな波長帯のデータが得られます。ここではあかりの4つの波長帯の写真を取得します。

それぞれの写真を表示してみます。

```python
import matplotlib.pyplot as plt

plt.figure(figsize=[7,7])
for i in range(4):
    plt.subplot(221+i,title=survey[i])
    plt.imshow(hdu[i][0].data,origin='lower',cmap='hot')
plt.show()
```

![q01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/805b8e8c-2c3d-a3eb-82c5-cfa7bb147041.png)

あまりよく見えないのですね。天文写真は大体こんな感じで明るい部分と暗い部分の差が大きすぎて尺度をちゃんと調整しないと大事な部分ははっきり見えないのです。



# zscale

大事な部分を強調するために天文写真でよく使われる方法は「IRAF zscale」です。詳しくはこちらhttp://docs.astropy.org/en/stable/api/astropy.visualization.ZScaleInterval.html

zscaleに変換する関数はastropyに含まれているので、ここではこれを使います。

```python
from astropy.visualization import ZScaleInterval,ImageNormalize

plt.figure(figsize=[7,7])
for i in range(4):
    plt.subplot(221+i,title=survey[i])
    norm = ImageNormalize(hdu[i][0].data,interval=ZScaleInterval())
    plt.imshow(hdu[i][0].data,norm=norm,origin='lower',cmap='hot')
plt.show()
```

![q02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/bd51a4db-1d6a-c8d9-cd80-e31fb43e00a3.png)

はっきり見えるようになったようです。中心の明るい部分は干潟星雲で、右側の広い線は天の川の中心です。いて座は天の川銀河の中心の方向にあるから、干潟星雲などの有名な天体がたくさん集まっています。



# 欠落の部分の補間

さっきの結果を見てあかりの写真には一部のデータが欠けているとわかります。その部分は0になっています。

astropyには畳み込みの方法で欠落の部分を補間する[interpolate_replace_nans](http://docs.astropy.org/en/stable/api/astropy.convolution.interpolate_replace_nans.html)があるので、それを使ったら簡単です。

```python
from astropy.convolution import Gaussian2DKernel, interpolate_replace_nans
import numpy as np

akari = np.array([hdu[i][0].data for i in range(4)])
kernel = Gaussian2DKernel(5,5)
for i in range(4):
    akari[i][akari[i]<=0] = np.nan
    akari[i] = interpolate_replace_nans(akari[i],kernel)

plt.figure(figsize=[7,7])
for i in range(4):
    plt.subplot(221+i,title=survey[i])
    norm = ImageNormalize(akari[i],interval=ZScaleInterval())
    plt.imshow(akari[i],norm=norm,origin='lower',cmap='hot')
plt.show()
```

![q03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/cf8b4799-7109-4ecb-e119-4be0a67fac1d.png)

これでもともと暗かった部分は埋まって綺麗になりました。



# RGBの写真

次は3つの波長帯を合成してRGBの写真にします。全部4つあるので、4つの組む方法があります。通常順番は決まって一番波長の長いのは赤に、短いのが青になります。

```python
for i in range(4):
    akari[i] = ImageNormalize(akari[i],interval=ZScaleInterval())(akari[i])

plt.figure(figsize=[7,7])
for i in range(4):
    ii = [0,1,2,3]
    ii.remove(i)
    img = np.stack([akari[a] for a in ii],2)
    plt.subplot(221+i,title=r' | '.join(survey[a].split(' ')[1] for a in ii))
    plt.imshow(img,origin='lower')
plt.show()
```

![q04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/4e05ead4-eda7-7b93-ae19-9ae410a538f6.png)



# 座標を表示する

SkyViewの写真では主に赤道座標が使われています。座標のデータはfitsのhduのヘッダーから取得してmatplotlibで描く時に表示できます。同時に[銀河座標](https://ja.wikipedia.org/wiki/銀河座標)も一緒に入れて表示できます。

```python
from astropy.wcs import WCS

plt.figure(figsize=[7,7])
wcs = WCS(hdu[0][0].header) # wcs座標をヘッダーから取得
ax = plt.gca(projection=wcs) # wcsをプロットに使う
img = np.stack([akari[0],akari[2],akari[3]],2) # N160とWIDE-SとN60を3つ色に使う
plt.imshow(img)
plt.xlabel('赤經',fontname='AppleGothic')
plt.ylabel('赤緯',fontname='AppleGothic')
galactic = ax.get_coords_overlay('galactic') # 銀河座標
galactic[0].set_ticks(spacing=2*u.deg) # 銀經を描く
galactic[1].set_ticks(spacing=2*u.deg) # 銀緯を描く
galactic.grid(color='w',linestyle='--')
galactic[0].set_axislabel('銀經',fontname='AppleGothic')
galactic[1].set_axislabel('銀緯',fontname='AppleGothic')
plt.show()
```

![q05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/82ea220d-d9f7-711f-d1b1-8e155bdd988c.png)

天の川の線は銀緯0度のところに架かっています。銀経と銀緯は0度であることろは天の川の中心に向かう方向なので、この部分は天の川の中心に随分近いとわかります。



# あかりによる天体の写真

次は他に色んな天体を見ます。残念ながら、SkyViewでのあかりの写真の解像度があまり高くないようだからちゃんと見える天体は多くないです。ここではその一部を挙げてみます。

```python
import numpy as np
import matplotlib.pyplot as plt
from astroquery.skyview import SkyView
from astropy import units as u
from astropy.visualization import ZScaleInterval,ImageNormalize
from astropy.convolution import Gaussian2DKernel, interpolate_replace_nans
from astropy.wcs import WCS

for m in [8,16,17,20,31,33,42,78]:
    survey = ['AKARI N160','AKARI WIDE-S','AKARI N60']
    hdu = SkyView.get_images('M%d'%m,survey=survey,radius=1*u.deg,pixels=640)
    img = np.stack([hdu[0][0].data,hdu[1][0].data,hdu[2][0].data],2)
    wcs = WCS(hdu[0][0].header)
    
    # 畳み込みで補間
    img[img==0] = np.nan
    akari = np.empty_like(img)*np.nan
    for k in range(1,21):
        if(np.any(np.isnan(akari))):
            kernel = Gaussian2DKernel(k,k)
            for i in range(3):
                akari[:,:,i] = interpolate_replace_nans(img[:,:,i],kernel)
        else:
            break
    
    # zscaleへ変換
    for i in range(3):
        akari[:,:,i] = ImageNormalize(akari[:,:,i],interval=ZScaleInterval())(akari[:,:,i])
    
    plt.figure(figsize=[5,5])
    plt.gca(projection=wcs,title='M%d'%m)
    plt.imshow(akari)
    plt.savefig('akari_M%d.jpg'%m)
    plt.close()
```

M8 [いて座](https://ja.wikipedia.org/wiki/いて座)の「干潟星雲」

![akari_M8.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/689671db-d0b6-4792-6c80-91f55ac55521.jpeg)

M16 [へび座](https://ja.wikipedia.org/wiki/へび座)の「[わし星雲](https://ja.wikipedia.org/wiki/わし星雲)」

![akari_M16.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/6d7b728a-d1ea-be6b-69e2-c1d34de0e4d9.jpeg)

M17 いて座の「[オメガ星雲](https://ja.wikipedia.org/wiki/オメガ星雲)」

![akari_M17.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/ef447e9a-361c-1c2d-aca6-15aa21a311cc.jpeg)

M20 いて座の「[三裂星雲](https://ja.wikipedia.org/wiki/三裂星雲)」

![akari_M20.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/8ebb7029-4403-860b-1c88-28ed52eb2dcb.jpeg)

M31 [アンドロメダ座](https://ja.wikipedia.org/wiki/アンドロメダ座)の「[アンドロメダ銀河](https://ja.wikipedia.org/wiki/アンドロメダ銀河)」

![akari_M31.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/a9b1de40-e56b-ceb0-6154-1f8d2de6c56e.jpeg)

M33 「[さんかく座銀河](https://ja.wikipedia.org/wiki/さんかく座銀河)」

![akari_M33.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/65bd1d7b-d99a-6541-d38f-68b02ceb89cf.jpeg)

M42 [オリオン座](https://ja.wikipedia.org/wiki/オリオン座)の「[オリオン大星雲](https://ja.wikipedia.org/wiki/オリオン大星雲)」

![akari_M42.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/53c4baf1-9359-7039-1840-09f613b40dbb.jpeg)

[M78](https://ja.wikipedia.org/wiki/M78) オリオン座にある星雲であり、ウルトラマンの故郷

![akari_M78.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/b681ade6-2ddb-43e7-f05e-d8e9ff658c4c.jpeg)



# 2MASSによる天体の写真

[2MASS](https://ja.wikipedia.org/wiki/2MASS)はあかりと同じく赤外線領域で撮影する計画です。波長帯は3つ

- 2MASS-J 1.25 μ
- 2MASS-H 1.65 μ
- 2MASS-K 2.17 μ

SkyViewではあかりより2MASSの写真の方が解像度が高いのでもっと性質高い写真ができます。

2MASSの写真を使ってみます。

```python
survey = ['2MASS-K','2MASS-H','2MASS-J']
for m in [17,43,51,78,83,104,109]:
    hdu = SkyView.get_images('M%d'%m,survey=survey,radius=20*u.arcmin,pixels=1280)    
    plt.figure(figsize=[6,6]).gca(projection=WCS(hdu[0][0].header),title='M%d'%m)
    plt.imshow(np.stack([ImageNormalize(hdu[i][0].data,interval=ZScaleInterval())(hdu[i][0].data) for i in range(3)],2))
    plt.savefig('2mass_M%d.jpg'%m)
    plt.close()
```

M17

![2mass_M17.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/5afdee13-c7c3-8ce0-c592-e5261e4a5c03.jpeg)

M42のすぐ北の方にある星雲である[M43](https://ja.wikipedia.org/wiki/M43)。右下でもっと明るいのはM42

![2mass_M43.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/a0b1214a-b51d-d1a4-2283-bcff4ca01b3b.jpeg)

M51 [りょうけん座](https://ja.wikipedia.org/wiki/りょうけん座)の[子持ち銀河](https://ja.wikipedia.org/wiki/子持ち銀河)

![2mass_M51.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/dd042537-4cc6-9b41-229d-c6de392ecbe1.jpeg)

ウルトラマンの故郷であるM78

![2mass_M78.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/5b308bbf-2703-9f6b-874e-8d0bef3a50bb.jpeg)

M83 [うみへび座](https://ja.wikipedia.org/wiki/うみへび座)にある「[南の回転花火銀河](https://ja.wikipedia.org/wiki/M83_(天体))」

![2mass_M83.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/b3a05558-7a4d-3fb8-cf35-bc344181cb12.jpeg)

M104 [おとめ座](https://ja.wikipedia.org/wiki/おとめ座)にある「[ソンブレロ銀河](https://ja.wikipedia.org/wiki/ソンブレロ銀河)」

![2mass_M104.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/525d8d0d-5862-83f9-8985-713a245741fb.jpeg)

[M109](https://ja.wikipedia.org/wiki/M109_(天体)) [おおぐま座](https://ja.wikipedia.org/wiki/おおぐま座)にある[棒渦巻銀河](https://ja.wikipedia.org/wiki/棒渦巻銀河)

![2mass_M109.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/58de86ea-b8a0-c112-9bdb-13481de2673c.jpeg)



# その他に

SkyViewには他にも色んな望遠鏡や計画の写真があります。ここでは一部の写真を挙げてみます。コードは大体2massの写真を描くのと同じなので省略します。

### UKIDSS

**UKIDSS-K, UKIDSS-H, UKIDSS-J**

UKIDSSは[イギリス赤外線望遠鏡](https://ja.wikipedia.org/wiki/イギリス赤外線望遠鏡)(UKIRT)による掃天観測計画

[M2](https://ja.wikipedia.org/wiki/M2_(天体)) [みずがめ座](https://ja.wikipedia.org/wiki/みずがめ座)にある[球状星団](https://ja.wikipedia.org/wiki/球状星団)

![UKIDSS_M2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/a06dcb69-7492-8d8d-0f40-0023692e3218.jpeg)

### Mellinger

**Mellinger Red, Mellinger Green, Mellinger Blue**

[Mellinger](https://skyview.gsfc.nasa.gov/blog/index.php/2009/12/16/mellinger-optical-survey)による可視光の波長帯の写真

M8

![Mellinger008.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/86bbe5cc-0be4-47b6-ed65-1e0c3f445913.jpeg)

M45 [おうし座](https://ja.wikipedia.org/wiki/おうし座)の[昴星団](https://ja.wikipedia.org/wiki/プレアデス星団)

![Mellinger045.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/f8332484-44d1-e4ec-baf3-6b28387d1e04.jpeg)

### WISE

**WISE 12, WISE 4.6, WISE 3.4**

[広域赤外線探査衛星](https://ja.wikipedia.org/wiki/広域赤外線探査衛星)（WISE）による写真

M27 [こぎつね座](https://ja.wikipedia.org/wiki/こぎつね座)にある[亜鈴状星雲](https://ja.wikipedia.org/wiki/亜鈴状星雲)

![wise027.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/7d8d4cb6-f148-b30d-c5d3-6f6d204428a2.jpeg)

M51

![wise051.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/2e4c8b53-46cb-f12b-e172-66ddcdabf88f.jpeg)

[M81](https://ja.wikipedia.org/wiki/M81_(天体)) おおぐま座にある渦巻銀河

![wise081.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/53c2396a-dfd8-eb7a-5272-6b08b8075fd2.jpeg)

[M96](https://ja.wikipedia.org/wiki/M96_(天体)) [しし座](https://ja.wikipedia.org/wiki/しし座)にある渦巻銀河

![wise096.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/17895a55-4c36-e482-efc6-c95b99821214.jpeg)

M101 おおぐま座にある[回転花火銀河](https://ja.wikipedia.org/wiki/M101_(天体))

![wise101.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/202903/9e0d3426-2fa1-8fd8-4053-d4f1f01f9f07.jpeg)



# 終わりに

こうやってSkyViewの色んな天体の写真を取得して3色合成図を作ることができました。こんな簡単なようにできるようにするastropyとastroqueryは本当に天文学者にとってすごく便利なpythonモジュールです。