Mrk307撮像データの整約
練習用データは馬渡健のサブページ/データ倉庫にある。
必要なデータが入ってるディレクトリ：IRAF_im_practice
中はとりあえず
$ ls
FLAT_Im BIAS_Im M307_Im mrk307_imaging_BVRI.tar
☆事前準備
データのband毎に分類、リスト化
◎M307_Im内にて
・下調べ
cl> hselect FCSA* "frameid,object,filter01,exptime" yes > FCSA_data
FCSA00065498 Mrk307 SCFCFLBI01 120.0
FCSA00065500 Mrk307 SCFCFLBI01 120.0
FCSA00065502 Mrk307 SCFCFLBI01 120.0
FCSA00065504 Mrk307 SCFCFLBV01 120.0
FCSA00065506 Mrk307 SCFCFLBV01 120.0
FCSA00065508 Mrk307 SCFCFLBV01 120.0
FCSA00065692 Mrk307 SCFCFLBB01 120.0
FCSA00065694 Mrk307 SCFCFLBB01 120.0
FCSA00065696 Mrk307 SCFCFLBB01 120.0
FCSA00065698 Mrk307 SCFCFLBR01 120.0
FCSA00065700 Mrk307 SCFCFLBR01 120.0
FCSA00065702 Mrk307 SCFCFLBR01 120.0
→これより
FCSA00065498〜502：Iバンド
FCSA00065504〜508：Vバンド
FCSA00065692〜696：Bバンド
FCSA00065698〜702：Rバンド
・バンド毎にリスト化
$ awk '$3=="SCFCFLBI01"{print $1".fits"}' FCSA_data > FCSA_I_list
$ awk '$3=="SCFCFLBV01"{print $1".fits"}' FCSA_data > FCSA_V_list
$ awk '$3=="SCFCFLBB01"{print $1".fits"}' FCSA_data > FCSA_B_list
$ awk '$3=="SCFCFLBR01"{print $1".fits"}' FCSA_data > FCSA_R_list
◎FLAT_Im内のデータ
・下調べ
cl> hselect FCSA* "frameid,object,filter01,exptime" yes > lawFCSA_data
FCSA00064890 DOMEFLAT SCFCFLBB01 2.5
FCSA00064892 DOMEFLAT SCFCFLBB01 2.5
FCSA00064894 DOMEFLAT SCFCFLBB01 2.5
FCSA00064896 DOMEFLAT SCFCFLBB01 2.5
FCSA00064898 DOMEFLAT SCFCFLBV01 9.0
FCSA00064900 DOMEFLAT SCFCFLBV01 9.0
FCSA00064902 DOMEFLAT SCFCFLBV01 9.0
FCSA00064904 DOMEFLAT SCFCFLBR01 9.0
FCSA00064906 DOMEFLAT SCFCFLBR01 9.0
FCSA00064908 DOMEFLAT SCFCFLBR01 9.0
FCSA00064910 DOMEFLAT SCFCFLBI01 1.2
FCSA00064912 DOMEFLAT SCFCFLBI01 1.2
FCSA00064914 DOMEFLAT SCFCFLBI01 1.2
FCSA00064916 DOMEFLAT SCFCFLBI01 1.2
→FCSA00064890〜896：Bバンド
FCSA00064898〜902：Vバンド
FCSA00064904〜908：Rバンド
FCSA00064910〜916：Iバンド
・リスト化
$ ls FCSA* > lawFCSA_list
☆BIAS引き
$ cd BIAS_im
$ ls
FCSA00064940.fits FCSA00064944.fits FCSA00065962.fits FCSA00065966.fits
FCSA00064942.fits FCSA00065960.fits FCSA00065964.fits FCSA00065968.fits
◎準備
cl> wcsreset * wcs="world"
◎これから平均BIAS画像をつくる
・リスト化
$ ls > biaslist
・imcombie
cl> epar imcombine
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatch
TASK = imcombine
input = @biaslist List of images to combine
output = meanbias.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= average) Type of combine operation
(reject = none) 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
(:go)
◎天体画像からこの平均BIAS画像を引く
・meanbias.fitsを天体ディレクトリに移動
$ cp meanbias.fits /media/HD-CEU2/IRAF_im_practice/M307_Im/
（以下M307_Imディレクトリ内にて作業）
バイアス画像もリスト化（後でimarithするときのため）
$ awk '{print "meanbias.fits"}' lawdata_list > meanBIAS_list
・BIAS引き後のリスト作り
$ awk '{print "subBIAS"substr($1,10,20)}' lawdata_list > subBIAS_list
・imarithでBIAS引き
cl> imarith @lawdata_list - @meanBIAS_list @subBIAS_list
(※実は引く@lawdata_list - meanBIAS.fits @subBIAS_listでもイケる)
・バンド毎に分類
$ awk '{print "subBIAS"substr($1,10,20)}' FCSA_I_list > subBIAS_I_list
$ more subBIAS_I_list
subBIAS498.fits
subBIAS500.fits
subBIAS502.fits
てな感じ
他も同様に
$ awk '{print "subBIAS"substr($1,10,20)}' FCSA_V_list > subBIAS_V_list
$ awk '{print "subBIAS"substr($1,10,20)}' FCSA_B_list > subBIAS_B_list
$ awk '{print "subBIAS"substr($1,10,20)}' FCSA_R_list > subBIAS_R_list
☆FLAT fielding
ドームフラット画像がFLAT_Imディレクトリに入っているので、そこに移動
$ cd FLAT_Im
◎準備
cl> wcsreset * wcs="physical"
◎平均FLATフレーム(それぞれのフレーム内平均値は１)の作成
・まずはBIAS引き
$ awk '{print "subBIAS"$1}' lawFCSA_list > subBIASflat_list
cl> imarith @lawFCSA_list - meanBIAS.fits @subBIASflat_list
・このsubBIASflat_list内のデータをバンド毎に分ける
$ ls subBIAS*89[0,2,4,6].fits > subBIASflatB_list
以下どうにかしてsubBIASflatV_list,subBIASflatR_list,subBIASflatI_listをつくる
最終的には
$ more subBIASflatB_list
subBIASFCSA00064890.fits
subBIASFCSA00064892.fits
subBIASFCSA00064894.fits
subBIASFCSA00064896.fits
$ more subBIASflatV_list
subBIASFCSA00064898.fits
subBIASFCSA00064900.fits
subBIASFCSA00064902.fits
$ more subBIASflatR_list
subBIASFCSA00064904.fits
subBIASFCSA00064906.fits
subBIASFCSA00064908.fits
$ more subBIASflatI_list
subBIASFCSA00064910.fits
subBIASFCSA00064912.fits
subBIASFCSA00064914.fits
subBIASFCSA00064916.fits
となってればいい
（バンド毎にリストから平均FLATフレームをつくる）
・imcombine
cl> imcombine @subBIASflatB_list hogeFLAT_B.fits combine="median"
cl> imcombine @subBIASflatV_list hogeFLAT_V.fits combine="median"
cl> imcombine @subBIASflatR_list hogeFLAT_R.fits combine="median"
cl> imcombine @subBIASflatI_list hogeFLAT_I.fits combine="median"
（このままではFLATフレームのカウントが大きすぎるので、一枚の中の平均値が１になるようにしたい）
・平均値を求めて規格化
cl> imstat hoge* fields="image,mean"
# IMAGE MEAN
hogeFLAT_B.fits 7683.
hogeFLAT_I.fits 9217.
hogeFLAT_R.fits 10434.
hogeFLAT_V.fits 11456.
よって
cl> imarith hogeFLAT_B.fits / 7683 FLAT_B.fits
cl> imarith hogeFLAT_I.fits / 9217 FLAT_I.fits
cl> imarith hogeFLAT_R.fits / 10434 FLAT_R.fits
cl> imarith hogeFLAT_V.fits / 11456 FLAT_V.fits
（確認）
cl> imstat FLAT* fields="image,mean"
# IMAGE MEAN
FLAT_B.fits 0.9999
FLAT_I.fits 0.9999
FLAT_R.fits 1.
FLAT_V.fits 1.
◎FLAT fielding
M307_Imディレクトリに今作ったFLAT_?.fitsをコピー
$ cp FLAT* /media/HD-CEU2/IRAF_im_practice/M307_Im/
（以下M307_Imディレクトリで作業していく）
・種々のリスト化
$ awk '{print "fl"substr($1,8,20)}' subBIAS_B_list > fl_Blist
$ awk '{print "fl"substr($1,8,20)}' subBIAS_V_list > fl_Vlist
$ awk '{print "fl"substr($1,8,20)}' subBIAS_R_list > fl_Rlist
$ awk '{print "fl"substr($1,8,20)}' subBIAS_I_list > fl_Ilist
・subBIAS_?_listをFLAT_?listでわる。結果はfl_?listにつっこむ
cl> imarith @subBIAS_B_list / @FLAT_B.fits @fl_Blist
cl> imarith @subBIAS_V_list / @FLAT_V.fits @fl_Vlist
cl> imarith @subBIAS_R_list / @FLAT_R.fits @fl_Rlist
cl> imarith @subBIAS_I_list / @FLAT_I.fits @fl_Ilist
☆sky引き
◎安全な領域だけ切り取り
これから一枚の天体画像ごとにsky画像をつくり引いたり、天体検出をして位置合わせをしたりしていくが、その際に今の1024*2048の画像のままだと不都合が生じる。具体的には、sky作成の際にoverscan領域が邪魔だったり、[600,600]あたりの観測装置の反射を天体として検出してしまったりすることが考えられる。
以上のことより、Mrk307を含む都合の良い領域だけ切り出して、これ以降その画像を使うことにする。
都合の良い領域としては、例えば[100:990,650:1350]とか？
$ awk '{print $1"[100:990,650:1350]"}' fl_Blist > hogefl_Blist （imcopyインプット用）
$ awk '{print substr($1,1,5)"+.fits"}' fl_Blist > fl_B+list （imcopyアウトプット用）
ecl> imcopy @hogefl_Blist @fl_B+list
以下V,R,Iバンド画像についても同様に切り出し作業を行う
◎sky引き
今回はバンド毎に３枚しか画像がないので１枚１枚の画像からskyを二次元的に関数fittingし引く。という手段をとる。使うタスクはimsurfit
ecl> epar imsurfit
で見ればわかるがフィッティングする関数(func,xorder,yorder)や関数フィッティングからのrejection(lower,upper,niter)やフィッティングに使う領域(regions,)をいじれる。これを上手く使えば、天体をフィッティングから除去して関数フィッティングした画像＝skyが作れるはず。さらにtype_ouをresidualにすれば最初からskyを引いた後の画像が出力される。(fitは単に関数フィッティングしただけのsky画像)
ecl> epar imsurfit
I R A F
Image Reduction and Analysis Facility
PACKAGE = imfit
TASK = imsurfit
input = fl692+.fits Input images to be fit
output = subsky692.fits Output images
xorder = 2 Order of function in x
yorder = 2 Order of function in y #skyは大体一定になってるはずなので、一次式でフィッティング#
(type_ou= residual) Type of output
(fit,residual,response,clean) #clean(フィッティング上ならした外れピクセルを見せてくれる)などにもして確認した方が良い#
(functio= spline3) Function to be fit(legendre,chebyshev,spline3) #関数形を選べる#
(cross_t= yes) Include cross-terms for polynomials?
(xmedian= 1) X length of median box
(ymedian= 1) Y length of median box
(median_= 50.) Minimum fraction of pixels in median box
(lower = 2.) Lower limit for residuals
(upper = 2.) Upper limit for residuals #フィッティングの際にrejectionするlowσとhighσ#
(ngrow = 0) Radius of region growing circle
(niter = 5) Maximum number of rejection cycles #σclippingする回数#
(regions= border) Good regions (all,rows,columns,border,sections,c #フィッティングに使う領域。今回は中心銀河を避けたいので「へりだけ使う」border#
(rows = *) Rows to be fit
(columns= *) Columns to be fit
(border = 250) Width of border to be fit #何ピクセル分へりを使うか指定#
(section= ) File name for sections list
(circle = ) Circle specifications
(div_min= INDEF) Division minimum for response output
(mode = ql)
(:go)
→結果画像を見てimexamなどで天体の写ってない領域が大体0になっていれば成功！
他の画像も同様にやれば良いが、上手くいくかどうかはケースバイケースなのでリスト化して一気にやるより一枚一枚確認しながらの方が良いだろう。
☆位置合わせ
SExtractorを用いた天体検出→同一天体の対応づけ(xyxymatch)→位置合わせの変換規則作り(geomap)→画像変換(geotran)という手順で行う。
（以下とりあえずBバンドデータのみで行う）
◎天体検出
・パラメータファイル作り
SExtractor入っている場合は適当なパラメータファイルを以下のように作る。（最初は~/sextractor2.5.0/config/からdefault.sexやdefault.paramをコピーしてきて書き換えるのを推奨する）
→detection.sex,detection.param
DETECT_THRESH=6,DETECT_MINAREA=10は明らかにしばりがきついが,主鏡の反射を天体として検出したくない＋xyxymatchでこけないようにわざとこうしてる
・SExtractorの実行
$ sex subsky692.fits -c detection_test.sex
$ mv test.cat detec692.cat （結果ファイルの名前変更）
$ more detec692.cat
→検出天体数が多すぎたり少な過ぎたりしないか確認
ecl> display subsky692.fits 1
ecl> display check.fits 2
→Tab連打で天体検出され具合を確認
（同様に他の二枚694,696に対しても天体検出をかける）
$ sex subsky694.fits -c detection.sex
$ mv test.cat detec694.cat
$ sex subsky696.fits -c detection.sex
$ mv test.cat detec696.cat
◎同一天体の対応表作り(xyxymatch)
・事前準備
今回は天体の検出数が30個くらいしかないので全部使うが、上部の説明書きを削ったファイルを作っとく。
$ awk '$1>0 && $1<100{print}' detec692.cat > detec692+.cat
$ awk '$1>0 && $1<100{print}' detec694.cat > detec694+.cat
$ awk '$1>0 && $1<100{print}' detec696.cat > detec696+.cat
・xyxymatch
以下reference=692画像に694や696を合わせていく
ecl> epar xyxymatch
I R A F
Image Reduction and Analysis FacilityM
PACKAGE = immatch
TASK = xyxymatch
input = detec694+.cat The input lists
referenc= detec692+.cat The reference lists
output = match694_692.cat The output matched coordinate lists
toleranc= 3. The matching tolerance in pixels
(refpoin= ) Optional list of reference points
(xin = INDEF) X origin of input list
(yin = INDEF) Y origin of input list
(xmag = INDEF) X magnification required to match input to refer
(ymag = INDEF) Y magnification required to match input to refer
(xrotati= INDEF) X rotation required to match input to reference
(yrotati= INDEF) Y rotation required to match input to reference
(xref = INDEF) X origin of reference list
(yref = INDEF) Y origin of reference list
(xcolumn= 5) Input list column containing the x coordinate
(ycolumn= 6) Input list column containing the y coordinate
(xrcolum= 5) Reference list column containing the x coordinat
(yrcolum= 6) Reference list column containing the y coordinat ＃インプットとリファレンスファイルのx,y座標に対応する列を指定＃
(separat= 9.) The minimum object separation
(matchin= triangles) The matching algorithm
(nmatch = 30) The maximum number of points for triangles algor
(ratio = 10.) The maximum ratio of longest to shortest side of
(nreject= 10) The maximum number of rejection iterations #triangleアルゴリズム内でのrejection回数。detec696+.catにxyxymatchかけるときはここを50にしたら上手くいった#
(xformat= %13.3f) The format of the output x coordinate
(yformat= %13.3f) The format of the output y coordinate
(interac= no) Interactive mode ?
(verbose= yes) Verbose mode ?
(icomman= ) The image display cursor
(mode = ql)
(:go)
→694と692内の同一天体の対応表match694_692.catができたはず
$ more match694_692.cat
でみて、今回ならばx方向のズレしかなければ成功。
同様に
ecl> xyxymatch detec696+.cat detec692+.cat match696_692.cat 3 nreject=50
などとして696と692の対応表も作っとく。
◎位置合わせの変換規則作り(geomap)
(xyxymatchさえ上手くいっていれば問題なく進めるはず)
ecl> geomap match694_692.cat db694_692.cat 1 891 1 701 xxorder=2 xyorder=2 yxorder=2 yyorder=2
("1 891 1 701"はreferenceのフォーマットを示している)
→subsky694.fitsをsubsky692.fitsに位置合わせするための画像変換規則の情報入ったdb694_692.catというファイル出来たはず
同様に
ecl> geomap match696_692.cat db696_692.cat 1 891 1 701 xxorder=2 xyorder=2 yxorder=2 yyorder=2
◎画像幾何変換（geotran）
「subsky694.fits→692と天体が一致するような画像」、「subsky696.fits→692と天体が一致するような画像」の具体的な変換を行う
ecl> geotran subsky694.fits shift694.fits db694_692.cat match694_692.cat
ecl> geotran subsky696.fits shift696.fits db696_692.cat match696_692.cat
・確認
ecl> display subsky692.fits 1 zr- zs- z1=-10 z2=5000
ecl> display shift694.fits 2 zr- zs- z1=-10 z2=5000
ecl> display shift696.fits 3 zr- zs- z1=-10 z2=5000
match Flameで画像あわせといてから、Tab連打するとしっかり位置合わせされたか確認出来る。
◎３枚の画像の足しあわせ→バンド毎にS/N比の良い画像を作る
・位置合わせした３枚の画像(subsky692.fits,shift694.fits,shift696.fits)のリスト作成
$ ls shift*.fits > Bshifted_list
$ vi Bshifted_list で１行目にsubsky692.fitsを加える
どんな方法でも良いので、中身が
$ more Bshifted_list
subsky692.fits
shift694.fits
shift696.fits
となっていればOK
・足しあわせ
ecl> imcombine @Bshifted_list B.fits combine=median reject=sigclip lsigma=3 hsigma=3
結果見てみると
ecl> display B.fits 4 zr- zs- z1=-10 z2=5000 →よさげ
以上の天体検出から位置合わせした画像の足しあわせまでの一連の作業を他のバンド(V,R,I)画像にも行ってやれば、精度の良い画像B.fits,V.fits,R.fits,I.fitsが出来るはず
☆おまけ＝＞三色合成
せっかくB(blue),V(green),R(red)バンドがそろっているのだから色付けして綺麗な画像にしてみる。
※ただしこの作業から科学的な情報は得られない。科学的にはバンド毎の等級などが分かればカラー求まるので・・・
◎colorパッケージのインストール
三色合成を行うためのcolorというパッケージはirafにデフォルトで入っているものではないので、以下のページを見て下さい
パッケージのインストールメモ
◎三色合成
・colorパッケージに入る
ecl> color
・rgbsunというタスクを使う
color> epar rgbsun
I R A F
Image Reduction and Analysis Facility
PACKAGE = color
TASK = rgbsun
red = R.fits Red image #赤色に対応するバンド画像#
green = V.fits Green image #緑色#
blue = B.fits Blue image #青色#
rgb = Mrk307.ras Output RGB image #出来上がり。拡張子はrasにしとく#
(rz1 = 0.) Red z1 #red画像の明るさの下限値#
(rz2 = 100.) Red z2 #red画像の明るさの上限値#
(gz1 = 0.) Green z1
(gz2 = 100.) Green z2 #green画像の明るさの上限下限#
(bz1 = 0.) Blue z1
(bz2 = 100.) Blue z2 #blue画像の明るさの上限下限#
(logmap= no) Use logarithmic intensity mapping?
(swap = yes) Swap bytes in output rasterfiles? #必ずyesにしとく#
(mode = ql)
(:go)
※z1,z2について・・・各バンド画像の出力最小最大値をを指定。 z2 の値を大きくとると暗い画像に、小さくとると明るい画像になる。 z1 と z2 の値の開きを広くすると軟らかい画像に、開きを狭くすると硬い画像になる。更に色調が変わる。例えばbz2だけ小さくすれば、青色が強い画像になったりする。
→Mrk307.rasという色のついた画像ができる
◎.ras→.jpgの変換
rbgsunで作った画像はSun 24-bit raster fileというよく分からんファイル形式(拡張子ras)で出てくる。これをよく使うjpgやgimpにするには以下のようにgimp(GNU Image Mainpulation Program)から
変換するとラク
$ gimp &
現れたgimpの"ファイル/開く"からMrk307.rasを開き、"ファイル/名前をつけて保存"でJpegなどを選び,Mrk307.jpgなどと出力ファイル名を指定すればjpg画像になる。
→例えばこんな感じ