天文データ解析入門 その24 (FilFinder (fil_finder) の使い方)

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

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

本記事では、[FilFinder](https://fil-finder.readthedocs.io/en/latest/index.html) の使い方について記述します。

# 準備
今回は、比較の例としてHerschelのアーカイブデータも用います。データは
http://herschel.esac.esa.int/Science_Archive.shtml
から、フィラメントが顕著な領域である Taurus の、
4:19:52.0 +27:12:30.0 (FK5)
あたりにします。
観測 ID　は　1342202254　です。level2_5の、
hspireplw274_25pxmp_0421_p2748_1342190616_1342202090_1462391077877.fits
を使用します。

まず、FilFinder の install ですが、python のカーネルが 2.7 や 3.7 等以外だと失敗します。少なくとも 3.6 では失敗しました。カーネルの追加方法ですが、各々の環境によります。例えば pyenv を使用している方は、3.7.4 などを install し、ターミナルで作業するディレクトリに移動して
```
pyenv local 3.7.4
```
とした後に、
```
jupyter notebook
```
で jupyter を起動してください。

# 距離とビームサイズが既知で、デフォルトのパラメータで実行したい場合 (とりあえず回したい人)

まずは必要なものを import します。入っていないものがあったらこの機会に pip しておきましょう (※作業ディレクトリで pip しましょう)。
```python
from astropy.io import fits
import numpy as np
from matplotlib import pyplot as plt
import aplpy
from astropy.wcs import WCS
import astropy.units as u
from astropy.coordinates import SkyCoord
import fil_finder
from fil_finder import FilFinder2D, Filament2D
from astropy.io import fits
import tqdm
import copy
import random
from matplotlib.colors import rgb2hex, Normalize
from matplotlib.colorbar import ColorbarBase
import pandas as pd
# import warnings # warningが気になる方はこちらも
# warnings.simplefilter('ignore')
```

念のため python のバージョンを確認します。
```python
import sys
sys.version
# '3.7.4 (default,...
```

```python
hdu = fits.open("./hspireplw274_25pxmp_0421_p2748_1342190616_1342202090_1462391077877.fits")[1] # 普通は [0] です。
w = WCS(hdu)
```
見た目を確認します。
```python
fig = plt.figure(figsize=(8, 8))
f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(cmap="Greys", aspect="equal")
```
![ダウンロード (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/c444bffa-97e5-2805-b740-3b70a3d2cc69.png)

header も確認しておきます。
```python
hdu.header
# XTENSION= 'IMAGE   '           / Java FITS: Wed May 04 11:32:01 CEST 2016       
# BITPIX  =                  -64                                                  
# NAXIS   =                    2 / Dimensionality                                 
# NAXIS1  =                 1461                                                  
# NAXIS2  =                 1280  
# ...
# ...
```
距離やビームサイズが入っていないようですので、こちらで用意します。
```python
distance = 150 * u.pc
BMAJ = 30.0 * u.arcsec
```
では早速 FilFinder2D に突っ込みます。
```python
fil = FilFinder2D(hdu, distance=distance, beamwidth=BMAJ)
# ---------------------------------------------------------------------------
# TypeError 
# ...
# ...
# TypeError: Input data is not of an accepted form:['numpy.ndarray', 'astropy.io.fits.PrimaryHDU', 'spectral_cube.Projection', 'spectral_cube.Slice', 'astropy.units.Quantity']
```
Herschel の fits　が少しクセがあるので、以下のようにして回避します。
```python
fil = FilFinder2D(fits.PrimaryHDU(hdu.data, hdu.header), distance=distance, beamwidth=BMAJ)
```
通りました。
次に、実際にフィラメント同定をかける image を生成します。これは、おそらく低い強度のフィラメントも同定するためです。詳しくは[ドキュメント](https://fil-finder.readthedocs.io/en/latest/tutorial.html)を参照してください。skip_flatten=True をオプションに入れることでスキップもできます。
```python
fil.preprocess_image()
```
次に、 mask を生成します。今は全てデフォルトのパラメータでやっていますが、パラメータを変更すると結果はかなり変動します。詳しくは[ドキュメント](https://fil-finder.readthedocs.io/en/latest/tutorial.html)を参照してください。
```python
plt.figure(figsize=(8, 8))
fil.create_mask(verbose=True)
```
![ダウンロード (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/5f12c9b8-58d7-d401-ade9-3219cd946910.png)
スケルトン化します。
```python
plt.figure(figsize=(8, 8))
fil.medskel(verbose=True)
```
![ダウンロード (3).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/e2e401a1-f55c-6129-b119-05053121321b.png)
スケルトンを解析します。少し時間がかかります。また、パラメータも複数あり、いくつかは切り落とされます。詳しくは[ドキュメント](https://fil-finder.readthedocs.io/en/latest/tutorial.html)を参照してください。
```python
fil.analyze_skeletons()
fil.exec_rht()
```
残ったフィラメントの数を確認します。
```python
len(fil.filaments)
# 80
```
map で確認します。
```python
fil_skeleton= fil.skeleton
fil_skeleton = fil_skeleton.astype("int")

fig = plt.figure(figsize=(16, 16))
f = aplpy.FITSFigure(hdu, slices=[0], figure=fig)
f.show_colorscale(vmin = 1, vmax = 100, stretch="log", cmap="Greys")
f.add_colorbar()
f.show_contour(fits.PrimaryHDU(fil_skeleton, hdu.header),colors="r",linewidths=2,slices=[0], alpha=0.4)
```
![ダウンロード (4).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/b078830b-c150-6251-331b-68f57f7178b1.png)

今度は最長経路を plot してみます。
```python
fil_skeleton_longpath = fil.skeleton_longpath
fil_skeleton_longpath = fil_skeleton_longpath.astype("int")

fig = plt.figure(figsize=(16, 16))
f = aplpy.FITSFigure(hdu, slices=[0], figure=fig)
f.show_colorscale(vmin = 1, vmax = 100, stretch="log", cmap="Greys")
f.add_colorbar()
f.show_contour(fits.PrimaryHDU(fil_skeleton_longpath, hdu.header),colors="r",linewidths=2,slices=[0], alpha=0.4)
```
![ダウンロード (5).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/d3b38388-0c62-8ff1-d05a-a67228fdaa5a.png)

ここからの解析は、ページの後半で紹介します。

# 距離やビームサイズが不明の場合

上の例 (デフォルトパラメータ) でなぜフィラメントがそこそこ上手く同定されていたかと言うと、フィラメントの幅がおおよそ 0.1-0.2 pc 程度であると内部にコーディングされているからです。
距離やビームサイズが不明の場合について紹介します。
つまり、与える image は numpy の 2D array だけです。

まずは必要なものを import します。入っていないものがあったらこの機会に pip しておきましょう (※作業ディレクトリで pip しましょう)。
```python
from astropy.io import fits
import numpy as np
from matplotlib import pyplot as plt
import aplpy
from astropy.wcs import WCS
import astropy.units as u
from astropy.coordinates import SkyCoord
import fil_finder
from fil_finder import FilFinder2D, Filament2D
from astropy.io import fits
import tqdm
import copy
import random
from matplotlib.colors import rgb2hex, Normalize
from matplotlib.colorbar import ColorbarBase
import pandas as pd
# import warnings # warningが気になる方はこちらも
# warnings.simplefilter('ignore')
```

```python
fil = FilFinder2D(image=hdu.data)
```
同様に image を作成します。この例では平坦化をスキップしています。
```python
fil.preprocess_image(skip_flatten=True)
```

mask　を作成します。パラメータについては[ドキュメント](https://fil-finder.readthedocs.io/en/latest/tutorial.html)を参照してください。
```python
plt.figure(figsize=(8, 8))
fil.create_mask(verbose=True, adapt_thresh=5*u.pix,  smooth_size=1*u.pix, size_thresh=100*u.pix*u.pix,  glob_thresh=10.0, fill_hole_size=10*u.pix*u.pix)
```
![ダウンロード (6).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/6c54d3f5-82a1-a134-7622-b9508e81e1f7.png)

```python
plt.figure(figsize=(8, 8))
fil.medskel(verbose=True)
```
![ダウンロード (7).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/5f98ae81-46d0-5bcc-0f29-a01f2b26a4fa.png)

スケルトン解析をします。パラメータについては[ドキュメント](https://fil-finder.readthedocs.io/en/latest/tutorial.html)を参照してください
```python
fil.analyze_skeletons(skel_thresh=10*u.pix)
fil.exec_rht()
```
先ほどと同様、フィラメントの数を確認します。
```python
len(fil.filaments)
# 205
```
先ほどと同様、見た目も確認します。
```python
fil_skeleton= fil.skeleton
fil_skeleton = fil_skeleton.astype("int")

fig = plt.figure(figsize=(16, 16))
f = aplpy.FITSFigure(hdu, slices=[0], figure=fig)
f.show_colorscale(vmin = 1, vmax = 100, stretch="log", cmap="Greys")
f.add_colorbar()
f.show_contour(fits.PrimaryHDU(fil_skeleton, hdu.header),colors="r",linewidths=2,slices=[0], alpha=0.4)
```
一旦これでフィラメントを同定するまでは終わりました (パラメータチューニング頑張ってください)。

# 解析例

同定したフィラメントに対して、その物理量を取り出す最も簡単な方法は以下のような方法です。

```python
fil.lengths()
# [53.769553, 68.941125, 264.42136, …, 46.526912, 18.142136, 30.970563]pix
```
```python
fil.curvature
# [0.95357014, 0.65610733, 1.1758287, …, 0.71241929, 0.57067502, 0.86255477]rad
```
```python
fil.orientation
# [0.49783893, −0.22334848, −0.0044566261, …, 0.85370036, −0.58684442, 1.4698403]rad
```
branch に関しては以下のようにして情報を取り出せます。
```python
fil.branch_properties["length"]
# [<Quantity [0.5       , 4.24264069, 1.41421356, 2.82842712, 9.41421356,
# ...
# ...
```
```python
fil.branch_properties["intensity"] # おそらく Total ではない
# [array([10.44615376, 11.01044167, 10.68661623, 10.75544016, 10.37313803,
# ...
# ...
```
```python
fil.branch_properties["pixels"]
# [[array([[ 1, 19]]),
# ...
# ...
```
```python
fil.branch_properties["number"]
# array([ 20,  14,  43,   6,  34,   3,  70,  10,  21,   1,   3,  21,   9,
# ...
# ...
```
フィラメントの幅に関しては、以下を実行します。
```python
fil.find_widths()
```

```python
fil.widths()
# (<Quantity [ 9.32484571,  9.20587192, 11.01286974,  7.03626262,  8.4960545 ,
# ...
# ...
```
0番目が FWHM (pix), 1番目がそのエラー (pix) です。

# 解析例 (branch ごと)

branch ごとに解析したい場合もあるかと思います。その場合、以下のようにします。
(※このあたりからドキュメントの範疇を超えているので、間違っていたらごめんなさい。)

まず、パラメータをセットします。
```python
min_length = 10 * u.pix # この長さを超えていなければフィラメントとして認めない
max_dist = 15 * u.pix # フィラメントの幅を見積もるために強度をサーチしに行く最大距離
```
以下を実行します。少し時間がかかります。
```python
ori_branches_array = []
total_intensity_branches_array = []
length_branches_array = []
width_branches_array = []
width_err_branches_array = []
pos_x_branches_array, pos_y_branches_array = [], []
for fil_tempo in tqdm.tqdm(fil.filaments):
    y_start, x_start = fil_tempo.pixel_extents[0][0], fil_tempo.pixel_extents[0][1]
    fil_tempo.rht_branch_analysis()
    ori_branches_tempo = np.rad2deg(fil_tempo.orientation_branches.value)
    for i in range(len(fil_tempo.branch_pts())):
        pixels_tempo = fil_tempo.branch_properties["pixels"][i]
        y_pix_all = pixels_tempo[:,0]
        x_pix_all = pixels_tempo[:,1]
        pos_x_tempo = int(np.nanmax(x_pix_all)/2-np.nanmin(x_pix_all)/2)+np.nanmin(x_pix_all)+x_start
        pos_y_tempo = int(np.nanmax(y_pix_all)/2-np.nanmin(y_pix_all)/2)+np.nanmin(y_pix_all)+y_start
        ori_branch_tempo = ori_branches_tempo[i]
        pos_x_wcs_tempo, pos_y_wcs_tempo = w.wcs_pix2world(pos_x_tempo, pos_y_tempo, 0)
        pos_x_branches_array.append(pos_x_tempo)
        pos_y_branches_array.append(pos_y_tempo)
        if fil_tempo.branch_properties["length"][i]>min_length:
            ori_branches_array.append(ori_branch_tempo)
            branch_tempo_coords = (y_pix_all+y_start, x_pix_all+x_start)
            branch_tempo = Filament2D(branch_tempo_coords)
            branch_tempo.skeleton_analysis(hdu.data)
            branch_tempo.width_analysis(hdu.data, deconvolve_width=False, max_dist=max_dist)
            width_tempo, width_err_tempo = branch_tempo.radprof_fwhm()
            total_intensity_branches_array.append(branch_tempo.total_intensity())
            length_branches_array.append(fil_tempo.branch_properties["length"][i].value)
            width_branches_array.append(width_tempo.value)
            width_err_branches_array.append(width_err_tempo.value)
        else:
            ori_branches_array.append(np.nan)
            total_intensity_branches_array.append(np.nan)
            length_branches_array.append(np.nan)
            width_branches_array.append(np.nan)
            width_err_branches_array.append(np.nan)
ori_branches_array = np.array(ori_branches_array)
total_intensity_branches_array = np.array(total_intensity_branches_array)
length_branches_array = np.array(length_branches_array)
width_branches_array = np.array(width_branches_array)
width_err_branches_array = np.array(width_err_branches_array)
pos_x_branches_array = np.array(pos_x_branches_array)
pos_y_branches_array = np.array(pos_y_branches_array)
```
何をしているかと言うと、branch ごとに、branch をフィラメントオブジェクトと偽装して特殊召喚しています。
詳しくは中身を読み解いてください。

orientation を確認します。
```python
fig = plt.figure(figsize=(8, 8))
plt.hist(ori_branches_array)
plt.xlabel("Orientation to the North (CCW) (degree)")
plt.ylabel("Number of filaments")
```
![ダウンロード (8).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/473a35a9-c416-4c14-6c09-643cb8610c5b.png)

確か基準が北 (numpy 2D array での、上) で、反時計回りだったと思うのですが、曖昧なので目で見て確認します。

```python
vmin1, vmax1 = 1, 100
vmin2, vmax2 = -90, 90
plt.rcParams["font.family"] = "Serif"

fig = plt.figure(figsize=(16, 16))
f = aplpy.FITSFigure(hdu, slices=[0], figure=fig, subplot=[0.1, 0.1, 0.8, 0.8])
f.show_colorscale(vmin = vmin1, vmax = vmax1, stretch="log", cmap="Greys")
f.add_colorbar()
f.colorbar.show()
f.colorbar.set_width(0.2)
f.colorbar.set_font(size=15)
f.colorbar.set_axis_label_text("MJy/str")
f.colorbar.set_axis_label_font(size=15)

pos_x_wcs_list = []
pos_y_wcs_list = []
color_tempo_list = []
cmap = plt.get_cmap("jet")
for fil_tempo in tqdm.tqdm(fil.filaments):
    y_start, x_start = fil_tempo.pixel_extents[0][0], fil_tempo.pixel_extents[0][1]
    fil_tempo.rht_branch_analysis()
    ori_branches_tempo = np.rad2deg(fil_tempo.orientation_branches.value)
    for i in range(len(fil_tempo.branch_pts())):
        if fil_tempo.branch_properties["length"][i]>min_length:
            pixels_tempo = fil_tempo.branch_properties["pixels"][i]
            y_pix_all = pixels_tempo[:,0]
            x_pix_all = pixels_tempo[:,1]
            pos_x_tempo = int(np.nanmax(x_pix_all)/2-np.nanmin(x_pix_all)/2)+np.nanmin(x_pix_all)+x_start
            pos_y_tempo = int(np.nanmax(y_pix_all)/2-np.nanmin(y_pix_all)/2)+np.nanmin(y_pix_all)+y_start
            ori_branch_tempo = ori_branches_tempo[i]
            pos_x_wcs_tempo, pos_y_wcs_tempo = w.wcs_pix2world(pos_x_tempo, pos_y_tempo, 0)
            branch_tempo_coords = (y_pix_all+y_start, x_pix_all+x_start)
            branch_tempo = Filament2D(branch_tempo_coords)
            branch_tempo.skeleton_analysis(hdu.data)
            branch_tempo.width_analysis(hdu.data, deconvolve_width=False)
            pos_x_wcs_array, pos_y_wcs_array = w.wcs_pix2world(x_pix_all+x_start, y_pix_all+y_start, 0)
            pos_x_wcs_list += list(pos_x_wcs_array)
            pos_y_wcs_list += list(pos_y_wcs_array)
            color_tempo = rgb2hex(cmap(Normalize(vmin=vmin2, vmax=vmax2)(ori_branch_tempo)))
            color_tempo_list += [color_tempo]*len(pos_x_wcs_array)
f.show_markers(pos_x_wcs_list, pos_y_wcs_list, marker="s", s=5, facecolor=color_tempo_list, edgecolor="none", zorder=11)

ax = fig.add_subplot(1, 1, 1)
norm = Normalize(vmin=vmin2, vmax=vmax2)
cb = ColorbarBase(ax=ax, cmap="jet", orientation="horizontal", norm=norm)
cb.ax.set_position([0.1, 0.92, 0.8, 0.02])
cb.set_label("Orientation to the North (CCW) (degree)", labelpad=-70, size=15)

fig.patch.set_facecolor('w')
#f.save('./neko.png', dpi=200)
```

![ダウンロード (9).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/f596e53b-afc6-f85b-ab3b-03c4af5fd7c3.png)

正解のようです。
今は　color_tempo　のところに ori_branch_tempo を与えていますが、例えば ori_branch_tempo を branch_tempo.total_intensity() などに置き換えて、vmin2 と vmax2 を変更すれば、いろいろな map が描けます。

ちなみに最後の
```python
ax = fig.add_subplot(1, 1, 1)
norm = Normalize(vmin=vmin2, vmax=vmax2)
cb = ColorbarBase(ax=ax, cmap="jet", orientation="horizontal", norm=norm)
cb.ax.set_position([0.1, 0.92, 0.8, 0.02])
```
は、カラーバーのみを顕現させる小技なので、覚えておくと良いと思います。

以上です。

追記：
matplotlib のバージョンによっては、
```python
from matplotlib import cm
cb = ColorbarBase(ax=ax, cmap=cm.jet, orientation="horizontal", norm=norm)
```
にする必要があるかもしれません。

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