2009.10/6
一次処理：bias
DARK
flat-fielding
sky(background)引き
などがあるが、今回使ったSUBARU,MOIRCSの場合biasやDARKがよく分からないのでflatわりとsky引きだけを行う。更に今日は一枚に対してのみで行う。
☆其の一 flat割り(sky flat)
sky flatをとる際にskyを時間変化しない一様光として扱いたいが、どうしてもskyは時間変化(skyムラ)してしまう。一方今欲しいpixel毎のムラは時間変化してないとみなせるので、時間毎のskyのmedian値を一定値(1)に強引に合わせた上でpixel毎のムラをはかり、flat flameをつくる。
◎とりあえず必要なファイルのリスト作り
今天体も入ってる画像から天体を削ったsky flameを作るが、どうせ天体はカットするのでcl0024の画像にこだわる必要はない。なので一晩分のデータ(7/23に作ったcl0024_data_selectリストに情報入ってる)の中からexp timeがある程度長いものを全部使う。
$ awk '$4>100 && $4<1000{print $1".fits",$4}' cl0024_data_select | more
でまずちゃんとexp timeが100秒以上1000秒以下のものがはき出されたかどうか確認。
$ awk '$4>100 && $4<1000{print $1".fits",$4}' cl0024_data_select | cat -n | awk '$1%2==1{print $2}'| more
MCSA00038235.fits
MCSA00038237.fits
MCSA00038239.fits
MCSA00038241.fits
MCSA00038243.fits
MCSA00038245.fits
MCSA00038247.fits
・・・・・・・・・・
・・・・・・・・・・
・・・・・・・・・・
これで積分時間が長い奇数ファイルがはき出されたのでこれをobject1006oddとしてリスト化
$ awk '$4>100 && $4<1000{print $1".fits",$4}' cl0024_data_select | cat -n | awk '$1%2==1{print $2}' > object1006odd
同様に奇数リストも作る。
$ awk '$4>100 && $4<1000{print $1".fits",$4}' cl0024_data_select | cat -n | awk '$1%2==0{print $2}' > object1006even
◎sky時間変化の補正
・sky時間変動を可視化
cl> epar imstat
I R A F
Image Reduction and Analysis Facility
PACKAGE = imutil
TASK = imstatistics
images = @object1006odd List of input images //imstatにかける画像をリストとしてインプット//
(fields = midpt) Fields to be printed //median//
(lower = INDEF) Lower limit for pixel values
(upper = INDEF) Upper limit for pixel values
(nclip = 1) Number of clipping iterations //クリッピング回数 //
(lsigma = 3.) Lower side clipping factor in sigma //3-σ以下を除外、という意味//
(usigma = 3.) Upper side clipping factor in sigma //3-σ以上を除外。明るい天体はここでカットされる。//
(binwidt= 0.1) Bin width of histogram in sigma
(format = no) Format output and print column labels ?
(cache = no) Cache image in memory ?
(mode = ql)
:goで実行
時間毎のflameのmedian値(天体はカットされてるはず)が羅列される。
これをファイルに保存
cl> imstat @object1006odd nc=1 > med1006odd
偶数ファイルも同様に
cl> imstat @object1006even nc=1 > med1006even
これらをgnuplotでplotさせてやると→odd,even(epsファイル)
flame毎のskyの値が変動(時間変動)しているのがよくわかる。
・skyの高さ合わせ
今pixel毎のムラを知りたいので上記のsky時間変動を補正する。具体的には、sky値(flame毎のmedian値)を全部1にしてしまう。
とりあえず新しいリストのハコだけつくってから(中身はnormMCSA00038???.fits)
$ awk '{print "norm"$1}' object1006odd > norm1006odd
iraf上で元画像(object1006odd)をmedian値(med1006odd)で割ってやれば
cl>imarith @object1006odd / @med1006odd @norm1006odd
全flameのskyが１に規格化された画像リストnorm1006oddの完成
偶数ファイルも同様に
$ awk '{print "norm"$1}' object1006even > norm1006even
cl>imarith @object1006even / @med1006even @norm1006even
一応こちらもimstatでmedianとってやると→evenだけ(epsファイル)
◎flat flame作り→flat-fielding
これでnorm1006〜というリストの中の画像は全部skyの値は同じになってる。しかしまだ天体の画像は残ってる。ここから更に天体を除去してやればpixel毎のムラだけの画像=flat flameが出来るはず。
そこでimcombineというタスクを使って、norm1006〜内の全flameのmedian値をとった一枚の画像を作れば、天体は入り込まないより精度が良いflat flameができる！
（少し分かりにくい話だが、pixel毎にflame1とflame2とflame3,flame4,・・・のmedian値をとるということ）
cl> epar imcombine
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatch
TASK = imcombine
input = @norm1006odd List of images to combine //インプット//
output = flat1006odd.fits List of output images //アウトプットのfitsファイル名はここで決めとく//
(headers= ) List of header files (optional)
(bpmasks= ) List of bad pixel masks (optional)
(rejmask= ) List of rejection masks (optional)
(nrejmas= ) List of number rejected masks (optional)
(expmask= ) List of exposure masks (optional)
(sigmas = ) List of sigma images (optional)
(logfile= STDOUT) Log file
(combine= median) Type of combine operation //ここを今はmedianにしとく！//
(reject = sigclip) Type of rejection
(project= no) Project highest dimension of input images?
(outtype= real) Output image pixel datatype
(outlimi= ) Output limits (x1 x2 y1 y2 ...)
(offsets= none) Input image offsets
(masktyp= none) Mask type
more //この下にも続くが省略//
偶数チップも同様にflat1006even.fitsをつくる
cl> imcombine @norm1006even flat1006even.fits
結果をds9で見てみると
cl> display flat1006odd.fits 1 zs- zr- z1=0.8 z2=1.2→こんなの
cl> display flat1006even.fits 2 zs- zr- z1=0.8 z2=1.2→こんなの
けっこうpixel毎の感度ムラがあることが分かる。
・flat fielding
前回作ったcl0024の画像リストを今作ったflat flameで割ってやる
まずは出来上がりリストを用意(fl_field〜の中身はflMCSA00038???.fits)
$ awk '{print "fl"$1}' cl0024odd > fl_field_odd
$ awk '{print "fl"$1}' cl0024even > fl_field_even
それから
cl> imarith @cl0024odd / flat1006odd.fits @fl_field_odd
cl> imarith @cl0024even / flat1006even.fits @fl_field_even
これでflat fielding済みのcl0024の画像できた。
☆其の二 sky引き
とりあえず今日はflMCSA00038524.fits一枚に対してその前後５枚からsky flameを作ってsky引きしてみる。skyの時間変動の激しさを考えると、前後５枚ずつくらいがいいと思われる。
◎flMCSA00038524.fitsに対するsky flameの作成
（さっきと同じ手順でmedian値の高さを１にしたskyを作って、flMCSA00038524.fitsのmedian値倍してやれば良い）
$ vi 524
とうつとviエディタで524というファイル開くので挿入モードになりfl_field_evenからflMCSA00038524.fitsの前後５枚ずつをコピペしたものを貼り付ける。
flMCSA00038514.fits
flMCSA00038516.fits
flMCSA00038518.fits
flMCSA00038520.fits
flMCSA00038522.fits
flMCSA00038524.fits
flMCSA00038526.fits
flMCSA00038528.fits
flMCSA00038530.fits
flMCSA00038532.fits
flMCSA00038534.fits
~
更に挿入モードを抜け(esc)flMCSA00038524.fitsだけ消して(dd)保存(ZZ)
出来上がりはこうなったはず。
$ more 524
flMCSA00038514.fits
flMCSA00038516.fits
flMCSA00038518.fits
flMCSA00038520.fits
flMCSA00038522.fits
flMCSA00038526.fits
flMCSA00038528.fits
flMCSA00038530.fits
flMCSA00038532.fits
flMCSA00038534.fits
・次にこの１０枚それぞれのflameのmedian値を求める
cl> imstat @524 nc=1 > 524med
先と同様に各flameのmedian値を全部１に規格化してやる。
$ awk '{print "norm"$1}' 524 > norm524med //先に出来上がりのリストだけ作っとく//
cl> imarith @524 / @524med @norm524med
やはり先と同様にimcombineして高さが１のskyを作ってやる。
cl> epar imcombine
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatch
TASK = imcombine
input = @norm524med List of images to combine
output = norm524med.fits List of output images
(headers= ) List of header files (optional)
(bpmasks= ) List of bad pixel masks (optional)
(rejmask= ) List of rejection masks (optional)
(nrejmas= ) List of number rejected masks (optional)
(expmask= ) List of exposure masks (optional)
(sigmas = ) List of sigma images (optional)
(logfile= STDOUT) Log file
(combine= median) Type of combine operation
More
:go
一応作ったnorm524med.fitsの高さが１であることを確認
cl> imstat norm524med.fits
1.00035
・更に高さを今考えてるflMCSA00038524.fitsにそろえてやれば、flMCSA00038524.fits用のsky flameの完成！
cl> imstat flMCSA00038524.fits
36523.95 //flMCSA00038524.fitsのmedian値を調べる//
cl> imarith norm524med.fits * 36523.95 sky524.fits //高さ合わせ//
cl> imarith flMCSA00038524.fits - sky524.fits fl524sky.fits
これでflMCSA00038524.fitsからsky flame=sky524.fitsを引いた画像fl524sky.fitsの完成！→こんなの
