本記事では、pycupid の基本的な使い方について記述します。pycupid には GaussClump, ClumpFind, FellWalker, Reinhold の４つの構造同定アルゴリズムが入っています。

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

インストール

pycupid (ver. 0.1.4) は2021年6月現在、Python2.7, 3.4, 3.5, 3.6 でしか動作しません。

pip install pycupid

もしくは

pip2.7 install pycupid

などでインストールしましょう。

基本的な使い方

今回は簡単のため 2D fits での使い方について紹介します。3D (cube) でも基本的には同じです。
まずは例によって必要なものを import し、fitsを読み込みます。

from astropy.io import fits
import numpy as np
from matplotlib import pyplot as plt
import aplpy
from astropy.wcs import WCS
from pycupid import clumpfind, fellwalker, reinhold, gaussclumps

hdu = fits.open("~/your/fits/dir/FGN_03100+0000_2x2_13CO_v1.00_cube.fits")[0] # 3D
w = WCS(hdu.header)

2D fitsにするため、以下のような関数を定義します。

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

以下のように 2D fits を作成します。

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

このデータに clumpfind をかけます。一番シンプルな方法は、

cf = clumpfind(integ_hdu.data, 20).astype("float32")

これだけです。cf は同定された構造の ID が入った numpy array です。引数の 20 (今はてきとう) は r.m.s. を意味しています。デフォルトでは、この値の2倍が最小コントアレベルとコントアレベル間隔に設定されます。

同定されなかった場所 (ピクセル) にはマイナスの大きな値が入っているので、NaN にします。

cf[cf==cf.min()] = np.nan

ID map を plot してみます。

fig = plt.figure(figsize=(8, 8))
plt.imshow(cf, cmap="jet", origin="lower left")
plt.show()

何かがおかしいです。
どうやらバグ (仕様?) のようで、以下を実行すると元の fits と同じになります (3D でも同様の仕様を確認)。
(pycupid (ver. 0.1.4) 2021年6月現在)

cf_2 = cf.reshape(cf.shape[::-1]).T

もう一度 ID map を plot してみます。

fig = plt.figure(figsize=(8, 8))
plt.imshow(cf_2, cmap="jet", origin="lower left")
plt.show()

正しくなりました。

同定した構造の数は

print(np.nanmax(cf_2))
#1335.0

これでわかります。

特定の構造 (例えば ID=1000) を取ってきたい場合は、

ind = np.where(cf_2==1000)
print(ind)
#(array([452, 452, 452, 452, 453, 453, 453, 453, 454, 454, 454, 454, 455,
#       455, 455, 456, 456, 456]), array([337, 338, 339, 340, 337, 338, 339, 340, 337, 338, #339, 340, 337,
#       338, 339, 337, 338, 339]))

とすると、tuple でその構造の index が返ってきます。ここまで来れば、他の記事を参考にして自由に好きな解析ができると思います。

しかしながら、np.where は非常に遅いので、以下のようなコードを使うことをお勧めします (参考)。

from scipy.sparse import csr_matrix

def compute_M(data):
    cols = np.arange(data.size)
    return csr_matrix((cols, (data.ravel(), cols)),
                      shape=(data.max() + 1, data.size))

def get_indices_sparse(data):
    M = compute_M(data)
    return [np.unravel_index(row.data, data.shape) for row in M]

ind_list = get_indices_sparse(cf_2.astype("int"))

また、例えば、コントアを引きたい場合は、

fig = plt.figure(figsize=(8, 8))
f = aplpy.FITSFigure(integ_hdu, slices=[0], convention='wells', figure=fig)
f.show_colorscale(stretch='linear', cmap="Greys", aspect="equal")
f.show_contour(cf_hdu, levels=np.arange(1, np.nanmax(cf_2)+1), colors="c")
f.recenter(30.75, 0.0, width=0.5, height=0.5)
plt.show()

などでできます。

同定パラメータを変更する

辞書形式で詳細な設定を渡すことができます。詳しくはドキュメントを参照してください。
例えば、

config = {}
config["DELTAT"] = "3*RMS"
config["MINPIX"] = 25
config["LEVEL1"] = 30
config["LEVEL2"] = 50
config["LEVEL3"] = 100
config["LEVEL4"] = 150
config["LEVEL5"] = 200
config["LEVEL6"] = 250
config["LEVEL7"] = 300

cf_3 = clumpfind(integ_hdu.data, 20, config=config).astype("float32")

などといった感じで使います。

ClumpFind 以外のアルゴリズム (GaussClump, FellWalker, Reinhold) も概ね同様かと思います。

以上です。

リンク
目次