Go to Qiita Advent Calendar 2024 Top
LoginSignup
89
67

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

天文学Advent Calendar 2019

Day 3
@phyblas(雛星 φυβλας)

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

Last updated at Posted 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()

q01.png

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

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

q02.png

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

欠落の部分の補間

さっきの結果を見てあかりの写真には一部のデータが欠けているとわかります。その部分は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()

q03.png

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

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

q04.png

座標を表示する

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

q05.png

天の川の線は銀緯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 いて座の「干潟星雲」

akari_M8.jpg

M16 へび座の「わし星雲」

akari_M16.jpg

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

akari_M17.jpg

M20 いて座の「三裂星雲」

akari_M20.jpg

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

akari_M31.jpg

M33 「さんかく座銀河」

akari_M33.jpg

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

akari_M42.jpg

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

akari_M78.jpg

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

M17

2mass_M17.jpg

M42のすぐ北の方にある星雲であるM43。右下でもっと明るいのはM42

2mass_M43.jpg

M51 りょうけん座の子持ち銀河

2mass_M51.jpg

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

2mass_M78.jpg

M83 うみへび座にある「南の回転花火銀河」

2mass_M83.jpg

M104 おとめ座にある「ソンブレロ銀河」

2mass_M104.jpg

M109 おおぐま座にある棒渦巻銀河

2mass_M109.jpg

その他に

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

UKIDSS

UKIDSS-K, UKIDSS-H, UKIDSS-J

UKIDSSはイギリス赤外線望遠鏡(UKIRT)による掃天観測計画

M2 みずがめ座にある球状星団

UKIDSS_M2.jpg

Mellinger

Mellinger Red, Mellinger Green, Mellinger Blue

Mellingerによる可視光の波長帯の写真

M8

Mellinger008.jpg

M45 おうし座の昴星団

Mellinger045.jpg

WISE

WISE 12, WISE 4.6, WISE 3.4

広域赤外線探査衛星（WISE）による写真

M27 こぎつね座にある亜鈴状星雲

wise027.jpg

M51

wise051.jpg

M81 おおぐま座にある渦巻銀河

wise081.jpg

M96 しし座にある渦巻銀河

wise096.jpg

M101 おおぐま座にある回転花火銀河

wise101.jpg

終わりに

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

89
67
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
89
67

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?