銀河

近傍銀河「M100 CO 分子輝線データ」

ここでは作成されたイメージデータをもとにサイエンスを行うための例として、幾つかの典型的な解析を紹介します。


Jupyter Notebook 版スクリプト
ここで紹介する解析について Jupyter Notebook 形式の Python スクリプトも用意しました。 Jupyter Notebook 形式のファイルは以下からダウンロードできます。
ipynb ファイル (2024/02/09 改訂)
Jupyter Notebook 版スクリプトを実行するためには Jupyter Notebook がインストールされている必要があります。 また、Python への CASA モジュールのインストールも必要になります (tarball でインストールされた CASA には対応しません)。 Jupyter Notebook 版スクリプトは例えば Google Colaboratory (※) 上で実行する事も出来ます (Google Colab 環境かを識別して 自動的に CASA モジュールのインストールを行います)。
なお Jupyter Notebook 版スクリプトでは、以降で説明する解析とタスク viewer を使用しない等の点で異なります。
※ Google Colaboratory (略称: Colab) は、Google の提供するブラウザから Python を記述、実行できるサービスです。 Jupyter Notebook 形式で書かれたスクリプトを読み込み、実行することが出来ます。

以降の解析は CASA 6.4 (Python 3) で動作確認しています。 使用している CASA 関係のコマンドは CASA 5 (Python 2) でも問題ないはずですが、一部の Python コマンドが Python 2 で問題が起こるかもしれません。

イメージデータ

ここで使用するデータは、ALMA Band 3 で観測された近傍渦巻銀河「M100」の CO J=1-0 分子輝線キューブイメージデータです。
Science Verification Data の一つとして取得されたデータで (「4. M100 Band 3」を参照)、12 m アレイ, 7 m アレイ, TP (単一鏡) データを結合するデモンストレーションとして作成されたイメージを使用します。 また、初期解析及びイメージングについては CASA ガイドが参考になります (CASA Guide: M100 Band3)。

JVO のサイトからも イメージ FITS ALMA00000188_00_00_00.fits を直接ダウンロードできます。

FITS イメージ表示

タスク viewer を使用します。 
# In CASA
viewer
または、 
# In command line
% casaviewer &
で、viewer が起動します。

"Viewer Display Panel" の他に、"Data Manager" というウインドウが開きます。 ここで開きたいデータを選択します。 CASA 形式のイメージデータの他に、イメージ FITS ファイルも直接開くことが出来ます。

"Viewer Display Panel" に、選択したイメージデータの画像が表示されます。

ここで扱うイメージデータはキューブデータですので、周波数 (速度) 方向の情報を持ちます。 周波数方向のチャンネルを変える場合は以下の様に操作してください。
  • チャンネル移動 (ムービー)
  • チャンネル移動 (1ch のみ)
  • ムービーの停止
  • 最小・最大チャンネルへの移動
スライダー をマウスでドラッグしてもチャンネルを変えることが出来ます。

"Viewer Display Panel" 上部のパネルの各アイコンをクリックすると各種操作ができます。
  • 拡大・縮小
  • 移動 (画面上でマウスをドラッグ)
  • カラートーン変更 (画面上でマウスをドラッグ)


右図は一例として、28 チャンネル目 (1540 km/s) を表示させた M100 の CO 輝線イメージです。

viewer の扱い方については、CASA docs の "Viewing Images" のページも参考にして下さい。

イメージ FITS 読み込み

# In CASA
orgimfits='ALMA00000188_00_00_00.fits'
orgimgdat='M100_combine_CO_cube.image'
importfits(fitsimage=orgimfits, 
           imagename=orgimgdat, overwrite=True)
タスク importfits でイメージ FITS を CASA 形式の画像ファイルに変換しています。 fitsimage には入力イメージ FITS 名を、imagename には出力される CASA 形式のファイル名 (実体はディレクトリ) を指定します。

ヘッダーから合成ビームサイズ等を確認

# In CASA
orgimgdat='M100_combine_CO_cube.image'
imhd_A_res = imhead(orgimgdat)
bmaj  =imhd_A_res['restoringbeam']['major']['value']
bmajut=imhd_A_res['restoringbeam']['major']['unit']
bmin  =imhd_A_res['restoringbeam']['minor']['value']
bminut=imhd_A_res['restoringbeam']['minor']['unit']
bmpa  =imhd_A_res['restoringbeam']['positionangle']['value']
bmpaut=imhd_A_res['restoringbeam']['positionangle']['unit']
if bmpaut == 'rad':
    bmpav=bmpa*180/math.pi; bmpautv='deg'
else:
    bmpav=bmpa; bmpautv=bmpaut
print('Synthesized Beam:')
print(f'   Major Axis = {bmaj:.7f} {bmajut}')
print(f'   Minor Axis = {bmin:.7f} {bminut}')
print(f'   Pos. Angle = {bmpa:.7f} {bmpaut}')

Synthesized Beam:
   Major Axis = 3.8690665 arcsec
   Minor Axis = 2.5323896 arcsec
   Pos. Angle = -89.5059891 deg
イメージデータのヘッダから合成ビームサイズが確認できます。 タスク imhead を使用し、Python の辞書オブジェクトとして様々な情報が読み込まれます。 これにより、ビームサイズの他に、イメージデータの各軸の座標情報・周波数情報なども確認できます。 

# In CASA
orgimgdat='M100_combine_CO_cube.image'
imhd_B_res = imhead(orgimgdat, mode='get', hdkey='restfreq')
rfrq  =imhd_B_res['value']
rfrqut=imhd_B_res['unit']
if rfrqut == 'Hz':
    rfrqv=rfrq/10**9; rfrqutv='GHz'
else:
    rfrqv=rfrq; rfrqutv=rfrqut
print(f'Rest Frequency = {rfrqv:.9f} {rfrqutv}')

Rest Frequency = 115.271201800 GHz
ここでは観測分子輝線の静止周波数を確認しています。 対象分子輝線の静止周波数は後で使用するので (ガス質量算出時)、期待される周波数かどうか注意して下さい。

サブイメージ作成

# In CASA
orgimgdat='M100_combine_CO_cube.image'
myimages ='M100_combine_CO_cube.subim'
imsubimage(imagename=orgimgdat, outfile=imname, dropdeg=True,
           region="centerbox[[413pix,357pix], [500pix,570pix]]")
イメージの余白部分を除くために、タスク imsubimage で中央部を切り出します。 領域 "region" で、切り出す四角形の範囲の中心位置及び縦横の幅を pixel で指定しています (region による指定方法は後述)。

切り出されたイメージについて、参考に 28 チャンネル目 (1540 km/s) を表示させたイメージです。

平均スペクトルを表示

タスク viewer を使って、キューブデータから平均スペクトルを作成することが出来ます。
viewer でキューブデータを開いた後、 アイコンのどれかを使って見たい方向に領域 (Region) を指定し、上部パネルの アイコンをクリックして "Spectral Profile Tool" を開いて下さい。 縦軸は Mean, Median, Sum, Flux Density、横軸は channel, radio velocity [km/s], frequency [GHz] などから選んで変更できます。

中心部の 300pix × 300pix の四角形を "region" としてしています。この範囲で平均されたスペクトルを表示します。

上記で指定された "region" 内の平均スペクトルです。ダブルピークの輝線が検出されていることが分かります。

ラインのある範囲とない範囲がどのあたりになるかに注意してください。

イメージ RMS (及び輝線ピーク) を確認

# In CASA
imname = 'M100_combine_CO_cube.subim'
imhd_A_res = imhead(imname)
xspwd = 200 ; yspwd = 250 # pixels
cetpx = [str(int((imhd_A_res['shape'][0]-1)/2)), str(int((imhd_A_res['shape'][1]-1)/2))]
ofpx = 0 ; wdpx = 5
nchan = imhd_A_res['shape'][2]
stats1=imstat(imagename=imname, 
      region='centerbox[['+cetpx[0]+'pix,'+cetpx[1]+'pix],['+str(xspwd)+'pix,'+str(yspwd)+'pix]],\
      range=['+str(ofpx)+'pix,'+str(ofpx+wdpx)+'pix]')
stats2=imstat(imagename=imname,
      region='centerbox[['+cetpx[0]+'pix,'+cetpx[1]+'pix],['+str(xspwd)+'pix,'+str(yspwd)+'pix]],\
      range=['+str(nchan-1-ofpx-wdpx)+'pix,'+str(nchan-1-ofpx)+'pix]')
stats3=imstat(imagename=imname,
      region='centerbox[['+cetpx[0]+'pix,'+cetpx[1]+'pix],['+str(xspwd)+'pix,'+str(yspwd)+'pix]],\
      range=['+str(ofpx+wdpx)+'pix,'+str(nchan-1-ofpx-wdpx)+'pix]')
imgrms=(stats1['rms'][0]+stats2['rms'][0])/2.0
imgpek=stats3['max'][0]
print(f'RMS of {imname}: {imgrms*1000:.9f} [mJy/beam]')
print(f'Peak of {imname}: {imgpek*1000:.9f} [mJy/beam]')

RMS of M100_combine_CO_cube.subim: 11.701785736 [mJy/beam]
Peak of M100_combine_CO_cube.subim: 673.587620258 [mJy/beam]
タスク imstat で、パラメータ region で指定された範囲内のイメージ RMS やピーク値などの統計量を計算しています。
領域 "region" として、
  • "centerbox" は中心位置と縦横方向の長さで四角形の領域を指定します。
  • "range" は周波数 (速度) チャンネルの範囲 (始点, 終点) を指定します。
上記の場合、"centerbox" でイメージ中心 (ヘッダーから参照) から 200pix × 250pix の四角形領域を、"range" で速度 (周波数) 方向の範囲 (RMS は両バンド端 5pix、Peak はバンド端 5pix を除いたその内側) を指定しています。

region の指定方法について、詳しくは CASA docs の "Region File Format" をご覧ください。

イメージ RMS (ノイズの大きさ) を測るときは、天体からの放射のない領域で測定する様にします。 上記では、range で輝線の検出されていない両端の周波数チャンネルを指定して RMS を求めています。

求められた統計量は Python の辞書オブジェクト形式として読み込まれ、RMS の他にフラックス密度の最大値と最小値 (及びその位置)・合計・平均なども知る事が出来ます。
全ての統計量を表示したい場合は、先のコマンド実行後に 
# In CASA
print(stats3)
等、としてください。

これらの統計量は viewer 上でも、領域 "region" の指定に制限がありますが、簡単に確認することが出来ます。

まず、調べたい周波数 (速度) チャンネルを選んで表示します。 イメージ RMS を調べるのなら、輝線のないチャンネルで測定するのが良いでしょう。 四角形領域選択ボタン ( の左側) をクリックし、画像上でマスクのドラッグによって region を指定できます。 詳細な領域設定は右側の "Regions" から "Properties" タブを開いて指定することもできます (中心位置・X 方向 Y 方向の幅を指定)。

region を指定した後に、 "Regions" の "Statistics" タブを開くと、選択した領域内の RMS を含む統計量を確認することが出来ます。

天体の輝線が検出されている領域を取り囲む region を指定すれば、ピークのフラックス密度などが分かります (右図は 49 チャンネルの場合)。

だたし、viewer で統計量を求める場合、周波数 (速度) チャンネルは表示されている 1 チャンネルのみに対応し、領域設定は四角形・楕円形 (いずれも回転はなし) 領域と、ポリゴンでの領域指定のみが可能です。

モーメントマップ作成

スペクトル線データから作成されるイメージは、2 つの空間軸と 1 つの周波数 (又は速度) 軸の 3 次元のイメージキューブになり、これらの 3 つのイメージ平面を操作してモーメントマップが作成できます。 

# In CASA
imname = 'M100_combine_CO_cube.subim'
imname_m0 = imname+'.m0'
imname_m1 = imname+'.m1'
imname_m2 = imname+'.m2'
imname_m8 = imname+'.m8'
imgrms = 11.701785736/10**3   # [Jy/beam] : RMS
cutrms = 5.0     # intensity cut off for moment1 & moment2
mochs = [8, 62]  # channel ranges for moment maps
os.system('rm -rf '+imname+'.m*')
# moment0 & moment8
immoments(imagename = imname, outfile = imname_m0, moments = 0,
       axis = 'spectral', chans=str(mochs[0])+'~'+str(mochs[1]))
immoments(imagename = imname, outfile = imname_m8, moments = 8,
       axis = 'spectral', chans=str(mochs[0])+'~'+str(mochs[1]))
# moment1 & moment2
immoments(imagename = imname, outfile = imname_m1, moments = 1,
       axis = 'spectral', chans=str(mochs[0])+'~'+str(mochs[1]),
       includepix=[imgrms*cutrms, 9999.])
immoments(imagename = imname, outfile = imname_m2, moments = 2,
       axis = 'spectral', chans=str(mochs[0])+'~'+str(mochs[1]),
       includepix=[imgrms*cutrms, 9999.])
ここで作成したモーメントマップは、moment 0 は積分強度マップ、moment 1 は速度場マップ、moment 2 は速度分散マップ、moment 8 はピーク強度マップになります。
上記では、"mochs" でモーメントマップを作成する際に使用する速度 (周波数) チャンネルを設定しています (タスク immoments のパラメータ chans で範囲指定)。
また、モーメント 1 及びモーメント 2 マップについてはノイズの影響を避けるため 5 σ ("cutrms" で指定) 以上のみのデータを使って作成しています。 "imgrms" で先に測定したイメージ RMS を指定して下さい。

作成されたモーメントマップを viewer で表示するとこのようになります。 モーメントマップを複数開きアニメーション表示することが出来ます。

モーメントマップには、他に以下の様な種類があります。 それぞれタスク immoments でパラメータ moments に以下の値を指定すれば作ることが出来ます。
  • moments=-1 - mean value of the spectrum
  • moments=0 - integrated value of the spectrum
  • moments=1 - intensity weighted coordinate; traditionally used to get 'velocity fields'
  • moments=2 - intensity weighted dispersion of the coordinate; traditionally used to get 'velocity dispersion'
  • moments=3 - median of I
  • moments=4 - median coordinate
  • moments=5 - standard deviation about the mean of the spectrum
  • moments=6 - root mean square of the spectrum
  • moments=7 - absolute mean deviation of the spectrum
  • moments=8 - maximum value of the spectrum
  • moments=9 - coordinate of the maximum value of the spectrum
  • moments=10 - minimum value of the spectrum
  • moments=11 - coordinate of the minimum value of the spectrum
モーメントマップについて、詳しくは以下を参考にして下さい。

プロファイルマップ作成

プロファイルマップ (イメージ領域を格子状に分割し、それぞれの位置での平均スペクトルを表示) を作成します。 
# In CASA
imname = 'M100_combine_CO_cube.subim'
plotprofilemap(imagename=imname, numpanels='6,7', spectralaxis='velocity',
               showtick=True, showticklabel=True, showaxislabel=True,
               figfile=imname+'.profilemap.png', overwrite=True)
パラメータ numpanels で縦横方向に何分割するか指定できます。ここでは RA (横) 方向に 6 分割、Dec (縦) 方向に 7 分割しています。

作成されたプロファイルマップは、別ウインドウが開いて表示されるとともに、figfile で指定された名前の PNG ファイルにも出力されます。

PV マップ作成

イメージキューブデータは、空間方向 XY の 2 成分と周波数 (速度) 成分の 3 次元の情報を含みます。 通常、viewer 等を使ってキューブデータを表示させると空間方向 2 成分が表示されます。 一方、任意の空間方向と周波数 (速度) 方向の強度分布を、位置速度図 (PV マップ; Position-Velocity map) と言います。 銀河の中心を基準にカットした PV マップから、銀河のガス ダイナミクスが分かります。

PV マップは、viewer を使って以下の様に簡単に作成できます。
キューブデータを表示した後、"P/V ツール" ( アイコン) をクリックし、PV マップを作成したい方向に始点から終点までマウスをドラッグして下さい。 右下の "Region" に "pV" タブが開きますので、"Generate P/V" をクリックすると、新しいウインドウが開き PV マップが表示されます。 また、この "pV" タブで、PV マップを作るためのパラメータを任意に指定することも出来ます。 例えば "averaging width" では、スライスする方向と垂直な方向へ平均する幅をピクセル単位で指定できます。 この値を大きくすると、PV マップの SN が良くなるはずです。

銀河の長軸に近いと思われる方向にスライスした PV マップを作成します。 右図では、カラーでキューブイメージを表示 (チャンネル 28) した上で、参考にコントアでモーメント 0 マップを重ねて表示しています。
赤線が PV マップを作成するためにスライスする位置を示しています。

新しいウインドウが開き、作成された PV マップが表示されます。 なお、周波数 (速度) 方向にはキューブの全範囲が指定されて作成されます。

拡大するとこのように見えます。 横軸が位置 (オフセット)、縦軸が速度になります。

また、PV マップをタスク impv を使って作成すると以下の様になります。 
# In CASA
imname = 'M100_combine_CO_cube.subim'

# search peak position in moment 0
imname_m0 = imname+'.m0'
imhd_A_res = imhead(imname_m0)
xspwd = 400 ; yspwd = 400 # pixels
cetpx = [str(int(imhd_A_res['shape'][0]/2)), str(int(imhd_A_res['shape'][1]/2))]
stats4=imstat(imagename=imname_m0, 
     region='centerbox[['+cetpx[0]+'pix,'+cetpx[1]+'pix],['+str(xspwd)+'pix,'+str(yspwd)+'pix]]')
pakpos=stats4['maxpos'][0:2]

# pv maps
pvname1 = imname+'.pv01'
ctx=pakpos[0] ; cty=pakpos[1] # peak pixel of moment 0
pvpa  =  -30  # PA (deg)
pleng =  200  # arcsec
vlrag = '0~69' # vel.ch.range [ch]
os.system('rm -rf '+pvname1+'*')
impv(imagename = imname, outfile = pvname1,
     mode = 'length', center = [ctx, cty],
     length = str(pleng)+'arcsec', pa = str(pvpa)+'deg',
     width = 7, chans = vlrag)
上記では、まずタスク imstat でモーメント 0 マップのピークのピクセル位置を探します。 そのピーク位置を基準に、タスク impv で指定した位置角 (PA) "pvpa (deg)" の方向へスライスして、長さ "pleng (arcsec)" の PV マップを作成します。 また、速度方向へは "vlrag (channels)" の範囲で作成するように制限出来ます。 これで、M100_combine_CO_cube.subim.pv01 という PV マップが作成されます。

右図はここで作成された PV マップです (詳細は Jupyter Notebook 版スクリプトを参照)。

総分子ガス質量の算出

CO 分子輝線強度から総分子ガス質量を推定します。 

# In Python
imname = 'M100_combine_CO_cube.subim'
imname2 = imname+'.thrd'
rfrq = 115.271201800*10**9   # [Hz] : rest frequency for CO line
imgrms = 11.701785736/10**3  # [Jy/beam] : RMS
nsigma = 5.0   # σ threshold
threshold = imgrms*nsigma
os.system('rm -rf '+imname2)
immath(imagename=[imname], outfile=imname2,
       expr='IM0[IM0>%f]'%threshold)
imhd2_res = imhead(imname2)

xspwd = 400 ; yspwd = 480 # pixels
cetpx = [str(int((imhd2_res['shape'][0]-1)/2)), str(int((imhd2_res['shape'][1]-1)/2))]
ofpx = 0 ; wdpx = 5
nchan = imhd2_res['shape'][2]
stats=imstat(imagename=imname2, 
      region='centerbox[['+cetpx[0]+'pix,'+cetpx[1]+'pix],['+str(xspwd)+'pix,'+str(yspwd)+'pix]],\
              range=['+str(ofpx+wdpx)+'pix,'+str(nchan-1-ofpx-wdpx)+'pix]')
imgflx=stats['flux'][0]
print(f'Total flux = {imgflx:.7f} Jy km/s')

cent_freq=(imhd2_res['refval'][2]+((imhd2_res['shape'][2]/2)-imhd2_res['refpix'][2])*imhd2_res['incr'][2])
sysfreq=cent_freq
zzz = rfrq/sysfreq-1
sysfreq=sysfreq*10**(-9)
spunit='GHz' if imhd2_res['axisunits'][2]=='Hz' else 'no match'
print(f'center freq = {sysfreq:.7f} {spunit}')
print(f'redshift    = {zzz:.7f}')

ddd=15.2   # [Mpc] : luminosity distance for M100 (Sun et al.2018)
alpha=4.35 # [Mo pc-2 (K km/s)-1] (Sun et al.2018)
luminosity=(3.25*10**7)*imgflx*sysfreq**(-2)*ddd**2*(1+zzz)**(-3) # (Solomon & Vanden Bout 2005)
print(f'CO luminosity = {luminosity:.7e} K kms-1 pc2')
mmass=alpha*luminosity
print(f'Total molecular gass mass = {mmass:.7e} Mo')

Total flux = 1358.9713282 Jy km/s
center freq = 114.6656090 GHz
redshift    = 0.0052814
CO luminosity = 7.6392645e+08 K kms-1 pc2
Total molecular gass mass = 3.3230801e+09 Mo
キューブイメージから閾値 5 σ ("nsigma") 以下のピクセルを除外し、指定した範囲内のピクセル強度からトータルフラックスを算出します。 "imgrms" で先に測定したイメージ RMS を指定して下さい。

右図は閾値以下を除外したキューブから作成した積分強度図 (モーメント 0 マップ) で、赤い四角形で囲った範囲でトータルフラックスを求めています。


CO 光度 (luminosity) を求めるにあたり、バンド中心の周波数及び CO J=1-0 輝線の静止周波数 (前述) から赤方偏移 (redshift) も計算しています。

分子ガス H2 質量 M(H2) は、CO 分子輝線積分強度 (total flux) SCO Δv より、以下の式で示されるとして計算しています (Solomon & Vanden Bout 2005)。
なお、それぞれの単位は (total CO luminosity), (integrated intensity), (observing frequency), (luminosity distance), (H2 molecular mass) になります。
また、変換係数 (CO-to-H2 conversion factor) αCOは
を採用しています (Sun et al.2018)。

チャンネルマップ

チャンネルマップとは、一定の視線速度の範囲 (チャンネル) ごとの輝線スペクトルの平均強度の分布図で、天体の速度構造を把握するのに有用です。

CASA には現状チャンネルマップを直接きれいに作るタスクがないので、Python で実現するために別途モジュールなどが必要になります。
Jupyter Notebook 版スクリプトに参考例を記述しています。 このスクリプトでは Kapteyn Package を使用してチャンネルマップを作成しています。 Kapteyn Package でのチャンネルマップについては、Kapteyn Package のホームページの "Tutorials" の "Mosaics of plots" を参考にして下さい (サンプルスクリプト)。 Kapteyn Package には他にも様々な機能があります。

チャンネルマップを作成する際、速度 (周波数) 方向のチャンネル (pixel) 数を減らしたい場合があります。 タスク imrebin でキューブデータを速度方向に平滑化する事が出来ます。 以下では、3 軸目の速度 (周波数) 方向を 2 チャンネル毎に平均化して出力しています。 
# In CASA
myimages2='M100_combine_CO_cube.subim'
myimages3='M100_combine_CO_cube.subim.bin'
os.system('rm -rf '+myimages3)
imrebin(imagename=myimages2, outfile=myimages3, factor=[1,1,2])
exportfits(imagename=imname, fitsimage=imname+'.fits', 
           dropstokes=True, overwrite=True)

右図は、Jupyter Notebook 版スクリプトで作成されたチャンネルマップの出力例です。
元の M100_combine_CO_cube.subim というキューブデータを 2 チャンネル毎に平均化して、12 ch～60 ch に相当する範囲 (平均化後は 6 ch～30 ch) を、1 チャンネル間隔で 25 枚のチャンネルイメージを表示しています。


ページ先頭に戻る

Last Update: 2024.02.15