More than 3 years have passed since last update.

天文学Advent Calendar 2019

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

Last updated at 2020-04-23Posted at 2019-05-11

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

これはastroqueryを使ってSkyViewというサイトから望遠鏡によって撮られた天文写真を取得してpythonへ持ち込んで用いる方法です。

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

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

SkyViewの天文写真

SkyViewには色んな掃天観測計画のデータは格納されています。その中には2006から2011年まで日本によって打ち上げられた赤外線天文衛星「あかり」も含まれています。

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

from astroquery.skyview import SkyView
SkyView.list_surveys()

(略)
'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

まずはあかりの撮った干潟星雲(通称M8)という輝線星雲の写真を見てみます。get_images関数を使ったらすぐ写真がダウンロードされます。

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)

[[<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つの波長帯の写真を取得します。

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

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()

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

zscale

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

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

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()

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

欠落の部分の補間

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

astropyには畳み込みの方法で欠落の部分を補間するinterpolate_replace_nansがあるので、それを使ったら簡単です。

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()

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

RGBの写真

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

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()

座標を表示する

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

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()

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

あかりによる天体の写真

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

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 いて座の「干潟星雲」

M16 へび座の「わし星雲」

M17 いて座の「オメガ星雲」

M20 いて座の「三裂星雲」

M31 アンドロメダ座の「アンドロメダ銀河」

M33 「さんかく座銀河」

M42 オリオン座の「オリオン大星雲」

M78 オリオン座にある星雲であり、ウルトラマンの故郷

2MASSによる天体の写真

2MASSはあかりと同じく赤外線領域で撮影する計画です。波長帯は3つ

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

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

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

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()