HSC Pipeline バイナリ配布所で配布しているパイプラインの 具体的な使い方を説明していきます。 ここを読む前にまずはインストールを済ませてください。
データを置くフォルダを好きな場所に作ります。
ここでは仮に ~/HSC として進みます。
このフォルダを作るには mkdir ではなく create-rootdir.sh をお使いください。
このスクリプトは単にフォルダを作るだけでなく、必要なファイルをその中に設置します。
[user@example ~]$ # hscPipe をセットアップ [user@example ~]$ setup-hscpipe [user@example ~]$ create-rootdir.sh ~/HSC /path/to/ps1_pv3_3pi_20170110 /path/to/STRAY_LIGHT
↑ 第二引数 ps1_pv3_3pi_20170110 はレファレンスカタログフォルダへのパスです。
このパス (REFDIR) は、 REFDIR/{number}.fits と
REFDIR/config.py が存在するような REFDIR です。
また、このフォルダ REFDIR の名前は ps1_pv3_3pi_20170110 でなくてはいけません。
パイプラインのデフォルト設定で、この名前が仮定されています。
↑ 第三引数 STRAY_LIGHT は y バンド迷光ファイルへのパスです。
このチュートリアルの入力データには y バンドのデータがありませんので、ここでは省略して構いません。
さらに、次のコマンドでトランスミッションカーブをインストールします。
[user@example ~]$ installTransmissionCurves.py ~/HSC
バイナリ配布所のデータリストの中から、 チュートリアル用画像をダウンロードしてください。 ダウンロードしたら、適当なディレクトリに展開してください。
[user@example ~]$ tar xvaf rawdata.tar.xz
ダウンロードした画像を、~/HSC の中に配置します。
[user@example ~]$ ## まだのときはここで setup-hscpipe してください [user@example ~]$ # setup-hscpipe [user@example ~]$ cd rawdata/ [user@example rawdata]$ # 画像を移動 [user@example rawdata]$ hscIngestImages.py ~/HSC */*.fits --mode=move
↑ 移動先 (~/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 です。
sqlite3 コマンドについてさらに知りたい人は、 SQLiteとかSQLとかでググってください。
ドームフラットからフラットデータを作ります。 (ドームフラットの visit 番号については前述)
[user@example ~]$ constructFlat.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id visit=907002..907010:2 $(: もしくは --id field=DOMEFLAT filter=HSC-G) \ --config \ isr.doBias=False \ isr.doDark=False \ --batch-type=smp --cores=4
↑ "--cores=4" は「4スレッドで」の意味です。スレッドの数は、搭載メモリを 4GB で割った値 (と搭載CPUコア数とのうち小さい方) を目安にしてください。
--rerun=... には rerun 名(任意)を指定します。
データ解析は、Aを処理してBができ、Bを処理してCができ、…と
順々に進んでいきます。この一連の流れ A→B→C→… をrerun と呼びます。
その他のオプションの意味については、コマンドラインオプションを ご覧ください。
平均化されたフラットデータは --rerun= に指定したディレクトリの下にできています。具体的には
[user@example ~]$ ls $HOME/HSC/rerun/tutorial/FLAT/2014-02-03/HSC-G FLAT-2014-02-03-HSC-G-000.fits FLAT-2014-02-03-HSC-G-056.fits FLAT-2014-02-03-HSC-G-001.fits FLAT-2014-02-03-HSC-G-057.fits FLAT-2014-02-03-HSC-G-002.fits FLAT-2014-02-03-HSC-G-058.fits : :
本来なら FITS 画像を一つ一つ閲覧すべきでしょうが、 ここでは CCD 112 枚分そろっているかだけ最低限見ておきます
[user@example ~]$ ls $HOME/HSC/rerun/tutorial/FLAT/2014-02-03/HSC-G/FLAT-*.fits | wc -l 112
また、$HOME/HSC/rerun/tutorial には、この最終出力以外の中間ファイルも書きだされています。
興味があればご覧ください。
フラットを作ったら、登録しなくてはいけません。 以降のステージで解析パイプラインがそのフラットを見つけられるようにです。
[user@example ~]$ ingestCalibs.py $HOME/HSC \
@<(find $HOME/HSC/rerun/tutorial/FLAT -name "*.fits") \
--calib=$HOME/HSC/CALIB \
--validity=100 \
--mode=link
--validity=100 とは、フラットの日付の前後100日間に撮られた撮像データに対して
このフラットが使用可能であることを意味します。
フラット作成はこれで終了です。ここではバイアス・ダーク・フリンジを無視しましたが、
同じようなコマンドで作れます(constructBias.py, constructDark.py, constructFringe.py)。
大きな銀河の周りでは、スカイのレベルを高く見積りがちです。 高く見積もりすぎたスカイを画像から引くと、このような銀河の周りに人工的な堀ができます。 スタック時にこの堀を埋めるのですが、そのときにスカイレベルを非常に大きなメッシュで 見積もりなおす必要があります。このときに必要になる、平均のスカイパターンを事前に作っておきます。 このスカイパターンを「グローバルスカイパターン」と呼んでいます。
グローバルスカイパターンは、オブジェクトフレーム (実際の星空を写した画像) から作ります。 以下のプログラムを実行してください。 (visit の範囲については前述)
[user@example ~]$ constructSky.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id visit=905496..905504:2 $(: もしくは --id field=DEN_C)\ --config \ isr.doBias=False \ isr.doDark=False \ --batch-type smp --cores=4
スカイを作ったら、登録しなくてはいけません。 以降のステージで解析パイプラインがそのスカイを見つけられるようにです。
[user@example ~]$ ingestCalibs.py $HOME/HSC \
@<(find $HOME/HSC/rerun/tutorial/SKY -name "*.fits") \
--calib=$HOME/HSC/CALIB \
--validity=100 \
--mode=link
古いバージョンのパイプラインをご存知の方へ:
レファレンスカタログをセットアップする必要はなくなりました。
create-rootdir.sh がすでにレファレンスカタログへのシンボリックリンクを作ってくれてあるからです。
CCD画像ごとの解析プログラムを走らせます。 (visit の範囲については前述)
[user@example ~]$ singleFrameDriver.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id visit=905496..905504:2 $(: もしくは --id field=DEN_C)\ --config \ processCcd.isr.doBias=False \ processCcd.isr.doDark=False \ --batch-type=smp --cores=4
↑ もしこのサイトで配布されているバイナリをご使用で、pthread_create が
うんちゃらかんちゃらというエラーメッセージがでて動作しない場合は、
env OMP_NUM_THREADS=1 singleFrameDriver.py $HOME/HSC...
のようにして起動してみてください。
オプションの意味については、コマンドラインオプションを ご覧ください。
終了したとき、$HOME/HSC/rerun/tutorial の下にファイルが出力されています。
画像は↓
[user@example ~]$ cd $HOME/HSC/rerun/tutorial/00762/HSC-G/corr/ [user@example corr]$ ls CORR-* CORR-0905496-000.fits CORR-0905498-070.fits CORR-0905502-036.fits CORR-0905496-001.fits CORR-0905498-071.fits CORR-0905502-037.fits CORR-0905496-002.fits CORR-0905498-072.fits CORR-0905502-038.fits : : :
天体カタログ (FITS BinTable) は↓
[user@example ~]$ cd $HOME/HSC/rerun/tutorial/00762/HSC-G/output/ [user@example output]$ ls SRC-* SRC-0905496-000.fits SRC-0905498-070.fits SRC-0905502-036.fits SRC-0905496-001.fits SRC-0905498-071.fits SRC-0905502-037.fits SRC-0905496-002.fits SRC-0905498-072.fits SRC-0905502-038.fits : : :
これも実際に目視してみたほうがよいですが、 ここではとりあえず個数だけ確かめておきます。
[user@example ~]$ ls $HOME/HSC/rerun/tutorial/00762/HSC-G/output/SRC-*.fits | wc -l 520
520 = 5×104 枚しかありません。
1ショットあたりのCCD数は、フォーカス合わせ用のものを除けば 112 ではなく 104 枚だけなのです。ですからこれで問題ありません。
ここでスカイレベルを見積もりなおします。 このセクションで説明したとおりです。
[user@example ~]$ skyCorrection.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id visit=905496..905504:2 $(: もしくは --id field=DEN_C)\ --batch-type=smp --cores=4
このコマンドからの出力が、スタックの入力に使われます。
モザイクの前に "tract" を作ります。 天球をさいの目にぶつ切りにした各々の領域を tract と言います。 一つの tract は一つの平面に投影されます (tract が違えば WCS も違うということです)。 また tract の集合を "sky map" と呼びます。
広い領域のサーベイではこれで良いのですが、 一視野くらいの領域をディザーして撮っただけのデータから モザイクを作るには、単に天球をぶつ切りにした tract を使うのは不便です。 一つの領域がいくつかのトラクトに分かれてしまう可能性が高いからです。
そこで、モザイクしたい手持ちの画像を全て含むような tract を 特別に定義します。このようなトラクトから成る sky map を "Discrete sky map" と呼びます。
[user@example ~]$ makeDiscreteSkyMap.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id visit=905496..905504:2 $(: もしくは --id field=DEN_C)
↑これで最後に
makeDiscreteSkyMap INFO: tract 0 has corners (134.970, 22.961), (132.538, 22.961), (132.517, 25.201), (134.992, 25.201) (RA, Dec deg) and 12 x 12 patches
などと表示されればOKです。トラクト 0 番が定義されました。
トラクト 0 番をモザイクします。 モザイクは座標を決めるだけで実際に画像を作りはしません。
[user@example ~]$ mosaic.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id visit=905496..905504:2 tract=0 $(: or --id field=DEN_C filter=HSC-G tract=0)
モザイクによって、各CCD画像の等級ゼロ点やWCSなどが決まります。 これを各CCD画像に反映します(較正されたCCD画像が不要ならスキップ可能です)。
[user@example ~]$ calibrateExposure.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id visit=905496..905504:2 tract=0 $(: or --id field=DEN_C tract=0) \ -j 4 # 4 スレッドで
↑この際、“NoResults: No locations for get (視野の端のほうにあるデータ)” という警告がでるかもしれません。 端の方がモザイクに失敗して位置が決まっていないということです。 無視してかまいません。
較正されたCCD画像は↓にできます
[user@example ~]$ ls $HOME/HSC/rerun/tutorial/00762/HSC-G/corr/0000/CALEXP-*.fits
同様に、各カタログにモザイクの結果を反映します(較正されたカタログが不要ならスキップ可能です)。
[user@example ~]$ calibrateCatalog.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id visit=905496..905504:2 tract=0 $(: or --id field=DEN_C tract=0) \ --config doApplyCalib=False \ -j 4 # 4 スレッドで
較正されたカタログは↓にできます
[user@example ~]$ ls $HOME/HSC/rerun/tutorial/00762/HSC-G/output/0000/CALSRC-*.fits
モザイクが終わった時点では、 まだ実際にスタックされた画像はできていませんので、 ここでスタックして Coadd画像を作ります。
[user@example ~]$ coaddDriver.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id filter=HSC-G tract=0 \ --selectId visit=905496..905504:2 ccd=0..103 $(: もしくは --selectId field=DEN_C filter=HSC-G) \ --batch-type=smp --cores=4
Coadd画像は↓にできます。
[user@example ~]$ cd $HOME/HSC/rerun/tutorial/deepCoadd-results/HSC-G/0 [user@example 0]$ ls */calexp-*.fits 1,3/calexp-HSC-G-0-1,3.fits 5,3/calexp-HSC-G-0-5,3.fits 1,4/calexp-HSC-G-0-1,4.fits 5,4/calexp-HSC-G-0-5,4.fits 1,5/calexp-HSC-G-0-1,5.fits 5,5/calexp-HSC-G-0-5,5.fits : :
トラクト0番を一枚スタックしただけなのに、こんなにたくさんファイルができています。 これはなぜかというと、トラクト一枚が広すぎてメモリに載せるのも一苦労ですので 細かくパッチに分割してあるのです。一つの大きな画像を分割してあるだけですので、 例えば ds9 を使えば、下のようにして 全画像を一度に繋げて表示することもできます(ただし16GBを越えるメモリを消費します)。
[user@example 0]$ ds9 -mosaic wcs */calexp-*.fits
トラクト全体はこんなふうに見えるはずです: stacked image (8MB; 102倍に縮小)
"det-*.fits" というファイルが画像 "calexp-*.fits"
の横にあることに気づかれると思います。 ↓
[user@example ~]$ cd $HOME/HSC/rerun/tutorial/deepCoadd-results/HSC-G/0 [user@example 0]$ ls */det-*.fits 1,3/det-HSC-G-0-1,3.fits 5,3/det-HSC-G-0-5,3.fits 1,4/det-HSC-G-0-1,4.fits 5,4/det-HSC-G-0-5,4.fits 1,5/det-HSC-G-0-1,5.fits 5,5/det-HSC-G-0-5,5.fits : :
これらのファイルには、スタック画像の中の天体の「検出」情報が入っています。 まだ測定はされていません。 測定は次の「マルチバンド解析」で行われます。
チュートリアルデータは 1色のみなのですが、これはチュートリアルですので 無理やりマルチバンド解析に進みます。 各々のバンドで作った天体の検出情報を一つにマージし、さらにそれを測ります。
[user@example ~]$ multiBand.py $HOME/HSC \ --rerun=tutorial \ --calib=$HOME/HSC/CALIB \ --id tract=0 filter=HSC-G \ $(: 他のバンドもある場合 filter=HSC-G^HSC-I^HSC-Z) \ --batch-type=smp --cores=4
マージされた共通のカタログを用いて各々のバンドを測定した結果は↓
[user@example ~]$ ls $HOME/HSC/rerun/tutorial/deepCoadd-results/HSC-G/0/*/meas-*.fits
天体の中心位置などを共通の値に固定して各々のバンドを測定した結果は↓
[user@example ~]$ ls $HOME/HSC/rerun/tutorial/deepCoadd-results/HSC-G/0/*/forced_src-*.fits
以上で終了です。
例えばパッチごとに細切れになった絵ではなく、全体の大きな絵が欲しい、と思ったとしてください。 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-results/HSC-G/0/*/calexp-*.fits
↑16GB以上メモリを積んでいるマシンで実行してください。 そうでないと Linux が生きることをあきらめてしまうことがあります。 出来上がった画像を閲覧するためにもそれくらいメモリが必要です。
注意: ↑上のスクリプトはマスクレイヤーの扱いに難があります。 マスク値の各ビットに振られる意味(「サチった」「宇宙線が入った」など) は パッチによって異な(ることがあ)ります。 どのビットがどういう意味を持つのかは FITS ヘッダ (EXTNAME=MASK のエクステンション) に記述されていますので、 それを読んでビットを入れ替えるなどの操作が本来は必要です。 (チュートリアルの範囲を超えているので割愛します)
hscPipe のコマンド群は、どれもだいたい共通のコマンドラインを持ちます。
singleFrameDriver.py を例にとりますと:
[user@example ~]$ singleFrameDriver.py $HOME/HSC \
--rerun=tutorial \
--calib=$HOME/HSC/CALIB \
--id visit=905496..905504:2 \
--config \
processCcd.isr.doBias=False \
processCcd.isr.doDark=False \
--batch-type=smp --cores=4
まず解説しなければいけないのが最後の部分です。
**Driver.py という名のコマンドは、本来、ジョブキュー (Torque PBSなど) にジョブを投げるコマンドです。
しかしそういうシステムを持っていない場合は、--batch-type=smp オプションによって
ローカルのコンピュータ上で処理を実行できます。
その際のスレッド数を --cores= に指定します。
[user@example ~]$ singleFrameDriver.py $HOME/HSC \
--rerun=tutorial \
--calib=$HOME/HSC/CALIB \
--id visit=905496..905504:2 \
--config \
processCcd.isr.doBias=False \
processCcd.isr.doDark=False \
--batch-type=smp --cores=4
次に--configオプションです。ここで解析時のコンフィグを変えられます。
指定できるパラメータとデフォルト値を知りたい場合は、コマンドラインの後ろに
--show configと付記して起動してください。
コンフィグを変更するとプログラムは動作しなくなります。
コンフィグを変更するとプログラムは動作しなくなります。
大事なことなので二回言いました。 同じ rerun の中に違う解析パラメータでできたデータが混在することは好ましくないからです。
もしプログラムが動作しなくなった場合は、
--clobber-config をコマンドラインの後ろに付記してください。
キャッシュされている過去のコンフィグを上書きします。
[user@example ~]$ singleFrameDriver.py $HOME/HSC \
--rerun=tutorial \
--calib=$HOME/HSC/CALIB \
--id visit=905496..905504:2 \
--config \
processCcd.isr.doBias=False \
processCcd.isr.doDark=False \
--batch-type=smp --cores=4
最後に--idについて。入力データを検索するための条件をここに入れます。
複数指定すると AND で連結されます。たとえば--id field="DEN_C" filter="HSC-G"
と指定すれば、DEN_Cフィールドかつ gバンドなデータのすべてが入力になります。
"field" や "filter" などは レジストリの中のテーブルのコラム名です。
なお、visit で直接指定せずに field などで指定する方法は直感的である一方で、 仮にFITSヘッダが間違っていたりするとどのデータが入力されどのデータがされなかったのか わかりにくくなります。ご注意ください。