天文データ解析入門 その5 (位置–速度図の作り方)

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

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

pip install matplotlib==3.5.2

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

この記事では位置–速度図の作り方を紹介します。

前回 は積分強度図の作り方を紹介しました。これは ３D データ (cube) を速度方向に足し合わせていました。今回はそれを座標軸方向 (空間 grid 方向) に足し合わせたものになります。
(斜め方向の位置－速度図に関してはまた別の記事で紹介します。)

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

全座標積分

まずは例によって必要なものを import し、fitsを読み込みます。

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

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

data の shape を確認します。

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

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

yv_map = np.nansum(hdu.data, axis=2)*np.abs(hdu.header["CDELT1"]) # X方向に足し合わせ
xv_map = np.nansum(hdu.data, axis=1)*np.abs(hdu.header["CDELT2"]) # Y方向に足し合わせ
fig = plt.figure(figsize=(12, 6))
ax1 = fig.add_subplot(1, 2, 1)
ax1.imshow(yv_map, vmin=0, vmax=6)
ax1.set_xlabel("Y [pix]")
ax1.set_ylabel("V [pix]")
ax2 = fig.add_subplot(1, 2, 2)
ax2.imshow(xv_map, vmin=0, vmax=6)
ax2.set_xlabel("X [pix]")
ax2.set_ylabel("V [pix]")
plt.show()

軸を座標などにしたい場合は、前回 のようにやるか、または以下の指定範囲積分で hdu を作っているので、そちらで行いましょう。

指定範囲積分

大抵は指定した座標の範囲を積分したいかと思います。
まずは位置–速度の hdu を作る関数を定義します。Y 方向に積分するものと X 方向に積分する２つです。plot するのは 2D ですが、aplpy の仕様の都合上、あえて 3D にしています。

def make_xv_hdu(hdu, y_start_wcs, y_end_wcs, w):
    import copy
    data = hdu.data
    header = copy.deepcopy(hdu.header)
    x_tempo_wcs, y_tempo_wcs, v_tempo_wcs = w.wcs_pix2world(0, 0, 0, 0)
    x_start_ch, y_start_ch, v_start_ch = w.wcs_world2pix(x_tempo_wcs, y_start_wcs, v_tempo_wcs, 0)
    x_end_ch, y_end_ch, v_end_ch = w.wcs_world2pix(x_tempo_wcs, y_end_wcs, v_tempo_wcs, 0)
    y_start_ch, y_end_ch = int(round(float(y_start_ch))), int(round(float(y_end_ch)))
    
    data_xv = np.nansum(data[:, y_start_ch:y_end_ch+1, :], axis=1)*np.abs(header["CDELT2"])
    data_xv = data_xv.reshape(header["NAXIS3"], 1, header["NAXIS1"])

    header["CRVAL3"] = header["CRVAL3"]/1000.0
    header["CDELT3"] = header["CDELT3"]/1000.0
    header["CUNIT3"] = "km/s"
    header["CTYPE3"] = "VELOCITY"
    header["NAXIS2"] = 1
    
    new_hdu = fits.PrimaryHDU(data_xv, header)
    
    return new_hdu


def make_yv_hdu(hdu, x_start_wcs, x_end_wcs, w):
    import copy
    data = hdu.data
    header = copy.deepcopy(hdu.header)
    x_tempo_wcs, y_tempo_wcs, v_tempo_wcs = w.wcs_pix2world(0, 0, 0, 0)
    x_start_ch, y_start_ch, v_start_ch = w.wcs_world2pix(x_start_wcs, y_tempo_wcs, v_tempo_wcs, 0)
    x_end_ch, y_end_ch, v_end_ch = w.wcs_world2pix(x_end_wcs, y_tempo_wcs, v_tempo_wcs, 0)
    x_start_ch, x_end_ch = int(round(float(x_start_ch))), int(round(float(x_end_ch)))

    data_yv = np.nansum(data[:, :, x_start_ch:x_end_ch+1], axis=2)*np.abs(header["CDELT1"])
    data_yv = data_yv.reshape(header["NAXIS3"], header["NAXIS2"], 1)

    header["CRVAL3"] = header["CRVAL3"]/1000.0
    header["CDELT3"] = header["CDELT3"]/1000.0
    header["CUNIT3"] = "km/s"
    header["CTYPE3"] = "VELOCITY"
    header["CTYPE1"] = "GLON-TAN" # 注意
    header["CTYPE2"] = "GLAT-TAN" #
    header["NAXIS1"] = 1

    new_hdu = fits.PrimaryHDU(data_yv, header)
    
    return new_hdu

yv ですが、projection が SFL だとバグるようですので (2021年6月現在)、TAN に変更しています。RADec の場合は消してください。
以下のように使います。

xv_hdu = make_xv_hdu(hdu, -0.9, 0.9, w) # Y 座標 (degree) を下, 上の順で入れる
yv_hdu = make_yv_hdu(hdu, 31.9, 30.1, w) # X 座標 (degree) を左, 右の順で入れる

fig = plt.figure(figsize=(12, 6))
f = aplpy.FITSFigure(xv_hdu, slices=[0], figure=fig, convention='wells', dimensions=[0, 2])
f.show_colorscale(vmin=0.0, vmax=10.0, aspect='auto', stretch="linear", cmap="nipy_spectral")
f.add_colorbar()
f.axis_labels.set_xtext("Galactic Longitude")
f.axis_labels.set_ytext("$V_{\mathrm {LSR}}$ (km s$^{-1}$)")
plt.show()

fig = plt.figure(figsize=(6, 12))
f = aplpy.FITSFigure(yv_hdu, slices=[0], figure=fig, convention='wells', dimensions=[2, 1], aspect='auto')
f.show_colorscale(vmin=0.0, vmax=10.0, aspect='auto', stretch="linear", cmap="nipy_spectral")
f.add_colorbar()
f.axis_labels.set_xtext("$V_{\mathrm {LSR}}$ (km s$^{-1}$)")
f.axis_labels.set_ytext("Galactic Latitude")
f.tick_labels.set_yposition("left")
plt.show()

Python2系の時は位置–速度図でも recenter が使えたのですが、仕様が変わって機能しなくなりました。もしどうしても使いたい場合は 2系 で実行するか、matplotlib.pyplot で plot しましょう。show_markers 等に関しては、coords_frame="pixel" にして pixel 指定をすれば一応可能です。

また、projection の関係で下のように、変な第二軸が描かれることがあります。

このような際は

gca = plt.gca()
gca.tick_params(axis=1,# vyの場合は0
               labelbottom=False,
               labelleft=False,
               labelright=False,
               labeltop=False,
               bottom=False,
               left=False,
               right=False,
               top=False)

を f.show_colorscale(...) の直後に入れてください。

以上です。

リンク
目次