天文データ解析入門 その21 (fitsをムービーで保存する)

本記事では、fits 画像をムービー形式で保存する (速度チャンネルマップをムービーにする) 方法について紹介します。

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

まず cv2 をインストールします。

pip install opencv-python

これで大抵はインストールできる気がしますが、うまくいかない方は、自身の環境にあった方法を google 検索しましょう。

必要なものを import します。

from astropy.io import fits
from astropy.wcs import WCS
import aplpy
import numpy as np
import matplotlib.pyplot as plt
import os
import shutil
import glob
import cv2

積分強度図を for 文で作って、それらをムービーにします。
必要な関数を定義します。

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)*header["CDELT3"]/1000.0
    header = del_header_key(header, ["CRVAL3", "CRPIX3", "CRVAL3", "CDELT3", "CUNIT3", "CTYPE3", "CROTA3", "NAXIS3"])
    header["NAXIS"] = 2
    new_hdu = fits.PrimaryHDU(new_data, header)
    return new_hdu

def make_figure_integ():
    fig = plt.figure(figsize=(16, 16))
    f = aplpy.FITSFigure(hdu_integ, slices=[0],convention='wells', figure=fig)
    f.show_colorscale(vmin=color_min, vmax=color_max, stretch=stretch, cmap=cmap_integ, aspect="equal")
    f.add_colorbar()
    f.colorbar.show()
    f.colorbar.set_width(0.2)
    f.colorbar.set_axis_label_text("K km s$^{-1}$")
    f.show_contour(hdu_integ, levels=levels, colors="k", linewidths=0.5)
    f.set_title(title, size=20, family='serif')
    plt.tight_layout()
    f.save(os.path.join(save_name)+'.jpeg', dpi=200)
    plt.clf()

def create_movie():
    vw = cv2.VideoWriter(movie_name, fmt, frame_rate, (shape_x, shape_y))
    for jpeg in jpeg_list:
        im = cv2.imread(jpeg)
        vw.write(im)
    vw.release()

fits を読み込みます。

fits_name = "./your/fits/dir/FGN_03100+0000_2x2_13CO_v1.00_cube.fits"
hdu = fits.open(fits_name)[0]
h = hdu.header
w = WCS(h)

積分強度図を作ります。

save_dir = "integs"
if os.path.exists(save_dir):
    shutil.rmtree(save_dir)
    os.mkdir(save_dir)
else:
    os.mkdir(save_dir)
start_v = 0.0 # km/s (最初の channel の速度)
interval_ch = 2 # channel (積分する channel 数)

color_min, color_max = 0, 50
cmap_integ = "gist_earth_r"
stretch = "linear"

levels = np.array([1,2,3,4,5,6,7,8])*3.0

for i in range(100): # 多い場合は multi process 推奨
    start_v_i = start_v + i*(int(interval_ch)*h["CDELT3"]/1000.0)
    end_v_i = start_v_i + int(interval_ch)*h["CDELT3"]/1000.0
    hdu_integ=  make_new_hdu_integ(hdu, start_v_i, end_v_i, w)
    save_name = os.path.join(save_dir, "%s_%.2f_%.2f"%(str(i).zfill(5), start_v_i, end_v_i))
    title = "\n\n%.2f $-$ %.2f km/s"%(start_v_i - h["CDELT3"]/1000.0/2.0, end_v_i - h["CDELT3"]/1000.0/2.0)
    make_figure_integ()

ムービーの設定をします。

movie_name = "integs.mp4" # 保存するムービーの名前
fmt = cv2.VideoWriter_fourcc(*"mp4v") # ムービーのフォーマット 
jpeg_list = sorted(glob.glob(os.path.join(save_dir, "*.jpeg"))) # ムービーにする jpeg のリスト
shape_x, shape_y = cv2.imread(jpeg_list[0]).shape[1], cv2.imread(jpeg_list[0]).shape[0] # ムービーの shape
frame_rate = 1 # フレームレート

実行します。

create_movie()

integs.mp4 というムービーが生成されているはずです。

以上です。

リンク
目次