2013年1月にMagellan/MEGACAMデータの処理を行った際のメモ。
バイアス/ダーク引き、フラット割り、宇宙線処理、illumination補正、スカイ引き、(フリンジ処理)などまではパイプラインにより行われている事を想定。
大まかには
・Multi Extension fits => Single Extension fits への変換
・重ね合わせ
の２処理のみだが、重ね合わせのスケーリング則を求めたりしなければならない。
使用pcは、MacPro(mountain lion OS X 10.8.2)@CfA/OIR と CF linuxマシン(mars)
!!中途半端にCfAの内部コンピュータの環境を利用しているため、今の所自分用のメモ。!!
☆参考：
Matt AshbyによるMEGACAM data reduction cookbook
CfA MEGACAM Pipeline Wiki
０、Multi Extension fits の扱い方
自身最初は良く分からなかったので
◎ds9での見方
・ds9表示させた状態から
Fileタブ
↓
Open Other
↓
Open Mosaic IRAF
・コマンド端末から
% ds9 -mosaic iraf ファイル名 &
◎iraf上での操作
・普通の状態
ファイル名の最後に何番目のextensionかを付ける
(例)
vocl> imhe MEXhoge.fits[3]
これはMEXhoge.fitsというMulti Extension fitsの３つ目を扱うという意味。
・mscredパッケージのタスクを使う場合
何番目のextensionかは気にせず普通に操作
(例)
mscred> mscarith MEXhoge.fits * 100 MEXaho.fits
1.funpack
データはfpackで圧縮されてくるので、funpackで解凍してやる。
"MEX????.fits.fz"というような名前の圧縮Multi Extension fitsファイルが複数毎あ ると考えて以下、進める。
% funpack MEX????.fits
でfzファイルが解凍されてfitsになる。
fzファイルを残さなくて良い場合は
% funpack -D MEX????.fits
など。
詳しくは
% funpack
% funpack -H
と打てば色々出てくる。
2.Multi Extension fits => Single fits 変換
◎WCS投影方法をZpnに直す
たぶん何もしなくてもそうなっていると思うが。headerに書かれている
直す場合は以下のようにする。
・clスクリプト化
% ls MEX*.fits > MEX.lis
% sed 's/MEX/zpn MEX/g' MEX.lis | sed 's/fits/fits; flpr/g' > run_zpn.macro
run_zpn.macroの中身は以下のようなもの
zpn MEX???1.fits; flpr
zpn MEX???2.fits; flpr
zpn MEX???3.fits; flpr
zpn MEX???4.fits; flpr
・・・・・
・・・・・・・・・・
・・・・・・・・・・・・・・・
・実行
megacam> cl< run_zpn.macro
◎SWarpによりSingle fitsに変換
・configurationファイルの作成
% swarp -d > swarp_one.def
基本的にはこのデフォルトのままでOKだが、以下の部分は変える
% emacs swarp_one.def
COMBINE_TYPE AVERAGE #ほとんどのケースで重要ではない
PROJECTION_TYPE TAN #絶対
CENTER_TYPE MANUAL #絶対
CENTER hh:mm:ss.s, +dd:mm:ss.s #自分で画像のセンターのRADec を決める事
PIXELSCALE_TYPE MANUAL #絶対ではないと思う
PIXEL_SCALE 0.16 #PIXELSCALE_TYPEをMANUALにした場合はここの値が適用されるので、MEGACAMのピクセルスケールを入れる
IMAGE_SIZE 13000 #画像サイズの指定。あまり大きいと後で困るかも
FSCALASTRO_TYPE FIXED #絶対
SUBTRACT_BACK Y #スカイ引きをこの段階で行うかどうかは状況に応じて
BACK_SIZE 128 #Matt推奨値
BACK_FILTERSIZE 3 #Matt推奨値
BACK_FILTTHRESH 5 #Matt推奨値
・shスクリプト化
% swarp MEX????.fits -c swarp_one.def -IMAGEOUT_NAME swarp????.fits-WEIGHTOUT_NAME swarp????.weight.fits
を複数毎分、一気に行うようなスクリプト組む
例えば、
% ls MEX*.fits > MEX.lis
% awk '{print "swarp",$1,"-c swarp_one.def -IMAGEOUT_NAME","swarp"substr($1,4,20),"-WEIGHTOUT_NAME","swarp"substr($1,4,4)".weight.fits"}' MEX.lis > swarp_all.sh
など
・実行
% sh swarp_all.sh
◎いらない部分の切り取り
SWarpは正方形の画像しか作れない事もあり、実際には余白があったりする。MEGACAMのデータはとにかくサイズが大きく、少しでもサイズを小さくして次に進んだ方が時間の節約にもなるし、余計なエラーを出さずに済む。
・必要な部分を調べる 色々方法はあるが、自分は一回適当に重ね合わせて、出来上がりの１枚を見て決める事にした。
% ls swarp*.fits > swarp.lis
vocl> imcomb @swarp.lis hoge_coadd.fits combine=average
% ds9 hoge_coadd.fits
↓
必要な[x1:x2,y1:y2]を決める
・imcopyのclスクリプト化
% awk '{printf "!mv %s tmp_%s\nimcopy tmp_%s[x1:x2,y1:y2] %s\n!rm tmp_%s\n",$1,$1,$1,$1,$1}' swarp.lis > swarp_imcopy.macro
(weight mapも)
% ls swarp*.weight.fits > weight.lis
% awk '{printf "!mv %s tmp_%s\nimcopy tmp_%s[x1:x2,y1:y2] %s\n!rm tmp_%s\n",$1,$1,$1,$1,$1}' weight.lis > weight_imcopy.macro
・実行
vocl> cl< swarp_imcopy.macro
vocl> cl< weight_imcopy.macro
3.Exposure毎の相対的なフラックス比(relative zeropoint)の推定
◎原理
airmass、天候(seeing)、積分時間など様々なファクターにより、各フレームにおける１カウントが持つ物理情報(フラックス)は異なる。同じ星を見ていても天気が悪い時のフレームの方がカウントが小さかったりするわけである。こうした情報はzero点と呼ばれる以下の値に全て含まれ、大気の透過率が悪かったり積分時間が短いとzero点は小さくなる。
mag = -2.5 × log(Count) + ZERO
重ね合わせる時に、フレーム毎のzero点の違いを考慮しないと物理的に意味が無くなってしまう。 原理的には、異なるフレームで同一の星を測光し、カウント値の違いから補正量<=>相対的なzero点を求める。
例えばフレーム１と対等なフレームになるようにフレーム2の補正量を求める事を考える。ある等級の星をそれぞれのフレームで見た場合、
mag = -2.5 × log(C1) + Z1 = -2.5 × log(C2) + Z2
補正後のフレーム2(2')ではzero点がフレーム1と同じになってほしいので、
mag = -2.5 × log(C1) +Z1 = -2.5 × log[C2×10^(-0.4ΔZ2)] + Z1
ここでΔZ2は相対的なzero点であり、ΔZ2=(Z2-Z1)である。
この式が意味する所は、フレーム2に10^(-0.4×ΔZ2)をかけると、フレーム1と等価な「カウント<=>フラックス」関係のフレームになるという事である。
各フレームのzero点の絶対値(Z1とZ2)が分からなくても、同じ星を見た時にC1=C2×10^(-0.4ΔZ2)になるという要求から相対的なzero点を求める事ができる。
Megared内にはこれを全自動で行ってくれる"zptmatch"というタスクがあるので、以下それを使う。
◎zptmatch(CFマシン内での作業 )
(MegaredのzptmatchタスクがLocalに動かないので・・・・)
・データの転送準備
CFマシンの割当容量は15GBまでと限られているので、小分けにして、更に転送時間が短くなるようにfpack
% fpack swarp*.fits
(小分け用ディレクトリ作って分配)
% mkdir mosaicking1
% mv swarp000?.fits mosaicking1
% mkdir mosaicking2
% mv swarp001?.fits mosaicking2
% mkdir ・・・・・
% mv ・・・・・・・
(リファレンスフレーム：番号が一番若いもの、例えばswarp0001.fitsを各小分けディレクトリに分配)
% cp mosaicking1/swarp0001.fits mosaicking2/
% cp mosaicking1/swarp0001.fits mosaicking3/
% cp mosaicking1/swarp0001.fits ・・・・・・・
% cp ・・・・・・・・・・・・・・・・・・・・・・・
!以下小分けにしたディレクトリ毎に!
・転送
% scp -r mosaicking1/ ユーザー名@mars.cfa.harvard.edu:~/obs/
・CFマシン内でfunpack
% ssh -X ユーザー名@mars.cfa.harvard.edu
mars> cd obs/
mars> cd mosaicking1
funpackは一枚ずつしか解凍できないので適当にスクリプト化する事
=> funpack.sh
mars> ls *fz > fz.lis
mars> sed 's/fits.fz/fits/g' fz.lis > fits.lis
mars> sh funpack.sh fits.lis
・zptmatchの実行
基本は
mars> zptmatch swarp????.fits
とやるだけで勝手にZpt.outというrelative zero点ファイルができる。
しかし全自動で本当に正しく動作しているか不安なのでチェックその他も合わせて以下のようにスクリプト化(時間は倍になるが・・・)
=> zptmatch.sh
中身は、毎回変えなければならず使い勝手は悪いが・・・
zptmatch swarp????.fits
#↑対象のファイル全てを含むように「swarp????.fits」の部分を工夫して書く必要有#
ls swarp????.fits > hoge.lis
#同上#
paste hoge.lis Zpt.out | awk '{print $1,$3}' > swarp_1.zeropoints
#ここは出来上がりの「swarp_1.zeropoints」を小分けにしたディレクトリに合わせて毎回違う名前をつける(今はmosaicking1)#
cp Zpt.dir/zptfit.0.cat ./zptfit_1.0.cat
#使った星のデータファイル、やはり毎回違う名前で保存#
zptmatch -m 10 swarp????.fits
#チェック用にパラメータをいじる。また敢えてリファレンスファイルを含めないように「swarp????.fits」の記述で指定#
ls swarp????.fits > aho.lis
paste aho.lis Zpt.out | awk '{print $1,$3}' > check_1.zeropoints
#チェック用の保存、毎回違う名前#
cp Zpt.dir/zptfit.0.cat ./check_1.0.cat
#毎回違う名前#
mars> sh zptmatch.sh
これでできたファイルのうち一番重要なのは swarp_1.zeropointsであり
「ファイル名 relative-zero点」
が中に書かれているはず
・出来たファイルをlocalマシンに転送
(localマシンから)
% scp ユーザー名@mars.cfa.harvard.edu:~/obs/swarp_1.zeropoints ./
% scp ユーザー名@mars.cfa.harvard.edu:~/obs/zptfit_1.0.cat ./
% scp ユーザー名@mars.cfa.harvard.edu:~/obs/ ./check_1.zeropoints
% scp ユーザー名@mars.cfa.harvard.edu:~/obs/ ./check_1.0.cat
!ここまでが１小分けディレクトリ分!
・小分けディレクトリ分だけ以上の作業を繰り返す
※(注意)常にreferenceとするファイルとして同じ物を使いたいので、一番目のファイル(swarp0001.fitsなど)は全小分けディレクトリにある事を確認してから行う事
◎localマシン内で一つのzero点ファイルにマージ
% cat swarp*.zeropoints > swarp.zeropoints
% emacs swarp.zeropoints
適当な絵エディタで開いて、不必要な行、つまり二つ目以降の「swarpECDFS1019_2535.fits 0」を消す。
中身は
swarp0001.fits 0
swarp0002.fits 0.0272265
swarp0003.fits 0.0246331
・・・・・・・・・・・・・・
・・・・・・・
(check用zero点ファイルも)
% cat check*.zeropoints > check.zeropoints
◎相対zero点の妥当性チェック
・本番zero点とチェック用zero点の整合性
本番zero点(swarp.zeropoints)はzptmatchの設定がデフォルトで、リファレンスはswarp0001.fitsに統一されたもの。check用(check.zeropoints)は(-m 10)でリファレンスを小分けディレクトリ毎に変えてある。
% cat -n swarp.zeropoints > Zpt_swarp.cat
% cat -n check.zeropoints > Zpt_check.cat
% cat -n check.zeropoints | awk '{if ($3==0) printf "\n\n\n%d %s %lf\n",$1,$2,$3; else print}' > Zpt_check.cat
gnuplot等で図示、比較
% gnuplot
gnuplot> set term x11
gnuplot> plot "Zpt_swarp.cat" u 1:3 pt 7 lc 1, "Zpt_swarp.cat" u 1:3 w l lc 1, "Zpt_check.cat" u 1:3 pt 7 lc 3, "Zpt_check.cat" u 1:3 w l lc 3
下図のように、赤線(本番)と青線(check用)で、絶対値の上下はあっても、小分けディレクトリ毎の変動パターンが大体一緒なら問題無し
・使用した星をチェック
zptfit????.0.catに使用した星の情報が入っている(9行目：x、10行目：y、12行目：フレーム番号)ので、画像上でどのような星が相対zero点推定に使われたかチェック。サチった星を使っているようだと問題
・時間変動
元のMulti Extension fitsのヘッダーから観測日時、SEEING、AIRMASS情報を入手
vocl> hselect MEX????.fits[1] fields=FILENAME,DATE-OBS,AIRMASS,SEEING expr+ > UT-AIRMASS-SEEING.data
中は以下のようになっている
MEGACAM/2012.1115/ECDFS2.3951 2012-11-15T07:16:33 1.221 0.692487
MEGACAM/2012.1115/ECDFS2.3952 2012-11-15T07:22:04 1.240 0.700011
MEGACAM/2012.1115/ECDFS2.3953 2012-11-15T07:27:36 1.259 0.688745
これを扱いやすくするために以下のように書き換える
% awk '{print substr($1,25,5),substr($2,1,4),substr($2,6,2),substr($2,9,2),substr($2,12,2),substr($2,15,2),substr($2,18,2),$3,$4}' ECDFS_UT-AIRMASS-SEEING.data > ECDFS_UT-AIRMASS-SEEING.cat
中は以下のようになっている
3951 2012 11 15 07 16 33 1.221 0.692487
3952 2012 11 15 07 22 04 1.240 0.700011
3953 2012 11 15 07 27 36 1.259 0.688745
UTをCLST(チリ時間)に変更し、更にswarp.zeropiontsの情報を加える
% paste ../swarp.zeropoints ECDFS_UT-AIRMASS-SEEING.cat| awk '{print $1,$4,$5,$6,$7-4,$8,$9,$10,$11,$2}' > ECDFS_CLST+ASZ.data
中身が「ファイル名、年、月、日、時、分、秒、AIRMASS、SEEING、相対zero点」になるように
zero点の時間変動を図示
% gnuplot
gnuplot> set term x11
gnuplot> plot "ECDFS_CLST+ASZ.data" u (($4==19)?($5+(($6+$7/60)/60)):1/0):10 w l lt 1 lc 0 t ""
など。例えば下図のような感じ。
・SEEINGやAIRMASSなどの観測条件の時間変動と比較
SEEINGやAIRMASSの時間変動のグラフも作ってみて、相対ゼロ点の時間変動が妥当であるか考察。
4.マスクがけ
3まででSingle fitsのフレームが出来ており、exposure無しの部分は0になっているのでこのまま重ね合わせてもいいのだが、画像サイズが大きすぎてirafがcrashしてしまう事があるらしい。原因不明だが、exposure無しのピクセルの部分がトリガーになっているらしい。
最近のirafではcrashしないかもしれないが、一応exposure無しのピクセルを-10000.0にする。また衛星・飛行機の通過跡などにも同様にマスクがけを行う。
◎exposure無領域のマスクがけ
SWarpで作成したweghtファイル(swarp????.weight.fits)があればimfillで簡単にできる。weightファイルの0の部分を-10000.0にしてしまえばいいだけ
・clスクリプト化
% ls swarp????.fits > swarplist
% awk '{printf "imarith %s.weight.fits * 10000 tmp_%s.weight.fits\nimfill %s tmp_%s.weight.fits \"mask .eq. 0\" -10000.0\n!rm tmp_%s.weight.fits\n",substr($1,1,9),substr($1,1,9),$1,substr($1,1,9),substr($1,1,9)}' swarplist > imfill.macro
中身は
imarith swarp0001.weight.fits * 10000 tmp_swarp0001.weight.fits
imfill swarp0001.fits tmp_swarp0001.weight.fits "mask .eq. 0" -10000.0
!rm tmp_swarp0001.weight.fits
imarith swarp0002.weight.fits * 10000 tmp_swarp0002.weight.fits
imfill swarp0002.fits tmp_swarp0002.weight.fits "mask .eq. 0" -10000.0
!rm tmp_swarp0002.weight.fits
imarith swarp0003.weight.fits * 10000 tmp_swarp0003.weight.fits
imfill swarp0003.fits tmp_swarp0003.weight.fits "mask .eq. 0" -10000.0
!rm tmp_swarp0003.weight.fits
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・
・実行
vocl> stsdas
stsdas> cl< imfill.macro
※imfillは入力ファイルを上書きするので慎重に
◎衛星等の通過跡(直線)のマスクがけ
・全フレームのvisual check
どのフレームのどこに直線状のマスクした方がいいfeatureがあるかをds9で画像見ながらチェック
面倒くさいが、CFマシンでzptmatchしてる間など手持ち無沙汰な時間にやると良い。
・satbgonのclスクリプト化
以下のようなclスクリプト(satbgon.macro)を組めばいい
satbgon swarp0004.fits x1 y1 x2 y2 width -10000.0
#↑x1,x2,y1,y2,widthにvisualチェックで調べた値を入れる#
satbgon swarp0010.fits x1 y1 x2 y2 width -10000.0
satbgon swarp0024.fits x1 y1 x2 y2 width -10000.0
satbgon swarp0031.fits x1 y1 x2 y2 width -10000.0
・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・
・実行
megacam> cl< satbgon.macro
※これもやはり入力ファイルを上書きするので慎重に。作業スペースに余裕があれば、上書きしないようにsatbgon.clを修正しても良いかも
5.重ね合わせ
imcombineするわけだが、その際の各フレームの標準偏差による重み付けやBad Pixel Maskなどの情報を付加する。
◎fitsヘッダーにBPM情報としてexposure無し領域(weightマップ)を登録
・作成
% ls swarp????.fits > swarplist
% ls swarp????.weight.fits > weightlist
% paste swarplist weightlist| awk '{printf "hedit %s BPM %s add+ delete- update+ verify-\n",$1,$2}' > hedit.macro
% awk '{printf "hedit %s BLANK -10000.0 add+ delete- update+ verify-\n",$1}' swarplist >> hedit.macro
・実行
cl< hedit.macro
ここまでやったら、もうSWarpで作ったweight.fitsは消してしまってもいい。
◎SEEING情報の調査
・header情報
既に調べてある(ECDFS_CLST+ASZ.dataなど)が、時々"-1"などいかれた値が書いてあって使いづらい。
・自前で測り直し
「SExtractorで天体のFWHM出力、flagやclass star,領域しばりで星っぽいものを抽出、FWHMの統計を調べる」プログラム作成
↓
est_seeing.cl
実行
% ls swarp????.fits > swarp.lis
cl> cl< est_seeing.cl
アウトプットのファイル：seeing.dataの中身は
# No name used_str_num mean midpt mode stddev min max
1 ../swarpECDFS1019_2535.fits 100 0.93986 0.94014 0.94639 0.03194 0.75840 1.07040
2 ../swarpECDFS1019_2536.fits 24 1.06647 1.07530 1.07565 0.04186 0.91520 1.15840
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・
・チェック
自分で測ったseeing(ECDFS_seeing.data)とheaderのseeingがconsistentかどうか図示して確認
% awk '$1!="#"{print}' ECDFS_seeing.data > ECDFS_seeing_NC.data
% paste ECDFS_CLST+ASZ.data ECDFS_seeing_NC.data| awk '{print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$14,$15,$16}' > ECDFS_CLST+ASZmyS.data
中身は「Name,yyyy,mm,dd,CLST/h,CLST/m,CLST/s,AIRMASS,SEEING,relaZPT,seeing(mean),seeing(midpt),seeing(mode)」
図示すると以下のようになる.
大体傾向としてはheader値とconsistent
なおケースバイケースだが、ここでは自分で測ったseeingとしてmedian値を採用している。
◎imcombineのスケーリング則、重み付け情報
・インプットリスト
% ls swarp????.fits > swarp.inlist
・スケーリング則
相対ゼロ点から
% awk '{printf "%9.7f\n",(10**(-0.4 * $2))}' swarp.zeropoints > swarp.scales
・重み付け情報
各フレームの標準偏差(σ)とシーイング情報(FWHM)を用いる => (1/scale/σ/FWHM)^2で重み付けする
vocl> imstat @swarp.inlist fields="image,npix,mean,midpt,mode,stddev,min,max" format- lower=-70.0 upper = 70.0 > swarp.stats
% paste swarp.scales swarp.stats seeing_NC.data | awk '{printf "%6.4f\n", (100.0/$1/$1/$7/$7/$14/$14)}' > swarp.weights
◎imcombine
パラメータ設定は以下のように
imcombine.hsigma=3.0
imcombine.lsigma=3.0
imcombine.offsets = "none"
imcombine.masktype = "none"
imcombine.lthresh = -300.0
imcombine.zero = ""
imcombine.weight = "@swarp.weights"
imcombine.scale = "@swarp.scales"
imcombine.nrejmasks="fin_coadd.nrej"
imcombine.expmasks="fin_coadd.cov"
imcombine.outlimits=""
imcombine.combine= "average"
imcombine.reject ="avsigclip"
imcombine.expname = "EXPTIME"
・実行
megacam> imcombine @swarp.inlist output="fin_coadd.fits" bpmask="fin_coadd.bpm"
megacam> imcopy fin_coadd.cov.pl fin_coadd.cov.fits
6.PSF size 測定
最終画像のPSFチェック
◎大体の目安をつける
・iraf/imexamやSExtractorなど
cl> imexam
"a"でPSFなどの情報出てくるが、自分の扱っている画像の場合、大体、FWHM~0.6''
各フレームを重ね合わせる時にSeeingで重み付けしているため、Seeingがいい所の重みが大きく、この値は妥当
◎プログラム
・先に作ったest_seeing.clなどを使う
※注意・・・SExtractorベースのPSFになるので、irafで測ったFWHMと多少はズレるかもしれない
% ls ECDFS_coadd_r_v2.fits > coadd_v2.lis
適当なSExtractorファイル(~.sex,~.param,~.nnw,~.conv)を同じディレクトリに置く。重要なパラメータはest_seeing.clの中で決められるので気にしなくてもいい。
est_seeing.clの中で大体のseeingを入れる欄がある(seeing_sex)のでそこにirafで測った大体の値を入れる。星か銀河かの区別をするためのパラメータになっており、出力の星数はこのパラメータに敏感なので、少し気を配った方がいい。
cl> cl< est_seeing.cl
・確認
出力はPSFデータの入ったファイルと使用した星の情報の入ったファイルの２種類
PSFファイルの中身は
% less ECDFSv2_PSF.data
# No name used_str_num mean midpt mode stddev min max
1 ECDFS_coadd_r_v2.fits 200 0.63385 0.61184 0.63480 0.08736 0.55360 1.03200
となっている。PSFはmean,midpt,mode(最頻値)の３種類が表示されているが、irafでの値と大きくはズレていない事を確認
次に星情報ファイル(str_1.cat)はSExtractorパラメータファイルで指定された情報が格納されており、「id,aperture-mag,kron-mag,kron-radius,background,ra,dec,ellipticity,FWHM(pix),star_class」となっている。
regionファイル化
% awk '{print $6,$7}' str_1.cat > str_1.reg
ds9上で確かに星を拾っているかなどを確認。またそれらに対してirafでimxamしてプログラムの出力値と同じ程度かなどを確認。(大体のケースでirafで測った方が少し小さめになる・・・)
使用した星のFWHMを全部図示
% gnuplot
gnuplot> plot "str_1.cat" u 2:($9*0.16), 0.6 lc 3
これを見ると暗い側の星がPSFを悪い側へ引っ張っている事が分かる
なのでこの場合は使用する星の数を200から100に減らしてもう一度est_seeing.clを回すなどした方がいい
mean,midpt,modeのどれを信じた方がいいかはケースバイケース。自分は星を100個にした場合の最頻値：0.57''を信じた。
7.WCS貼付け
既に貼付けてあると思うがWCS精度の確認も含めて。
◎USNOBカタログの入手
・USNO Archiveから自分のデータと重複するようにカタログを拾ってくる。
自分は「id,RA,Dec,sigmaRA,sigmaDEC,sRaEp,sDeEp,MuRa,MuDEC,flag,B1,B1flag,R1,R1flag,B2,B2flag,R2,R2flag,I2,I2flag」をとってきた。ファイル名はfchAYcFcm_usnob
・必要情報だけ抜粋
固有速度が10年で0.1''以下のものだけを選ぶ。それ以外は領域やflagで数を減らす程度
% sed 's/#/# /g' fchAYcFcm_usnob | awk '$1!="#" && sqrt($8*$8+$9*$9)<10 && $2>52.89583 && $2<53.32083 && $3>-27.82361 && $3<-27.54917 && $10==0{print $2,$3,$8,$9,$13,$14}' > USNOB_ECDFS.cat
% awk '{print $1,$2}' USNOB_ECDFS.cat > USNOB_ECDFS_radec.reg
◎自分の画像から星らしきものだけカタログ化
・SExtractor天体検出
以下のようにコンフィグレーションファイル、パラメータファイルを組む.
strdetect.sex、strdetect.param
% sex ECDFS_coadd_r_v2.fits -c strdetect.sex
・領域、flag、FWHM、ellipticity、CLASS_STARで星らしきものだけ抽出
% awk '$1!="#" && $8>52.89583 && $8<53.32083 && $9>-27.82361 && $9<-27.54917 && $12==0 && $10<0.1 && $13>0.9 && $11*0.16>0.55 && $11*0.16<0.65 && $2<22{print $6,$7,$8,$9}' test.cat > megaR_str.cat
これで自分の画像中の星の「x,y,ra,dec」ファイルができた。
x,yデータだけをregionファイル化
% awk '{print $1,$2}' megaR_str.cat > megaR_str_xy.reg
◎現在のmegacam画像の星ファイルとUSNOBカタログとのWCSのカタログマッチング
megaR_str.catとUSNOB_ECDFS.catとのカタログマッチングから、確からしい星の「画像中でのx座標、画像中でのy座標、USNOBでのRA、USNOBでのDec」カタログを作るねらい
・マッチング
matchcat_mega-usnob.cというプログラムを作り実行。ラフにWCSは貼られているはずなので、matching radiusは2''で十分のはず
% gcc -lm -o matchcat_mega-usnob matchcat_mega-usnob.c
% ./matchcat_mega-usnob
matching radius = 2.000000('')
MEGACAM catalog have 432 stars
USNOB catalog have 664 objects
the number of USNOB stars having 0 MEGACAM counterpart = 581
the number of USNOB stars having 1 MEGACAM counterpart = 83
the number of USNOB stars having >2 MEGACAM counterpart = 0
マッチングした星の数が100近くあれば十分
出力のMEGA_imRaDec.catは「自分画像x、自分画像y、USNOカタログRA、USNOカタログDEC」となっている。
・現在のWCSのズレの確認
matchcat_mega-usnob.cでcheck_MEGA_RaDec.catというチェック用ファイルも出力されており、それの7列目がoffset('')になっているので、offsetの統計を見る。
% awk '$1!="#"{print $7}' check_MEGA_RaDec.cat > MEGA_USNO_offset.data
cl> rtext MEGA_USNO_offset.data MEGA_USNO_offset.fits otype="real" header- pix+ nsk=0 dim=83
cl> imstat MEGA_USNO_offset.fits fields="image,npix,mean,midpt,mode,stddev,min,max"
これで出力されるoffsetの平均値や最大値は原理的には0.1''を超えてはいけない。もちろん装置による重心検出誤差はあるが、とにかくココを見て、この先自分でWCSをはるかどうか決断する。
◎WCS貼付け
・貼付ける用の画像を準備
cl> imcopy ECDFS_coadd_r_v2.fits wECDFS_coadd_r_v2.fits
・(x,y)=>(RA,Dec)への変換規則作り
cl> ccmap MEGA_imRaDec.cat megaR_WCSsol.db images=wECDFS_coadd_r_v2.fits xcol=1 ycol=2 lngcol=3 latcol=4 lngunit=degree latunit=degree project=tan fitgeom=general functio=polynomial xxorder=2 xyorder=2 yxorder=2 yyorder=2 inter-
基本的には線形変換になるはずなので、関数の次数は一次。
・貼付け
cl> ccsetwcs wECDFS_coadd_r_v2.fits megaR_WCSsol.db wECDFS_coadd_r_v2.fits project=tan
・exposure mapも同様に貼付け
cl> ccmap MEGA_imRaDec.cat megaR_WCSsol.cov.db images=wECDFS_coadd_r_v2.cov.fits xcol=1 ycol=2 lngcol=3 latcol=4 lngunit=degree latunit=degree project=tan fitgeom=general functio=polynomial xxorder=2 xyorder=2 yxorder=2 yyorder=2 inter-
cl> ccsetwcs wECDFS_coadd_r_v2.cov.fits megaR_WCSsol.db wECDFS_coadd_r_v2.cov.fits project=tan
◎変換後のWCSズレの確認
WCS貼付けられた画像に対して、再度、USNOカタログとのマッチング=>ズレの統計確認を行う
% sex wECDFS_coadd_r_v2.fits -c strdetect.sex
% awk '$1!="#" && $8>52.89583 && $8<53.32083 && $9>-27.82361 && $9<-27.54917 && $12==0 && $10<0.1 && $13>0.9 && $11*0.16>0.55 && $11*0.16<0.65 && $2<22{print $6,$7,$8,$9}' test.cat > wmegaR_str.cat
% gcc -lm -o matchcat_mega-usnob matchcat_mega-usnob.c
% ./matchcat_mega-usnob
% awk '$1!="#"{print $7}' check_wMEGA_RaDec.cat > wMEGA_USNO_offset.data
cl> rtext wMEGA_USNO_offset.data wMEGA_USNO_offset.fits otype="real" header- pix+ nsk=0 dim=85
cl> imstat wMEGA_USNO_offset.fits fields="image,npix,mean,midpt,mode,stddev,min,max"
以前に比べて良くなったかどうか確認
8.ゼロ点決定(flux calibration)
9.限界等級決め