【以下追記】
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()
同定したクラスタの番号と、ランダムに割り当てた色を定義します。
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()
clusts_structure = [d[c] for c in clusts]
と定義してしまえば、天文データ解析入門 その6 (astrodendroの使い方: 基本編) や 天文データ解析入門 その13 (astrodendroの使い方: 中級編) の方法で色々解析ができます。
以上です。
リンク
目次