天文データ解析入門 その14 (astrodendroの使い方: scimes編)

【以下追記】
numpy のバージョンが 1.23 以上だと、astrodendro 内のパラメータの計算の中で
```python
module 'numpy' has no attribute 'asscalar'
```
とエラーが出ます。
astrodendro はおそらくもうアップデートされないので、numpy のバージョンを下げる (もしくは自分でastrodendroの中身を修正する) 必要があります。例えば、一時的な解決策として、解析ディレクトリで
```
pip install 'numpy==1.22.4' -t ./
```
などとして対応してください。
【以上追記】
本記事では、以前の記事 ([天文データ解析入門 その6 (astrodendroの使い方: 基本編)](https://qiita.com/Shinji_Fujita/items/5b9b86aee8ed8cfb6960)) に関連して、[scimes](https://scimes.readthedocs.io/en/latest/) の使い方について紹介します。

# scimes の概要
簡単に言うと、astrodendro の Dendrogram.compute で生成した Dendrogram オブジェクトを入力すると、それに含まれる構造を k-means法 でクラスタリングしてくれるモジュールです。詳細は[ドキュメント](https://scimes.readthedocs.io/en/latest/algorithm.html) を参照してください。


# インストール方法
ドキュメントには

```
pip install scimes
```
とありますが、うまくいかないことが多いです。

[githubのページ](https://github.com/Astroua/SCIMES) に行って、「↓Code」から「DOWNLOAD ZIP」を選んでください。解凍した後、

```
cd SCIMES-master
python setup.py install
```
を実行します。
astropy_helpers がないと言われたら、

```
pip install astropy_helpers
```
を実行してください。

# 使用方法
必要なものを import し、fitsを読み込みます。fits は、上の SCIMES-master に入っている orion_12CO.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 astropy import units as u
from astrodendro import Dendrogram
from astrodendro import ppv_catalog
from scimes import SpectralCloudstering

hdu = fits.open("SCIMES-master/orion_12CO.fits")[0] # 3D
w = WCS(hdu.header)
```
普通に astrodendrogram を実行します。

```python
dendro_min_value = 0.3
dendro_min_delta = 0.6
dendro_min_npix = 4
d = Dendrogram.compute(hdu.data, min_value=dendro_min_value, min_delta=dendro_min_delta, min_npix=dendro_min_npix, verbose=False)
```
metadata を定義します。そしてカタログを作ります。

```python
metadata = {}
metadata['data_unit'] = u.Jy # u.K はなぜか対応していないので、とりあえず Jy にしておく。
cat = ppv_catalog(d, metadata)
```
以下でクラスタリングを実行します。

```python
dclust = SpectralCloudstering(d, cat, hdu.header)
```
デフォルトでは、"luminosity"(flux) と "volume" が使われるようです。どちらか片方にしたい場合などは criteria で指定します。

```python
dclust = SpectralCloudstering(d, cat, hdu.header, criteria=["luminosity"])
dclust = SpectralCloudstering(d, cat, hdu.header, criteria=["volume"])
```
また、デフォルトでは独立した leaf (つまり、trunk 兼 leaf) はクラスタとして認識しません。認識させたい場合は、SpectralCloudstering 実行時に save_isol_leaves=True を加えてください。
さらに、k-means のパラメータも調節可能です。例えば user_k で個数を、user_iter でイテレーション回数を指定できます (大抵思うようにいきませんが)。


以下を実行することで結果の樹形図が表示されます。

```python
dclust.showdendro()
```
![ダウンロード (4).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/12a844c7-20b7-1737-aa2e-62e824481541.png)
同定したクラスタの番号と、ランダムに割り当てた色を定義します。

```python
clusts = dclust.clusters
colors = dclust.colors
print(clusts)
print(colors)
#[11, 37, 43, ..., 428, 475, 755]
#['#1BA58D', '#5F7FE6', '#BDB864', ..., '#EAED72', '#77E1E4', '#B37499']
```
aplpy で plot します。

```python
hdu_integ = fits.open("./SCIMES-master/orion_12CO_mom0.fits")[0]

f = aplpy.FITSFigure(hdu_integ, figsize=(8, 6), convention='wells', slices=[0])
f.show_colorscale(cmap='gray_r', vmin=1, vmax=80, stretch = 'log')

count = 0
for c in clusts:
        mask = d[c].get_mask()
        mask_hdu = fits.PrimaryHDU(mask.astype('short'), hdu.header)
        mask_coll = np.amax(mask_hdu.data, axis = 0)
        mask_coll_hdu = fits.PrimaryHDU(mask_coll.astype('short'), hdu.header)        
        f.show_contour(mask_coll_hdu, colors=colors[count], linewidths=2, convention='wells', levels = [0], slices=[0])
        count += 1
f.add_colorbar()
```

![ダウンロード (5).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/00744188-0330-54ef-9196-abee9c609313.png)

```python
clusts_structure = [d[c] for c in clusts]
```
と定義してしまえば、[天文データ解析入門 その6 (astrodendroの使い方: 基本編)](https://qiita.com/Shinji_Fujita/items/5b9b86aee8ed8cfb6960) や [天文データ解析入門 その13 (astrodendroの使い方: 中級編)](https://qiita.com/Shinji_Fujita/items/c84ed54bbe0ced5121c8) の方法で色々解析ができます。

以上です。

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