2009.11/10
前回導入したSExtractorを使って,flat-fielding,sky引き,distortion補正済の複数枚のfitsファイルの位置合わせを行う。今日はとりあえず二枚の位置合わせを行う。
disMCSA38539.fitsをdisMCSA38537.fitsに合わせる感じ
・SExtractorの天体検出の仕方・・・sky値よりも?-σ以上高いcount値のpixelが?個以上つながってる部分を天体とみなして検出してくれる。?の部分は自分で設定可能。
・xyxymatch・・・画像ファイル間の同一天体を探してリスト化してくれるが、コケやすいので注意
・geomap・・・xyxymatchで作った天体対応表を元に二枚の画像内の天体が重なるような変換規則を算出してくれる。
☆其の一 天体検出
◎必要なsexファイル(パラメータファイル)つくる
・sextractorディレクトリのなかからdefaultファイルを持ってくる
ホームディレクトリ/sextractor-2.5.0/config/に入り
$ cd ~/sextractor-2.5.0/config/
$ ls
defaultから始まるファイルを作業スペースに移動
$ cp default* ~/irafwork/practice/
作業スペースに移り
$ cd ~/irafwork/practice/
$ emacs default.sex &
default.sexを右のように書き換えて1110.sexと名前変えて保存→1110.sex
ちなみに何を出力させるかのファイルdefault.paramも目的に応じて書き換えたりするが今回はdefaultのまま使う。→default.param
◎作ったパラメータファイルを適当な二枚にかける
・まずはdisMCSA38539.fitsにかける
コマンド端末上で
$ sex disMCSA38539.fits -c 1110.sex
----- SExtractor 2.5.0 started on 2010-03-01 at 16:37:37 with 1 thread
Measuring from: "CL0024" / 2048 x 2048 / 32 bits FLOATING POINT data
(M+D) Background: -9.2743 RMS: 87.2708 / Threshold: 174.542
Objects: detected 5744 / sextracted 4523
> All done (in 3 s)
で実行してくれる。
ちゃんと結果が出てるかチェックしてみる。
$ more test.cat //結果ファイルの確認//
# 1 NUMBER Running object number
# 2 FLUXERR_ISO RMS error for isophotal flux [count]
# 3 FLUX_AUTO Flux within a Kron-like elliptical aperture [count]
# 4 FLUXERR_AUTO RMS error for AUTO flux [count]
# 5 X_IMAGE Object position along x [pixel]
# 6 Y_IMAGE Object position along y [pixel]
# 7 FLAG SExtraction flags
1 3694.351 -1275657 9815.155 1809.857 118.680 3
2 2875.968 323287.6 6223.211 1805.493 48.211 19
3 2895.761 -1096936 6894.395 1897.971 54.518 19
4 3670.565 -464499.5 8764.963 1818.271 81.989 3
5 2893.13 -2101446 6904.33 1914.395 93.433 3
一列目から順に番号、Isophotal fluxのエラー,Auto flux,Auto fluxのエラー,x座標,y座標,flag(小さい方がいい？)が列挙されてる。今default.paramをそのまま使ったせいでこうなっている。
ds9でcheck.fitsをみてみる
cl> !ds9 &
cl> reset stdimage=imt4096
cl>display check.fits 1 zs- zr- z1=0.8 z2=1.2
cl>display disMCSA38539.fits 2
Tab連打で元のdisMCSA38539.fitsと比べてみる→上手く検出されてるっぽい
結果ファイルtest.catをtest38539.catという名で保存しとく。
$ mv test.cat test38539.cat
・同様にしてもう一枚disMCSA38537.fitsにかける
$ sex disMCSA38537.fits -c 1110.sex
$ mv test.cat test38537.cat
☆其の二 天体の位置合わせ
◎xyxymatchにかけるリスト作り
今天体リスト,test38539.catとtest38537.catがあるが、これをこのまま使って同一天体を探すことはできない。例えば、check.fitsを見るとへりの部分の明らかなゴミを天体として検出しているのが分かるし、また暗い天体として検出されてる天体は本当に天体なのか単なるsky noiseなのか疑わしい。なのでへりを除いた範囲で明るい天体だけのリストを作って、それをxyxymatchにかけた方が安全なのである。xyxymatchはコケやすいタスクなのでここは工夫のし所であるが、自分は以下のように行った。
$ awk '$5>150 && $5<1750 && $6>150 && $6<1750{print $1,$5,$6,$3}' test38539.cat | sort -n -k 4 -r | head -40 | tail -30 > test_luminous38539.cat
一応説明すると、test38539.catの天体リストの中から150<$5(x座標)<1750かつ150<$6(y座標)<1750である天体の番号,x,y,Auto fluxだけリスト化、更にAuto fluxの明るい順に並べかえ、明るいtop10からtop40までの30個をtest_luminous38539.catという名でリスト化したということ。
中見てみるとこんな感じ
$ more test_luminous38539.cat
3502 1099.852 1514.901 478996.5
2382 1619.069 1043.269 473133.7
3488 511.297 1524.058 460485.2
3217 1656.436 1603.225 459972.8
・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・
同様にtest_luminous38537.catというxyxymatch用天体リストを作る。
$ awk '$5>150 && $5<1750 && $6>150 && $6<1750{print $1,$5,$6,$3}' test38537.cat | sort -n -k 4 -r | head -40 | tail -30 > test_luminous38537.cat
◎同一天体の対応表を作る(xyxymatch)
cl> epar xyxymatch
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatch
TASK = xyxymatch
input = test_luminous38539.cat The input lists
referenc= test_luminous38537.cat The reference lists
output = shift1110.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= 2) Input list column containing the x coordinate //インプットの//
(ycolumn= 3) Input list column containing the y coordinate //$2,$3と//
(xrcolum= 2) Reference list column containing the x coordinat//リファレンスの$2,$3をx,y座標とみて//
(yrcolum= 3) Reference list column containing the y coordinat//対応させる//
(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= 50) The maximum number of rejection iterations //ここもdefaultから変えた//
(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)
→結果
$ more shift1110.cat
# Input: test_luminous38539.cat Reference: test_luminous38537.cat Number of ti
e points: 0
・・・・・・・・・・・・・・（略）・・・・・・・・・・・・・・・・・・
# Column definitions
# Column 1: X reference coordinate
# Column 2: Y reference coordinate
# Column 3: X input coordinate
# Column 4: Y input coordinate
# Column 5: Reference line number
# Column 6: Input line number
736.932 1515.176 1099.852 1514.901 4 1
1255.871 1043.472 1619.069 1043.269 6 2
1293.688 1603.740 1656.436 1603.225 8 4
992.239 337.174 1355.301 336.912 9 6
338.671 1579.012 701.232 1579.150 10 7
1244.448 386.024 1607.315 385.757 13 8
1307.324 1475.462 1669.822 1475.432 15 10
876.113 813.799 1239.063 813.594 16 11
365.222 715.919 728.109 715.587 18 17
1187.599 725.974 1550.757 725.676 20 19
570.399 658.062 933.88 657.743 21 21
1380.873 1100.369 1743.541 1099.966 22 22
1173.227 956.046 1536.269 955.635 27 25
ちょっと少ない気がするがこのまま行く。
◎disMCSA38539.fits→disMCSA38537.fitsの変換規則を作る(geomap)
cl> epar geomap
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatch
TASK = geomap
input = shift1110.cat The input coordinate files
database= database1110.cat The output database file //対応表shift1110.catから作る変換規則のDatabaseを指定//
xmin = 1. Minimum x reference coordinate value
xmax = 2048. Maximum x reference coordinate value
ymin = 1. Minimum y reference coordinate value
ymax = 2048. Maximum y reference coordinate value
(transfo= ) The output transform records names
(results= ) The optional results summary files
(fitgeom= general) Fitting geometry
(functio= polynomial) Surface type
(xxorder= 2) Order of x fit in x //一次変換したい場合は2にしとく？？//
(xyorder= 2) Order of x fit in y //同じく一次変換なので２？？//
(xxterms= half) X fit cross terms type
(yxorder= 2) Order of y fit in x
(yyorder= 2) Order of y fit in y
(yxterms= half) Y fit cross terms type
(maxiter= 50) Maximum number of rejection iterations
(reject = 3.) Rejection limit in sigma units
(calctyp= real) Computation type
(verbose= yes) Print messages about progress of task ?
(interac= no) Fit transformation interactively ?
(graphic= stdgraph) Default graphics device
(cursor = ) Graphics cursor
(mode = ql)
(:go)
Coordinate list: shift1110.cat Transform: shift1110.cat
Results file:
Coordinate mapping status
X fit ok. Y fit ok.
Xin and Yin fit rms: 0.2192054 0.1499203
Coordinate mapping parameters
Mean Xref and Yref: 978.662 993.0944
Mean Xin and Yin: 1341.589 992.8342
X and Y shift: 363.3746 -0.1867662 (xin yin)
X and Y scale: 0.999912 1.000093 (xin / xref yin / yref)
X and Y axis rotation: 0.010 359.979 (degrees degrees)
これでたぶん変換規則のDatabase=database1110.catが出来たはず。
◎disMCSA38539.fitsの画像幾何変換(geotran)
cl>epar geotran
I R A F
Image Reduction and Analysis Facility
PACKAGE = immatch
TASK = geotran
input = disMCSA38539.fits Input data
output = shift38539.fits Output data
database= database1110.cat Name of GEOMAP database file
transfor= shift1110.cat Names of coordinate transforms in database file //ここにはdatabase1110.catの中に書かれている変換式名shift1110.catを入れる//
(geometr= geometric) Transformation type (linear,geometric)
(xin = INDEF) X origin of input frame in pixels
(yin = INDEF) Y origin of input frame in pixels
(xshift = INDEF) X origin shift in pixels
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・（略）・・・・・・・・・・・・・・・・・・・・・・
・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・
◎結果の確認
cl> display shift38539.fits 1 →こんなの
cl> display disMCSA38537.fits 2 →こんなの
ds9上でTab連打してもらうと、天体の位置が合ってることが確認できるはず！
