天文データ解析入門 その4 (積分強度図の作り方)

[以下追記] matplotlib のバージョンが3.7.0等だと、aplpy でエラーが起きるようです。
```python
TypeError: WCSAxes.__init__() got multiple values for argument 'wcs'
```

```
pip install matplotlib==3.5.2
```
などとしてダウングレードしてください。
[以上追記]

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

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

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

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

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

data の shape を確認します。

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

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

```python
spe = hdu.data[: ,424, 424]
plt.plot(range(len(spe)), spe)
plt.show()
```
![ダウンロード (4).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/ff8733a4-f24d-4083-b2b1-e96d82fbadb5.png)
ここで「:」は「全て」を意味します。また、python では 424 と書くと実際には 424+1=425 ピクセル目をとってくることに注意です。

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

```python
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 関数を使用します。

```python
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で割り算します。

```python
spe = hdu.data[: ,424, 424]
plt.plot(v_array/1000.0, spe)
plt.xlabel("Velocity [km/s]")
plt.ylabel("Tmb [K]")
plt.show()
```
![ダウンロード (5).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/5c39a137-ccd0-6c0f-1333-92e1ab2527ee.png)
これでスペクトルが描けました。

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

```python
integ_map = np.nansum(hdu.data, axis=0)*hdu.header["CDELT3"]/1000.0
plt.imshow(integ_map, vmin=0, vmax=500)
plt.show()
```
![ダウンロード (6).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/bbf8b5d7-2539-6dda-5603-ceb452d4c03f.png)

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

```python
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 をスライスします。

```python
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()
```
![ダウンロード (7).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/fad683a8-3388-de60-dac2-0b9b4639c59c.png)
余計なチャンネルを除いたので、その分ノイズが減ってちょっと綺麗になりました。

# 積分強度図を fits として保存する
上でやった処理を効率化するために以下のような関数を定義します。

```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)*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 を以下で作ります。

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

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

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

![ダウンロード (8).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/509a5ee7-98af-9ea3-16c4-146500b69965.png)

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

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

以上です。

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