天文データ解析入門 その20 (pycupid: pythonでclumpfind等を使う)

本記事では、[pycupid](https://pycupid.readthedocs.io/en/latest/index.html) の基本的な使い方について記述します。pycupid には GaussClump, ClumpFind, FellWalker, Reinhold の４つの構造同定アルゴリズムが入っています。

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

```python
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にするため、以下のような関数を定義します。

```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
```

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

```python
integ_hdu = make_new_hdu_integ(hdu, 25.0, 125.0, w)
```
このデータに clumpfind をかけます。一番シンプルな方法は、

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

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

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

ID map を plot してみます。

```python
fig = plt.figure(figsize=(8, 8))
plt.imshow(cf, cmap="jet", origin="lower left")
plt.show()
```
![ダウンロード (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/cca590ea-8712-1520-50e0-73e95516e6f5.png)
何かがおかしいです。
どうやらバグ (仕様?) のようで、以下を実行すると元の fits と同じになります (3D でも同様の仕様を確認)。
(pycupid (ver. 0.1.4) 2021年6月現在) 

```python
cf_2 = cf.reshape(cf.shape[::-1]).T
```
もう一度 ID map を plot してみます。

```python
fig = plt.figure(figsize=(8, 8))
plt.imshow(cf_2, cmap="jet", origin="lower left")
plt.show()
```
![ダウンロード (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/eeb0e5b2-48b3-e3f3-7b35-8f13017116ac.png)
正しくなりました。

同定した構造の数は

```python
print(np.nanmax(cf_2))
#1335.0
```
これでわかります。

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

```python
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 が返ってきます。ここまで来れば、[他の記事](https://qiita.com/Shinji_Fujita/items/ecb869e5b6a9fa468483)を参考にして自由に好きな解析ができると思います。

しかしながら、np.where は非常に遅いので、以下のようなコードを使うことをお勧めします ([参考](https://base64.work/so/python/418382))。

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

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

```python
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()
```
などでできます。
![ダウンロード (4).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/a5b516b8-f67b-a498-611b-8ae0570c39e3.png)

# 同定パラメータを変更する
辞書形式で詳細な設定を渡すことができます。詳しくは[ドキュメント](https://pycupid.readthedocs.io/en/latest/parameters.html)を参照してください。
例えば、

```python
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) も概ね同様かと思います。

以上です。

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