LoginSignup
2
1

More than 1 year has passed since last update.

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

Last updated at Posted at 2021-06-15

【以下追記】
numpy のバージョンが 1.23 以上だと、astrodendro 内のパラメータの計算の中で

module 'numpy' has no attribute 'asscalar'

とエラーが出ます。
astrodendro はおそらくもうアップデートされないので、numpy のバージョンを下げる (もしくは自分でastrodendroの中身を修正する) 必要があります。例えば、一時的な解決策として、解析ディレクトリで

pip install 'numpy==1.22.4' -t ./

などとして対応してください。
【以上追記】
本記事では、以前の記事 (天文データ解析入門 その6 (astrodendroの使い方: 基本編)) に関連して、scimes の使い方について紹介します。

scimes の概要

簡単に言うと、astrodendro の Dendrogram.compute で生成した Dendrogram オブジェクトを入力すると、それに含まれる構造を k-means法 でクラスタリングしてくれるモジュールです。詳細はドキュメント を参照してください。

インストール方法

ドキュメントには

pip install scimes

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

githubのページ に行って、「↓Code」から「DOWNLOAD ZIP」を選んでください。解凍した後、

cd SCIMES-master
python setup.py install

を実行します。
astropy_helpers がないと言われたら、

pip install astropy_helpers

を実行してください。

使用方法

必要なものを import し、fitsを読み込みます。fits は、上の SCIMES-master に入っている orion_12CO.fits を用います。

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 を実行します。

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 を定義します。そしてカタログを作ります。

metadata = {}
metadata['data_unit'] = u.Jy # u.K はなぜか対応していないので、とりあえず Jy にしておく。
cat = ppv_catalog(d, metadata)

以下でクラスタリングを実行します。

dclust = SpectralCloudstering(d, cat, hdu.header)

デフォルトでは、"luminosity"(flux) と "volume" が使われるようです。どちらか片方にしたい場合などは criteria で指定します。

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 でイテレーション回数を指定できます (大抵思うようにいきませんが)。

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

dclust.showdendro()

ダウンロード (4).png
同定したクラスタの番号と、ランダムに割り当てた色を定義します。

clusts = dclust.clusters
colors = dclust.colors
print(clusts)
print(colors)
#[11, 37, 43, ..., 428, 475, 755]
#['#1BA58D', '#5F7FE6', '#BDB864', ..., '#EAED72', '#77E1E4', '#B37499']

aplpy で plot します。

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

clusts_structure = [d[c] for c in clusts]

と定義してしまえば、天文データ解析入門 その6 (astrodendroの使い方: 基本編)天文データ解析入門 その13 (astrodendroの使い方: 中級編) の方法で色々解析ができます。

以上です。

リンク
目次

2
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
1