Astropyで天文データ解析 [2] -convolution, filtering-
はじめに
本記事ではPythonを用いた天文解析の理解を深めることを目的に、Astropy公式チュートリアルのSynthetic Images from simulated dataとConvolution and Filteringを参考に記事にまとめる。なお、筆者の備忘録として活用できるように分かりやすい解説を付け加えるよう配慮した。
本記事のゴール
-
astropy.wcsを用いた天体画像へのWCS座標の割り当て -
astropy.modeling.modelを用いたPSF(Point Spread Function)の構築 -
astropy.convolutionを用いた生データとPSFの畳み込み - ストークスパラメータ
I Q U - 天体画像にquiverを重ねる
必要なライブラリ
Astropy関係で必要なライブラリは以下の通り。
| 必要なライブラリ | どんなライブラリ？ |
|---|---|
astropy.utils.data.download_file |
fitsデータのDL |
astropy.io.fits |
HDU(Header Data Unit)の取得 |
astropy.units |
単位取得 |
astropy.coordinates.SkyCoord |
天体座標の変換や取得 |
astropy.wcs.WCS |
WCS取得 |
astropy.convolution.Gaussian2DKernel |
画像を2Dガウシアンで畳み込む |
astropy.convolution.convolve_fft |
単位の取得 |
FITSファイルの読み込みと操作
astropy.utils.dataのdownload_fileを用いて公式チュートリアルのFITSファイルを読み込む。
FITSファイルの読み込みと操作
from astropy.utils.data import download_file
from astropy.io import fits
fits_file = download_file(
'http://data.astropy.org/tutorials/synthetic-images/synchrotron_i_lobe_0700_150MHz_sm.fits',
cache=True)
hdu_list = fits.open(fits_file)
.info()でFITSファイルの情報を取得できる。
hdu_list.info()
出力結果
Filename: \hoge\hogehoge\.astropy\cache\download\url\8da27de5aa6b0db633441e82715bedf3\contents
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 4 ()
1 NN_EMISSIVITY_I_LOBE_150.0MHZ 1 ImageHDU 23 (1024, 1452) float64
出力結果を見ると、PrimaryHDU以外にExtensionHDUのNN_EMISSIVITY_I_LOBE_150.0MHZが存在することがわかる。これらのHDUはhdulist[0], hdulist[1]のようにindexで取得することも可能だが、直接Nameを指定しても取得できる。
hdu = hdulist['NN_EMISSIVITY_I_LOBE_150.0MHZ']
print(hdu.header)
出力結果
XTENSION= 'IMAGE ' / Image extension
BITPIX = -64 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 1024
NAXIS2 = 1452
PCOUNT = 0 / number of parameters
GCOUNT = 1 / number of groups
EXTNAME = 'NN_EMISSIVITY_I_LOBE_150.0MHZ' / extension name
BTYPE = 'nn_emissivity_i_lobe_150.0MHz'
BUNIT = 'Jy/arcsec**2'
WCSAXES = 2
CRPIX1 = 512.0
CRPIX2 = 726.0
CDELT1 = 9.42382812499999E+19
CDELT2 = 9.42382812499999E+19
CUNIT1 = 'cm '
CUNIT2 = 'cm '
CTYPE1 = 'LINEAR '
CTYPE2 = 'LINEAR '
CRVAL1 = 0.0
CRVAL2 = 0.0
LATPOLE = 90.0
WCSNAME = 'yt '
読み込んだFITSキーワードのCUNIT1, CUNIT2を見るとWCSNAMEキーワードのytで指定された座標系の単位は、cmであることがわかる。これをRADEC(赤道座標)に変換したい。なお、各Keywordの詳細はFITSの手引きを参照。
まずは、FITSファイルのヒストグラムを取得してみる。
HDUからヒストグラムを取得
matplotlibとNumPyを用いてヒストグラムを描く。
import matplotlib.pyplot as plt
import numpy as np
#最大値と最小値を取得
print(f'最大値: {hdu.data.max()}')
print(f'最小値: {hdu.data.min()}')
出力結果
最大値: 129.7177858088622
最小値: 0.0
ヒストグラムを表示するポイントは、2次元配列データを１次元配列データに平坦化する(flatten)こと。平坦化はNumPyのflatten()やravel()を用いればよい。hdu.dataの型は、
print(type(hdu.data))
出力結果
<class 'numpy.ndarray'>
見ての通りndarrayオブジェクトであるため、HDUのデータにはflatten()やravel()を用いた平坦化が有効。これらを考慮してヒストグラムを描くと以下のようになる。
#真数条件のWarningを無視
np.seterr(divide='ignore')
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(np.log10(hdu.data.flatten()), range=(-3, 2), bins=100)
ax.set_xlabel('Pixel values (logarithmic values)')
ax.set_ylabel('Frequency')
出力結果
ヒストグラムを見ると、-1～1付近の頻度が高く、データが集中していることがわかる。
ヒストグラムによりデータ範囲が取得できれば、適切な範囲でFITS画像を可視化することが可能。
適切なデータ範囲でFITS画像を表示する
対数表記でプロットするとともに、非数を回避するために小さな値(1E-3)を加えておく。また、ヒストグラムを見ると、-1～1付近にデータが集中していることがわかるため、vmin=-1, vmax=1として適切な範囲でデータをプロットする。
from mpl_toolkits.axes_grid1 import make_axes_locatable
divider = make_axes_locatable(ax)
cax = divider.new_horizontal(size="2%", pad=0.05)
fig.add_axes(cax)
fig = plt.figure(figsize=(6,12))
ax = fig.add_subplot(111)
# 対数表記および非数を回避するために、np.log10(hdu.data+1E-3)とする
ax.imshow(np.log10(hdu.data+1E-3), cmap='jet', vmin=-1, vmax=1, origin='lower')
出力結果
なお、cmap='jet'とした。
物理座標を各座標系に変換
公式チュートリアルのFITSヘッダを見ると、縦横軸ともに単位がcmになっており、物理座標になっている。しかし、大抵観測データは赤道座標か銀河座標だから、これをastropyを用いて変換する。当然、物理座標を赤道座標に変換するためには、観測天体の距離や中心座標が必要。
観測天体の距離や中心座標を決める
まずは観測天体の距離と中心座標を設定。
# 天体までの距離
dist_obj = 200*u.Mpc
# 赤経をhh:mm:ss, 赤緯をdd:mm:ss形式で指定。あとで度表記になおす
ra_obj = '19h59m28.3566s'
dec_obj = '+40d44m02.096s'
ピクセル値をdist_objで割って、cmをdegreeに変換する。なお、ピクセル値はCDELTヘッダで取得する。このヘッダは参照点の位置の増分を表す。すなわち、画像の解像度(分解能)に対応するもの。
import astropy.units as u
cdelt1 = ((hdu.header['CDELT1']*u.cm/dist_obj.to('cm'))*u.rad).to('deg')
cdelt2 = ((hdu.header['CDELT2']*u.cm/dist_obj.to('cm'))*u.rad).to('deg')
print(f'CDELT1:{cdelt1} CDELT2:{cdelt2}')
出力結果
CDELT1:8.74922223983974e-06 deg CDELT2:8.74922223983974e-06 deg
では、次にastropy.wcsの.WCS()を用いて赤道座標用のFITSヘッダーを準備する。
赤道座標用のFITSヘッダを作成し天体画像を表示する
変数wに赤道座標用のFITSヘッダを格納する。
# データ配列の座標軸の本数を決める(２次元なら画像、３次元ならdata cube)
w = WCS(naxis=2)
# 参照点のピクセル座標
w.wcs.crpix = [hdu.data.shape[0]/2,hdu.data.shape[1]/2]
# 画素の大きさ(単位: 度)
w.wcs.cdelt = [-cdelt1.base, cdelt2.base]
# 度をRA.とDEC.に変換
c = SkyCoord(ra_obj, dec_obj)
w.wcs.crval = [c.ra.deg, c.dec.deg]
# 軸ラベルの単位を度に設定
w.wcs.cunit = ['deg', 'deg']
作成したwを既存のヘッダに反映する。
# to_header()でヘッダを生成
wcs_header = w.to_header()
# updateメソッドでhduヘッダを更新
hdu.header.update(wcs_header)
print(hdu.header)
出力結果
XTENSION= 'IMAGE ' / Image extension
BITPIX = -64 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 1024
NAXIS2 = 1452
PCOUNT = 0 / number of parameters
GCOUNT = 1 / number of groups
EXTNAME = 'NN_EMISSIVITY_I_LOBE_150.0MHZ' / extension name
BTYPE = 'nn_emissivity_i_lobe_150.0MHz'
BUNIT = 'Jy/arcsec**2'
WCSAXES = 2 / Number of coordinate axes
CRPIX1 = 726.0 / Pixel coordinate of reference point
CRPIX2 = 512.0 / Pixel coordinate of reference point
CDELT1 = -8.7492222398396E-06 / [deg] Coordinate increment at reference point
CDELT2 = 8.74922223983969E-06 / [deg] Coordinate increment at reference point
CUNIT1 = 'deg ' / Units of coordinate increment and value
CUNIT2 = 'deg ' / Units of coordinate increment and value
CTYPE1 = 'LINEAR '
CTYPE2 = 'LINEAR '
CRVAL1 = 299.8681525 / [deg] Coordinate value at reference point
CRVAL2 = 40.733915555556 / [deg] Coordinate value at reference point
LATPOLE = 90.0 / [deg] Native latitude of celestial pole
WCSNAME = 'yt '
MJDREF = 0.0 / [d] MJD of fiducial time
あとは、アップデート後のヘッダをprojectionに指定し、座標軸を変換すればよい。
wcs = WCS(hdu.header)
fig = plt.figure(figsize=(6,12))
ax = fig.add_subplot(111, projection=wcs)
ax.imshow(np.log10(hdu.data+1e-3), vmin=-1, vmax=1, origin='lower')
ax.set_xlabel('RA')
ax.set_ylabel('Dec')
出力結果
点像分布関数を作成
ここではastropy.convolutionのGaussian2DKernelを使用してPoint Spread Function (PSF(点像分布関数))を作成する。望遠鏡のresolutionを設定して、CDELTヘッダを使ってピクセル単位の
PSFを作成
# 望遠鏡の分解能を設定。ここでは１秒角とする
telescope_resolution = 1*u.arcsecond
# ピクセル単位のσの計算。.to('deg')でresolutionをdeg.に変換しCDELT2で割る
sigma = telescope_resolution.to('deg')/cdelt2
Gaussian2DKernelの引数にsigmaを渡してmatplotlibで表示する。
# PSFの作成
psf = Gaussian2DKernel(sigma)
fig = plt.figure()
ax = ax.add_subplot()
ax.imshow(psf.array.value)
出力結果
Discussion