HSC Pipeline バイナリ配布所で配布しているパイプラインの 具体的な使い方を説明していきます。 ここを読む前にまずはインストールを済ませてください。
データを置くフォルダを好きな場所に作ります。ここでは仮に ~/HSC として 進みます。
[user@example ~]$ mkdir ~/HSC
[user@example ~]$ # このフォルダを初期化するためのおまじない
[user@example ~]$ echo 'lsst.obs.hsc.HscMapper' > ~/HSC/_mapper
バイナリ配布所のデータリストの中から、 チュートリアル用画像をダウンロードしてください。 ダウンロードしたら、適当なディレクトリに展開してください。
[user@example ~]$ tar xvaf rawdata.tar.xz
ダウンロードした画像を、~/HSC の中に配置します。
[user@example ~]$ # hscPipe をセットアップ [user@example ~]$ setup-hscpipe [user@example ~]$ cd rawdata/ [user@example ~]$ # 画像を移動 [user@example ~]$ hscIngestImages.py ~/HSC */*.fits --mode=move --create
↑ 移動先 (~/HSC) は絶対パスで指定しないとうまく動きません。
~/HSC の下に移動した画像は、~/HSC/registry.sqlite3 に登録されます。 どのような画像が登録されたか見てみましょう。
[user@example ~]$ sqlite3 ~/HSC/registry.sqlite3
sqlite> .header on
sqlite> SELECT visit, filter, field, count(visit)
FROM raw
GROUP BY visit, filter, field;
sqlite> .q
次のような表示が出たら成功しています：
visit|filter|field|count(visit) 905496|HSC-G|DEN_C|112 905498|HSC-G|DEN_C|112 905500|HSC-G|DEN_C|112 905502|HSC-G|DEN_C|112 905504|HSC-G|DEN_C|112 907002|HSC-G|DOMEFLAT|112 907004|HSC-G|DOMEFLAT|112 907006|HSC-G|DOMEFLAT|112 907008|HSC-G|DOMEFLAT|112 907010|HSC-G|DOMEFLAT|112
↑もし右端のカラムの数が少ないならば、データが欠けています。
DOMEFLAT が visit=907002..907010:2 (始点..終点:間隔と表記します)、 DEN_C が visit=905496..905504:2 です。
ドームフラットからフラットデータを作ります。 (ドームフラットの visit 番号については前述)
[user@example ~]$ reduceFlat.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--output=$HOME/HSC/rerun/calib_tutorial \
--id visit=907002..907010:2 \
--detrendId calibVersion=calib_tutorial \
--config isr.doBias=False \
isr.doDark=False \
isr.doFringe=False \
--do-exec --mpiexec "-n 4"
↑重要: --calib や calibVersion の値は、正規表現 [A-Za-z0-9_+-]+ に マッチするものでなければいけません。例えば、ピリオドを含んではいけません。
↑ "-n 4" は「4スレッドで」の意味です。
↑ --config の後ろに指定できる引数とそのデフォルト値は --show config で表示できます。
データができたら、それを登録します。
[user@example ~]$ genCalibRegistry.py \
--root=$HOME/HSC/CALIB_tutorial \
--camera=HSC \
--create
ここではバイアス・ダーク・フリンジを無視しましたが、 同じようなコマンドで作れます(reduceBias.py, reduceDark.py, reduceFringe.py)。
まずアストロメトリのために、astrometry_net_data を指定します。
[user@example ~]$ setup astrometry_net_data sdss-dr8
次にCCDごとの解析プログラムを走らせます。 (visit の範囲については前述)
[user@example ~]$ reduceFrames.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id visit=905496..905504:2 \
--config processCcd.isr.doBias=False \
processCcd.isr.doDark=False \
processCcd.isr.doFringe=False \
doFocus=False $(: hscPipe 3.2.2 以降) \
--do-exec --mpiexec "-n 4"
↑ --rerun=RERUN に指定した名前を使って、$HOME/HSC/rerun/RERUN というディレクトリが作成されます (今の場合 $HOME/HSC/rerun/tutorial)。 このディレクトリに出力ファイルが書きだされます。 これらのファイルはまた、以降のコマンドの入力ファイルとなります。
★結果のみかた★(工事中)
モザイクの前に "tract" を作ります。 一つの WCS で張られる領域 (つまりだいたい 1 視野ぶん) を トラクトと呼びならわします。
[user@example ~]$ makeDiscreteSkyMap.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id visit=905496..905504:2
↑これで tract 0 が作成されます。
この時点でもし latex と dvipng がシステムにインストールされていない場合、インストールしてください。 例えば redhat 系ならば
[user@example ~]$ sudo yum install dvipng texlive-latex
そしてモザイクします。
[user@example ~]$ mosaic.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id visit=905496..905504:2 ccd=0..103 tract=0 \
--config outputDiag=True \
outputDir=$HOME/HSC/rerun/tutorial/output_mosaic \
--clobber-config
一つのトラクトは複数のパッチに分かれています (トラクト一個をまるまる扱うのは大きすぎるため)。 トラクトをスタックするには、スタックするパッチを指定する必要があります。 ではパッチの情報はどこにあるのかと言いますと、 (rerun ディレクトリ)/deepCoadd/skyMap.pickle に入っています。
実際に見てみましょう。
[user@example ~]$ cd $HOME/HSC/rerun/tutorial/deepCoadd
[user@example deepCoadd]$ python
>>> import pickle
>>> tracts = pickle.load(open("skyMap.pickle", "rb"))
>>> tract = tracts[0] # tract 番号 0 を選択
>>> print tract.getNumPatches()
↑この結果、下のように表示されるはずです。
(12, 12)
↑これは、12列 × 12行のパッチが存在する、という意味です(下図)
| (0,0) | (1,0) | ⋯ | (11,0) |
| (0,1) | (1,1) | ⋯ | (11,1) |
| ⋮ | ⋮ | ⋱ | ⋮ |
| (0,11) | (1,11) | ⋯ | (11,11) |
そのことがわかったら、もとのディレクトリに戻りましょう。
>>> exit() [user@example deepCoadd]$ cd ~ [user@example ~]$
さて閑話休題です。 スタックしてパッチ (0,0) を作りたければ patch=0,0 と指定することになります。 パッチ (0,0) と (0,1) を各々作りたければ、patch=0,0^0,1 と ^ で繋げて指定します。 一つのトラクトの中の、全てのパッチを作りたければ patch=0,0^0,1^0,2.... と書いていくわけです。
しかしこれを手で書くのはやりきれないですから、次のようにします。
[user@example ~]$ # 現在のフィールド分割子を保存 [user@example ~]$ # (IFS は普通は空白文字ですから、必ず引用符で囲んでください) [user@example ~]$ oldIFS="$IFS" [user@example ~]$ # (0,0 0,1 0,2 ... 11,11) と全て列挙するのと同じ [user@example ~]$ patches=({0..11},{0..11}) [user@example ~]$ # フィールド分割子を ^ にして [user@example ~]$ IFS='^' [user@example ~]$ # そのフィールド分割子で配列を結合する [user@example ~]$ patches="${patches[*]}" [user@example ~]$ # フィールド分割子を元に戻す [user@example ~]$ IFS="$oldIFS" [user@example ~]$ # (IFS は必ず元に戻してください)
↑この結果、$patches は次のようになっているはずです。
[user@example ~]$ echo $patches 0,0^0,1^0,2^0,3^(省略)^11,10^11,11
これを使ってスタックコマンドを発行します。
[user@example ~]$ stack.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id filter=HSC-G tract=0 patch=$patches \
--selectId visit=905496..905504:2 ccd=0..103 \
--config \
doBackgroundReference=False \
makeCoaddTempExp.bgSubtracted=True \
assembleCoadd.doMatchBackgrounds=False \
--clobber-config \
--do-exec --mpiexec "-n 4"
スタック画像の解析は forcedPhotCoadd.py で行うのですが、 現在、たくさんのパッチを一度に突っ込むと タイムアウトしてしまう不具合があります。 ですので patch=$patches と書く代わりに for ループを回し、一つ一つ処理します。
[user@example ~]$ for i in {0..11},{0..11}
do
forcedPhotCoadd.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id filter=HSC-G tract=0 patch=$i \
--config references.filter="HSC-G"
done
これも、visit=905496..905504:2 のように全部突っ込むとタイムアウトします。 ですので for ループを回して一つずつ処理します。
[user@example ~]$ # 905496..905504:2 ではないので注意
[user@example ~]$ for visit in {905496..905504..2}
do
forcedPhotCcd.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id visit=$visit ccd=0..103 tract=0 \
--config references.filter="HSC-G"
done
[user@example ~]$ calibrateCatalog.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id visit=905496..905504:2 ccd=0..103 tract=0 \
--config doApplyCalib=False \
-j 4 # 4 スレッドで
[user@example ~]$ calibrateExposure.py $HOME/HSC \
--calib=$HOME/HSC/CALIB_tutorial \
--rerun=tutorial \
--id visit=905496..905504:2 ccd=0..103 tract=0 \
-j 4 # 4 スレッドで
以上で終了です。
例えばパッチごとに細切れになった絵ではなく、全体の大きな絵が欲しい、と思ったとしてください。 hscPipe は Python 用ライブラリの巨大な集合体ですから、 Python スクリプトでそのライブラリを用いることができます。
パッチを貼り合わせて大きな絵を作るスクリプトは、例えば下のようになります。
#!/usr/bin/env python # encoding: utf-8 # AFW (Application FrameWork) のうち画像に関するライブラリ import lsst.afw.image as afwImage # AFW のうちジオメトリに関するライブラリ import lsst.afw.geom as afwGeom import sys def main(): """ パッチを並べて大きな画像を作り、"large.fits" に保存する make-large-image PATCHES... (e.g. make-large-image *,*.fits) """ if len(sys.argv) <= 1: print "一つ以上のパッチを指定してください" return 1 largeImage = makeLargeImageFromPatches(sys.argv[1:]) print "保存中 (しばらくお待ちください)..." largeImage.writeFits("large.fits") def makeLargeImageFromPatches(filenames): """ @param filenames ファイル名のリスト @return afwImage.ExposureF """ # 全ての画像をロードするとメモリ的な意味でヤバイので # とりあえずメトリックの情報だけ取得する bboxes = [] for filename in filenames: print "バウンディングボックスを取得:", filename, "..." bboxes.append( # 親 (= トラクト) におけるこのパッチの # バウンディング・ボックスを取得 afwImage.ExposureF(filename).getBBox(afwImage.PARENT) ) # 全ての BBox を含むような BBox を求める # getMin/Max 系のメソッドの戻り値はインクルーシヴ # (領域 = [min, max]) cf. getBegin/End (下で使用) x0 = min(bbox.getMinX() for bbox in bboxes) y0 = min(bbox.getMinY() for bbox in bboxes) x1 = max(bbox.getMaxX() for bbox in bboxes) y1 = max(bbox.getMaxY() for bbox in bboxes) width = x1 - x0 + 1 height = y1 - y0 + 1 # 大きな画像を用意する print "画像を用意: ({}, {})".format(width, height) exposure = afwImage.ExposureF(width, height) # トラクト内における画像の位置を設定 exposure.setXY0(afwGeom.Point2I(x0, y0)) # Exposure から画像部分だけ抜き出す largeImage = exposure.getMaskedImage() # largeImage にパッチを貼り付けていく for filename in filenames: print "パッチを取得:", filename, "..." # 再びパッチをロード patch = afwImage.ExposureF(filename) bbox = patch.getBBox(afwImage.PARENT) # このパッチの、largeImage の中における位置 # getBegin/End 系のメソッドの戻り値は one past the end # (領域 = [begin,end) ) cf. getMin/Max (上で使用) xp0 = bbox.getBeginX() - x0 yp0 = bbox.getBeginY() - y0 xp1 = bbox.getEndX () - x0 yp1 = bbox.getEndY () - y0 # パッチを張りつける ([]の中の記法については numpy のドキュメントなどを参考に) largeImage[xp0:xp1, yp0:yp1] = patch.getMaskedImage()[:] # 全てのパッチが同じ WCS その他の情報をもっているものだと盲信して # どれか一個を大きな画像にセットする if not exposure.hasWcs() and patch.hasWcs(): exposure.setCalib (patch.getCalib ()) exposure.setFilter (patch.getFilter ()) exposure.setMetadata(patch.getMetadata()) exposure.setWcs (patch.getWcs ()) return exposure if __name__ == "__main__": main()
↑のスクリプトを make-large-image.py として保存したとすると、 次のようにして使うことができます。
[user@example ~]$ setup-hscpipe
[user@example ~]$ python make-large-image.py \
$HOME/HSC/rerun/tutorial/deepCoadd/HSC-G/0/*.fits
↑16GB以上メモリを積んでいるマシンで実行してください。 そうでないと Linux が生きることをあきらめてしまうことがあります。 出来上がった画像を閲覧するためにもそれくらいメモリが必要です。