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

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

例として、国立天文台の [FUGIN](https://nro-fugin.github.io/) プロジェクトで得られた野辺山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 します。

```python
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 文で作って、それらをムービーにします。
必要な関数を定義します。

```python
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 を読み込みます。

```python
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)
```
積分強度図を作ります。

```python
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()
```
ムービーの設定をします。

```python
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 # フレームレート
```
実行します。

```python
create_movie()
```

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

以上です。

リンク
[目次](https://qiita.com/Shinji_Fujita/items/ecb869e5b6a9fa468483)