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

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

pip install matplotlib==3.5.2

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

本記事では、FilFinder の使い方について記述します。

準備

今回は、比較の例として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 しましょう)。

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 のバージョンを確認します。

import sys
sys.version
# '3.7.4 (default,...

hdu = fits.open("./hspireplw274_25pxmp_0421_p2748_1342190616_1342202090_1462391077877.fits")[1] # 普通は [0] です。
w = WCS(hdu)

見た目を確認します。

fig = plt.figure(figsize=(8, 8))
f = aplpy.FITSFigure(hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(cmap="Greys", aspect="equal")

header も確認しておきます。

hdu.header
# XTENSION= 'IMAGE   '           / Java FITS: Wed May 04 11:32:01 CEST 2016       
# BITPIX  =                  -64                                                  
# NAXIS   =                    2 / Dimensionality                                 
# NAXIS1  =                 1461                                                  
# NAXIS2  =                 1280  
# ...
# ...

距離やビームサイズが入っていないようですので、こちらで用意します。

distance = 150 * u.pc
BMAJ = 30.0 * u.arcsec

では早速 FilFinder2D に突っ込みます。

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　が少しクセがあるので、以下のようにして回避します。

fil = FilFinder2D(fits.PrimaryHDU(hdu.data, hdu.header), distance=distance, beamwidth=BMAJ)

通りました。
次に、実際にフィラメント同定をかける image を生成します。これは、おそらく低い強度のフィラメントも同定するためです。詳しくはドキュメントを参照してください。skip_flatten=True をオプションに入れることでスキップもできます。

fil.preprocess_image()

次に、 mask を生成します。今は全てデフォルトのパラメータでやっていますが、パラメータを変更すると結果はかなり変動します。詳しくはドキュメントを参照してください。

plt.figure(figsize=(8, 8))
fil.create_mask(verbose=True)

スケルトン化します。

plt.figure(figsize=(8, 8))
fil.medskel(verbose=True)

スケルトンを解析します。少し時間がかかります。また、パラメータも複数あり、いくつかは切り落とされます。詳しくはドキュメントを参照してください。

fil.analyze_skeletons()
fil.exec_rht()

残ったフィラメントの数を確認します。

len(fil.filaments)
# 80

map で確認します。

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)

今度は最長経路を plot してみます。

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)

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

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

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

まずは必要なものを import します。入っていないものがあったらこの機会に pip しておきましょう (※作業ディレクトリで pip しましょう)。

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')

fil = FilFinder2D(image=hdu.data)

同様に image を作成します。この例では平坦化をスキップしています。

fil.preprocess_image(skip_flatten=True)

mask　を作成します。パラメータについてはドキュメントを参照してください。

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)

plt.figure(figsize=(8, 8))
fil.medskel(verbose=True)

スケルトン解析をします。パラメータについてはドキュメントを参照してください

fil.analyze_skeletons(skel_thresh=10*u.pix)
fil.exec_rht()

先ほどと同様、フィラメントの数を確認します。

len(fil.filaments)
# 205

先ほどと同様、見た目も確認します。

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)

一旦これでフィラメントを同定するまでは終わりました (パラメータチューニング頑張ってください)。

解析例

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

fil.lengths()
# [53.769553, 68.941125, 264.42136, …, 46.526912, 18.142136, 30.970563]pix

fil.curvature
# [0.95357014, 0.65610733, 1.1758287, …, 0.71241929, 0.57067502, 0.86255477]rad

fil.orientation
# [0.49783893, −0.22334848, −0.0044566261, …, 0.85370036, −0.58684442, 1.4698403]rad

branch に関しては以下のようにして情報を取り出せます。

fil.branch_properties["length"]
# [<Quantity [0.5       , 4.24264069, 1.41421356, 2.82842712, 9.41421356,
# ...
# ...

fil.branch_properties["intensity"] # おそらく Total ではない
# [array([10.44615376, 11.01044167, 10.68661623, 10.75544016, 10.37313803,
# ...
# ...

fil.branch_properties["pixels"]
# [[array([[ 1, 19]]),
# ...
# ...

fil.branch_properties["number"]
# array([ 20,  14,  43,   6,  34,   3,  70,  10,  21,   1,   3,  21,   9,
# ...
# ...

フィラメントの幅に関しては、以下を実行します。

fil.find_widths()

fil.widths()
# (<Quantity [ 9.32484571,  9.20587192, 11.01286974,  7.03626262,  8.4960545 ,
# ...
# ...

0番目が FWHM (pix), 1番目がそのエラー (pix) です。

解析例 (branch ごと)

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

まず、パラメータをセットします。

min_length = 10 * u.pix # この長さを超えていなければフィラメントとして認めない
max_dist = 15 * u.pix # フィラメントの幅を見積もるために強度をサーチしに行く最大距離

以下を実行します。少し時間がかかります。

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 を確認します。

fig = plt.figure(figsize=(8, 8))
plt.hist(ori_branches_array)
plt.xlabel("Orientation to the North (CCW) (degree)")
plt.ylabel("Number of filaments")

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

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)

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

ちなみに最後の

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 のバージョンによっては、

from matplotlib import cm
cb = ColorbarBase(ax=ax, cmap=cm.jet, orientation="horizontal", norm=norm)

にする必要があるかもしれません。

リンク
目次