at 国立天文台

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

Jupyter notebookの使い方
講習0 : Pythonとnotebook
講習1 : IRAFを使ってみる
講習2 : IRAFで一次処理講習3 : 星の測光
講習4 : Numpyの基本
講習5 : astropy.io.ﬁtsの基本
講習6 : matplotlibの基本
講習7 : スクリプト作成等

PyRAF(IRAF)-

Python interactive + script + ノートブックで記録に- サンプルデータ(Suprime Cam)のデータ処理、測光を行いながら、基本的なスキルを身につける。
- 今回の資料を見ながら、それを思い出して、後に自分のデータ処理、解析に応用することができるようになる。
- 講習で作成したnotebookファイルを持ち帰って、自分のマシンで復習あるいは今後の参考に (USB?クラウド?)

IRAF Image Reduction and Analysis Facility
アメリカのNOAO (国立光学天文台)が作成 (1986~ )光赤外データ処理・解析のデファクトスタンダード
画像演算、測光などのタスク(プログラム)の集合
長期間、多くの人によって使われてきた「枯れたシステム」--- バグが淘汰されている
IRAFのコマンドラインから対話的データ処理
独自のプログラム言語 IRAF-CL

Python
今、最も広く使われているコンピュータ言語のひとつ
PyRAFが「仲介」することでPythonからIRAFのタスクを利用できる！
数多くのライブラリ ̶ 数値計算、統計、可視化、機械学習、、、、、
多くの教科書、ウェブ上の情報

Python
今、最も広く使われているコンピュータ言語のひとつ
PyRAFが「仲介」することでPythonからIRAFのタスクを利用できる！
数多くのライブラリ ̶ 数値計算、統計、可視化、機械学習、、、、、
多くの教科書、ウェブ上の情報

実験ノート(研究ノート)をとろう！

Jupyter notebookがひとつのソリューション

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

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

講習0 : Pythonとnotebook
notebook形式の資料の使い方
Pythonのforループ ~ 字下げが意味をもつ
Pythonのモジュール
notebookでの変数の読み込みipython の補完機能 ~ [tab]キー
最初の何文字か + [tab]候補のリスト•

定義済み変数•

使用可能な関数など

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

0文字なら全候補

使用可能なメソッド

講習1 : IRAFを使ってみる
IRAFとPythonをつなぐPyRAF
IRAFの基本タスク(コマンド)
display, imexam, imstat
IRAFタスクのパラメータ設定
演習1

演習1
data1/SUPA00317885.ﬁtsはtarget2の生データです。
新しいノートブックファイルを作成し、1. imexamでバックグラウンドの値とばらつき、星の特徴量、を調べる。
2. imstatでカウント値のmedian, mean, standard deviationを求める。
3. imstatのnclipを2以上にしてみてください。どうなりましたか？

mean

3-sigma clip

standard deviation = 1 sigma

3-sigma講習2 : IRAFで1次処理
生データ取得 + 一次処理に必要なデータ取得
天体以外の信号の除去
感度と光学系のムラの補正

一次処理

一次処理済みデータ例) 点光源の測光 ( + 標準星データ取得 )

測光値の較正可視光
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方向への依存は無しと仮定し)

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

1.1

0.9

0.9

1.0

1.20.81.1

0.7

1.1200204201220183181

× 198 400 202 = 198 480 162205

199200226139220

可視光
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方向への依存は無しと仮定し)

ここがオーバースキャン領域

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

フラット
raw(x, y) = flat(x, y) x object(x, y) + dark(x, y) + bias(x, y)object(x, y) = constant (一様光)を観測してやれば、
ﬂat(x, y) =raw(x, y) - bias(x, y)
constant

constant = median(raw(x, y) - bias(x, y))
規格化

1次処理

ドームフラットデータから ﬂat(x, y)
天体生データから bias(x, y)object(x, y) =raw(x, y) - bias(x, y)flat(x, y)iraf.imcombine()
im01.ﬁts

im02.ﬁts

im03.ﬁts

101 104 102

105 101

99

9899

100 101

10099

102

103 1021039710198

96

105 101 103

100

102 10599

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

98100

iraf.imcombine()
im01.ﬁts

im02.ﬁts

im03.ﬁts

101 104 102

105 101

99

98

99

100 101

100

99

102

103 102

103

97

101

98

96

105 101 103

100

102 105
99

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

演習2
2-1
target2を観測した、'./data1/SUPA00317885.ﬁts' について、バイアス引き+フラット割りの処理をしましょう。これは5番フレームです。フィルターも同じBバンドなので、フラット割りには、'bﬂatn5.ﬁts'が使えます。
この結果のフレームを'btarget2n5.ﬁts'と呼ぶことにします。(後の演習で利用します)

2-2
2番フレームの生データ'./data2/SUPA00317702.ﬁts' について、バイアス引き+フラット割りの処理をしましょう。
先ほどの5番フレームとは違い、これは2番フレームなので、2番フレームのためのフラットを作成する必要があります。

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

講習3 : 星の測光
視野の星を検出させ、 iraf.daoﬁnd()
それらをアパーチャ測光
測光値の較正 ~ 標準星
演習3

iraf.phot()星の測光星を含む円内のカウントを総計する。
明るさの分かっている星のカウントと比べて、明るさを求める。

標準星の観測が必要

アパーチャ測光

FWHM
skyレベルをココで計算

アパーチャ半径

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

測光対象ファイル
iraf.daoﬁnd( ʻxxxxxx.ﬁtsʼ, output=ʻout1.cooʼ )検出結果のテキストファイル

測光対象ファイル

daoﬁndの結果のファイル

iraf.phot(ʻxxxx.ﬁtsʼ, coords=ʻout1.cooʼ
output=ʻout1.magʼ)
測光結果を書き出す(テキスト)ファイル

iraf.photの測光結果ファイル
iraf.txdump(ʻout1.magʼ,ﬁelds=ʻxc, yc, mag, merrʼ, Stdout=ʻresult1.txtʼ)抽出する項目

この結果の書き出し先

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

ʻbtarget2n5.ﬁtsʼの視野の中には測光標準星は写っていません。ただし、上のtarget1と近い時間に観測したデータですので、等級ゼロ点は同じだと仮定し、上と同じ較正値(1.672)を使ってください。
target1とは積分時間が異なりますので注意してください。

iraf.obsutil.psfmeasure講習4 : Numpyの基本
配列のベクトル演算 ~ ndarray
numpyの関数をいくつか二次元配列の注意点
演習4

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

1011 1220 21 2230 31 32

100 110 120+

200 210 220
300 310 320(2) 下の3つの二次元配列のメジアンを、numpy.stack()とnumpy.median()を用いて求めてください。

101112100 110 120202224202122200 210 220404244303132300 310 320606264

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

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

5-2. astropy.io.ﬁtsとnumpyを使って以下の処理をしてください。
(1) ʻ./data2/SUPA00317502.ﬁts'から2番フレーム用のフラットを作成。
(2) './data2/SUPA00317882.ﬁts' について、バイアス引き+フラッ
ト割りの処理(3) trimmingして(2)の結果からオーバスキャン部をとりのぞく。

講習6 : matplotlibの基本
ヒストグラム (光度関数)
X-Y プロット (等級 vs 等級エラー)
FITS画像の表示
FITS画像にオーバープロット
演習6

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