at 国立天文台

Python+Jupyter notebookによる光赤外撮像データ処理入門
2018-08-30, 31 中島 康

Jupyter notebookの使い方
Jupyter notebookの起動
新規ノートブックファイルの作成
Pythonとの対話
対話を保存
Jupyter notebookの終了

メモの書き込み
Markdownで書き込み

講習0 : Pythonとnotebook
notebook形式の資料の使い方
Pythonについて少し
notebookでの変数の読み込み
ipythonの補完機能

ipython の補完機能 ~ [tab]キー
最初の何文字か + [tab]講習1 : IRAFタスクをPythonから使ってみる

候補のリスト•

定義済み変数•

使用可能な関数など

IRAFとPythonをつなぐPyRAF
IRAFの基本タスク(コマンド)
display, imexam, imstat, imarith
IRAFタスクのパラメータ設定

変数.何文字か + [tab]候補のリスト•

使用可能なメソッド

IRAFタスクから戻り値を取得

ドット0文字なら全候補

演習1

DS9にFITSを表示
ds9をあらかじめ起動しておく !
from pyraf import iraf
対象ファイルあるいは import pyraf.iraf as iraf

iraf.display( ʻxxxxxx.ﬁtsʼ, 1)
ds9での表示Frame番号(なしの場合、直前で使った番号)

imexamで画像を調べる
対象ファイルDS9で表示されている
Frame番号

imstatで画像の統計を調べる

iraf.imexam( ʻxxxxxx.ﬁtsʼ, 1)

対象ファイルm : 5x5ピクセルのカウント値の統計

iraf.imstat( ʻxxxxxx.ﬁtsʼ)a : 星の特徴量(ピーク値、FWHMなど)
r : 星のradial proﬁleq : 終了

irafタスクのパラメータ設定
iraf.epar( ʻimstatʼ)

パラメータのリセット
1. GUIで設定

iraf.imstat( ʻxxxxxx.ﬁtsʼ, ﬁelds=ʻmean, stddevʼ)

iraf.imstat.unlearn( )またはiraf.unlearn( ʻimstatʼ )2. 関数の引数として書き込む
iraf.imstat.ﬁelds = ʻmidpt, mean, stddevʼ
iraf.imstat.lower = 0
3. 変数に代入

戻り値を変数へ

画像の四則演算2w
out = iraf.imstat( ʻxxxxxx.ﬁtsʼ,Stdout=1)
標準出力へ出る内容が変数に(文字列として)数値 または ファイル
対象ファイル2w ʻ-ʻ, 100, ʻout.ﬁtsʼ)iraf.imarith( ʻxxxxxx.ﬁtsʼ,format=ʻnoʼ ヘッダ行を非表示
演算子

出力ファイル

IRAFが勝手にFITS拡張するのを禁止する

helpドキュメントを読む

iraf.reset( fkinit = ʻappend = noʼ )iraf.help( ʻimstatʼ )講習2 : IRAFで1次処理

演習1
data1/SUPA00317885.ﬁtsはtarget2の生データです。

生データ取得 + 一次処理に必要なデータ取得

新しいノートブックファイルを作成し、天体以外の信号の除去

1. imexamでバックグラウンドの値とばらつき、星の特徴量、一次処理感度と光学系のムラの補正を調べる。
2. imstatでカウント値のmedian, mean, standard deviationを一次処理済みデータ求める。
3. imarithで、全pixelからmedian値を引く。(元のファイルに上書きせずに、別名のﬁtsファイルに)例) 点光源の測光 ( + 標準星データ取得 )

測光値の較正可視光

raw(x, y) = flat(x, y) x object(x, y) + dark(x, y) + bias(x, y)raw(x, y) = flat(x, y) x object(x, y) + dark(x, y) + bias(x, y)flat(x, y)object(x, y)

ﬂat(x,y) : ピクセル間の感度むら, 光学系の透過むら
1.10.90.9ﬂat(x,y)×const. が検出される。

1.01.20.8

近年のCCDではdarkは無視できるレベル

1.10.7

1.1

200

204

201

220

183

181

object(x,y)=const. の光がやってきたとしてもbiasは、raw(x, y)のオーバースキャン領域から推測(x方向への依存は無しと仮定し)× 198 400 202 = 198 480 162
205

199

200

226139220可視光raw(x, y) = flat(x, y) x object(x, y) + dark(x, y) + bias(x, y)
ﬂat(x,y) : ピクセル間の感度むら, 光学系の透過むら

ここがオーバーobject(x,y)=const. の光がやってきたとしてもスキャン領域

ﬂat(x,y)×const. が検出される。
近年のCCDではdarkは無視できるレベル
biasは、raw(x, y)のオーバースキャン領域から推測(x方向への依存は無しと仮定し)mean

3-sigma clip

Y方向にも大部分で一様で、端のほうで 0.数% 増加
10030

standard deviation = 1 sigma

3-sigmaフラット
raw(x, y) = flat(x, y) x object(x, y) + dark(x, y) + bias(x, y)

9990

オーバースキャンのメジアン
iraf.imstat( ʻdata1/SUPA00317505.ﬁts[2049:2080, *]ʼ)object(x, y) = constant (一様光)を観測してやれば、
ﬂat(x, y) =raw(x, y) - bias(x, y)
constant

光が当たっている部分のメジアン
iraf.imstat( ʻdata1/SUPA00317505.ﬁts[1:2048, *]ʼ)constant = median(raw(x, y) - bias(x, y))
規格化iraf.imcombine()
im01.ﬁts
1次処理

im02.ﬁts

im03.ﬁts

ドームフラットデータから ﬂat(x, y)

101 104 102

105 101

99

98

天体生データから bias(x, y)

99

100 101

100

99

102

103 102

103

97

101

98

96

105 101 103

object(x, y) =raw(x, y) - bias(x, y)flat(x, y)

im02.ﬁts

99

averageでcombine

medianでcombine

101.3 102.3 102

101 102 102

100.7 100.3 100.7

100 100 101

103

103

98.7 99.7

98

100

演習2

iraf.imcombine()
im01.ﬁts

100

102 1052-1

im03.ﬁtstarget2を観測した、'./data1/SUPA00317885.ﬁts' について、バイアス引き+フ101 104 102

105 101

9998

102 10599

100 101

10099102103 102

103

971019896

105 101 103

ラット割りの処理をしましょう。これは5番フレームです。フィルターも同じBバンドなので、フラット割りには、'bﬂatn5.ﬁts'が使えます。

10099

この結果のフレームを'btarget2n5.ﬁts'と呼ぶことにします。(後の演習で利用します)

2-2
2番フレームの生データ'./data2/SUPA00317702.ﬁts' について、バイアス引き+フラット割りの処理をしましょう。

iraf.imcombine(input=ʻim01.ﬁts, im02.ﬁts, im03.ﬁtsʼ,
output=ʻtest.ﬁtsʼ, combine=ʻmedianʼ)combine=ʻaverageʼ

先ほどの5番フレームとは違い、これは2番フレームなので、2番フレームのためのフラットを作成する必要があります。

(1) './data2/SUPA00317502.ﬁts' を規格化したものをフラットとして作成する。
(2) './data2/SUPA003175[0-6]2.ﬁts' から平均のフラットを作成する。
(3) 上のどちらかのフラットを使って、バイアス引き後のフラット割りを行う。
注意 : 2番フレームはオーバースキャン領域が5番とは異なる。

at 国立天文台

講習3 : 星の測光

2日目

Python+Jupyter notebookによる光赤外撮像データ処理入門

視野の星を検出させ、 iraf.daoﬁnd()
それらをアパーチャ測光

iraf.phot()

測光値の較正 ~ 標準星
2018-08-31 中島 康

アパーチャ測光

アパーチャ半径星を含む円内のカウントを総計する。
明るさの分かっている星のカウントと比べて、明るさFWHMを求める。

標準星の観測が必要

mean3-sigma clip

standard deviation = 1 sigma

3-sigma

iraf.phot()のパラメータ設定

skyレベルをココで計算

星の測光

演習3

daoﬁndの結果のファイル

測光対象ファイル

測光対象ファイル

iraf.phot(ʻxxxx.ﬁtsʼ, coords=ʻout1.cooʼ
iraf.daoﬁnd( ʻxxxxxx.ﬁtsʼ, output=ʻout1.cooʼ )output=ʻout1.magʼ)
検出結果のテキストファイル

iraf.photの測光結果ファイル

測光結果を書き出す(テキスト)ファイル

演習3
演習2-1で処理をしたbtarget2n5.ﬁtsで測光をしてみましょう。このときも、オーバースキャン領域などの不要な部分を削除して行いましょう。

iraf.txdump(ʻout1.magʼ,ﬁelds=ʻxc, yc, mag, merrʼ, Stdout=ʻresult1.txtʼ)btarget2n5.ﬁts'btarget2n5.ﬁts'の視野の中には測光標準星は写っていません。ただし、上のtarget1と近い時間に観測したデータですので、等級ゼ抽出する項目

講習4 : Numpyの基本
配列のベクトル演算 ~ ndarray

この結果の書き出し先ロ点は同じだと仮定し、上と同じ較正値(1.672)を使ってください。

ndarrayの作成
np.array(［リスト])

二次元配列の注意点
numpyの関数をいくつかndarrayの演算(ベクトル演算)
np.array([1, 2, 3])+np.array([10, 20, 30])

=np.array([11, 22, 33])2次元のndarray は インデックスが(y, x)
np.array([ [7, 3, 8], [13,11,16], [105, 121, 153] ])np.array([ [7, 3, 8],[13,11,16],[105, 121, 153] ])矩形部分の抽出

離散部分の抽出

iraf.imstat( ʻdata1/SUPA00317505.ﬁts[2049:2080, *]ʼ)aa = np.array([0, 10, 20, 30, 40, 50])aa[ np.array([0, 2, 5]) ]array[:, 2048:2080]

または array[:, 2048: ]array([0, 20, 50])

1. 区間ワイルドカードは :2. IRAFのインデックスは1はじまりPythonのインデックスは0はじまり

aa = np.array([0, 10, 20, 30, 40, 50])aa[ aa > 30 ]array([40, 50])

3. numpyでは (y, x)

aa = np.array([0, 10, 20, 30, 40, 50])aa[ aa > 30 ]array([40, 50])

numpy.where( )aa = np.array([0, 10, 20, 30, 40, 50])x = np.where(aa > 30)xaa = np.array([0, 10, 20, 30, 40, 50])aa[ aa > 30 ] = -99
aaarray[(0, 10, 20, 30, -99, -99)]

numpyの関数をいくつか

aa[x](array([4, 5]), )array([40, 50])

numpy.stack( )numpy.median( [スタックされたndarray], axis=0)s1

aa = np.array([1, 2, 3])bb = np.array([4, 5, 6])

stack!

cc = np.array([10, 20, 30])
s1 = np.stack((aa, bb, cc))s1

array([[1,2,3],[4, 5, 6],[10, 20, 30]])

array([[1,2,3],[4, 5, 6],[10, 20, 30]])

np.median(s1, axis=0)

array([4., 5., 6.])

np.average(s1, axis=0)

array([5., 9., 13.,])

np.sum(s1, axis=0)

numpy.append( )axis=0array([15., 27., 39.,])numpy.empty( )aa = np.array([1, 2, 3])

aa = np.array([1, 2, 3])bb = np.array([4, 5, 6])bb = np.array([4, 5, 6])

cc = np.array([10, 20, 30])

cc = np.array([10, 20, 30])

0 x 3の2次元ndarrays3 = np.empty((0, 3))s3 = np.append(s3, aa[np.newaxis, :], axis=0)
s3 = np.append(s3, bb[np.newaxis, :], axis=0)
s3 = np.append(s3, cc[np.newaxis, :], axis=0)s3
array([[1,2,3],[4, 5, 6],[10, 20, 30]])

s2 = aa[np.newaxis, : ]s2 = np.append(s2, bb[np.newaxis, :], axis=0)
s2 = np.append(s2, cc[np.newaxis, :], axis=0)s2
array([[1,2,3],[4, 5, 6],[10, 20, 30]])

numpy.loadtxt( )演習4(1) 下の二次元配列の演算をnumpyを用いて行ってください。

数値(半角)だけのテーブル(サイズ)であれば、これを使うのが楽。(ヘッダ行に文字があっても、無視できる。) forループを回さなくて良い！
numpy.savetxt( )(2) 下の3つの二次元配列のメジアンを、numpy.stack()とnumpy.median()を用いて求めてください。

講習5 : astropy.io.ﬁtsの基本
IRAFなしでFITSの処理
getdata()でFITSの読み込み
numpyを使って処理
writeto()でFITSの書き出し

モジュールのよびだしfrom astropy.io import ﬁts(画像)データの読み込み
imdata = ﬁts.getdata(ʻxxxxx.ﬁtsʼ)
ヘッダの読み込み
header = ﬁts.getheader(ʻxxxxx.ﬁtsʼ)

フラットの作成
pyraf版

データの書き出し(新規)ﬁts.writeto(ʻnewname.ﬁtsʼ, imdata, header)
データの書き出し(更新)
ﬁts.update(ʻexsiting.ﬁtsʼ, imdata, header)

numpy版

演習5

講習6 : matplotlibの基本

1. 上のmy10x10.ﬁtsを作成し、DS9で見てxとyが反転していることを確かめてください。

ヒストグラム (光度関数)

2. astropy.io.ﬁtsとnumpyを使って以下の処理をしてください。

X-Y プロット (等級 vs 等級エラー)(1) './data2/SUPA003175[0-6]2.ﬁts'から2番フレーム用のフラットを作成。
(2) './data2/SUPA00317882.ﬁts' について、バイアス引き+フラット割りの処理(3) trimmingして(2)の結果からオーバスキャン部をとりのぞく。

FITS画像の表示
FITS画像にオーバープロットモジュールのよびだし
import matplotlib.pyplot as plt

ヒストグラムplt.hist( ndarray )put.show()基本

notebook内に描画するために
データプロット%matplotlib inline
plt.scatter( ndarray, ndarray )put.show()

FITSの表示基本
plt.ﬁgure( ﬁgsize=(7, 14) )plt.imshow(imdata, plt.cm.gray,origin=ʻlowerʼ, interpolation=ʻnoneʼ)put.show()

基本

演習6
演習3で行った'btarget2n5.ﬁts'の測光結果を用いて、(1) 「光度関数のヒストグラム」と 「等級vs等級エラーのプロット」を作成してください。
(2) FITS画像をnotebookに表示して、そこに測光した星をプロットしてください。