[以下追記] matplotlib のバージョンが3.7.0等だと、aplpy でエラーが起きるようです。

TypeError: WCSAxes.__init__() got multiple values for argument 'wcs'

pip install matplotlib==3.5.2

などとしてダウングレードしてください。
[以上追記]

この記事では電波輝線データの積分強度図の作り方について紹介します。

今回、例として国立天文台の FUGIN プロジェクトで得られた野辺山45m電波望遠鏡の CO 輝線のアーカイブデータを用います。データは
http://jvo.nao.ac.jp/portal/nobeyama/
から
FGN_03100+0000_2x2_12CO_v1.00_cube.fits
をダウンロードしました (重いです)。

全速度積分

積分強度とは、いわばスペクトルの面積のことです。まずは前回同様 fits を読み込みます。

from astropy.io import fits
import numpy as np
from matplotlib import pyplot as plt
import aplpy

hdu = fits.open("~/your/fits/dir/FGN_03100+0000_2x2_12CO_v1.00_cube.fits")[0]

data の shape を確認します。

print(hdu.data.shape)
# (462, 848, 848)

これは、Z (速度), Y (銀緯), X (銀経) のチャンネル数を表しています。つまり、848 x 848 のピクセルがあって、その一つ一つが462 チャンネルのスペクトルを持っているということです。

では適当な場所 (424, 424) のスペクトルを見てみましょう。python では以下のように書きます。

spe = hdu.data[: ,424, 424]
plt.plot(range(len(spe)), spe)
plt.show()

ここで「:」は「全て」を意味します。また、python では 424 と書くと実際には 424+1=425 ピクセル目をとってくることに注意です。

縦軸は [K] ですが、横軸はチャンネルです。これを速度に変換します。色々方法はありますが、今回はシンプルに、0 チャンネル目 (一番最初のチャンネル) と 461 チャンネル目 (一番最後のチャンネル) の速度を計算して、線形補完することにします。そのために、astropy.wcs を使用します。

from astropy.wcs import WCS

w = WCS("~/your/fits/dir/FGN_03100+0000_2x2_12CO_v1.00_cube.fits")
print(w)
#WCS Keywords
#
#Number of WCS axes: 3
#CTYPE : 'GLON-SFL'  'GLAT-SFL'  'VRAD'  
#CRVAL : 31.0  0.0  -99675.0  
#CRPIX : 424.5  424.5  1.0  
#NAXIS : 848  848  462

チャンネル数を与えると座標が返ってくる wcs_pix2world 関数を使用します。

x_start, y_start, v_start = w.wcs_pix2world(0, 0, 0, 0) #最後の 0 は気にしなくて良い
x_end, y_end, v_end = w.wcs_pix2world(hdu.data.shape[2]-1, hdu.data.shape[1]-1, hdu.data.shape[0]-1, 0)

v_array = np.linspace(v_start, v_end, num=hdu.data.shape[0], endpoint=True)

これを使って横軸を速度にします。v_array は m/s になっているので km/s にするために1000.0で割り算します。

spe = hdu.data[: ,424, 424]
plt.plot(v_array/1000.0, spe)
plt.xlabel("Velocity [km/s]")
plt.ylabel("Tmb [K]")
plt.show()

これでスペクトルが描けました。

積分強度図とはスペクトルの面積のことでした。つまり、速度軸 (pythonでいう 0 軸目) 方向に足し合わせて速度幅をかければ完成です。

integ_map = np.nansum(hdu.data, axis=0)*hdu.header["CDELT3"]/1000.0
plt.imshow(integ_map, vmin=0, vmax=500)
plt.show()

指定速度積分

指定した速度の範囲だけを積分したい場合があると思います。上の例では np.nansum に hdu.data、つまり data 全部を渡していました。これをスライスして渡してあげれば完了です。
スライスするチャンネルを計算するために、先程とは反対の wcs_world2pix を使用します。
例として、視線速度 25 km/sから 125 km/sを積分する場合を考えます。

v_start_wcs, v_end_wcs = 25.0*1000.0, 125.0*1000.0
x_tempo, y_tempo, v_tempo = w.wcs_pix2world(0, 0, 0, 0) #最後の 0 は気にしなくて良い
x_start_ch, y_start_ch, v_start_ch = w.wcs_world2pix(x_tempo, y_tempo, v_start_wcs, 0)
x_end_ch, y_end_ch, v_end_ch = w.wcs_world2pix(x_tempo, y_tempo, v_end_wcs, 0)
v_start_ch, v_end_ch = int(round(float(v_start_ch))), int(round(float(v_end_ch)))

この v_start_ch, v_end_ch を使って data をスライスします。

integ_map = np.nansum(hdu.data[v_start_ch:v_end_ch+1], axis=0)*hdu.header["CDELT3"]/1000.0
plt.imshow(integ_map, vmin=0, vmax=500)
plt.show()

余計なチャンネルを除いたので、その分ノイズが減ってちょっと綺麗になりました。

積分強度図を fits として保存する

上でやった処理を効率化するために以下のような関数を定義します。

def v2ch(v, w): # v(km/s)をchに変える
    x_tempo, y_tempo, v_tempo   = w.wcs_pix2world(0, 0, 0, 0)
    x_ch, y_ch, v_ch   = w.wcs_world2pix(x_tempo, y_tempo, v*1000.0, 0)
    v_ch = int(round(float(v_ch), 0))
    return v_ch

def del_header_key(header, keys): # headerのkeyを消す
    import copy
    h = copy.deepcopy(header)
    for k in keys:
        try:
            del h[k]
        except:
            pass
    return h
def make_new_hdu_integ(hdu, v_start_wcs, v_end_wcs, w): # 積分強度のhduを作る
    data = hdu.data
    header = hdu.header
    start_ch, end_ch = v2ch(v_start_wcs, w), v2ch(v_end_wcs, w)
    new_data = np.nansum(data[start_ch:end_ch+1], axis=0)*np.abs(header["CDELT3"])/1000.0
    header = del_header_key(header, ["CRVAL3", "CRPIX3", "CRVAL3", "CDELT3", "CUNIT3", "CTYPE3", "CROTA3", "NAXIS3", "PC1_3", "PC2_3", "PC3_3", "PC3_1", "PC3_2"])
    header["NAXIS"] = 2
    new_hdu = fits.PrimaryHDU(new_data, header)
    return new_hdu

25 km/s から 125 km/s の積分強度の hdu を以下で作ります。

integ_hdu = make_new_hdu_integ(hdu, 25.0, 125.0, w)

前回と同様、aplpy で plot してみます。

fig = plt.figure(figsize=(8, 8))
f = aplpy.FITSFigure(integ_hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(vmin=1, vmax=500, stretch='log', cmap="jet", aspect="equal")
plt.show()

名前をつけて fits で保存するには、以下を実行します。

integ_hdu.writeto("integ_25_125.fits", overwrite=True)

以上です。

リンク
目次