チュートリアル

English

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 rawdata]$ # 画像を移動
[user@example rawdata]$ 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 です。

sqlite3 コマンドについてさらに知りたい人は、 SQLiteとかSQLとかでググってください。

較正用データを作成

フラットデータの作成

ドームフラットからフラットデータを作ります。 (ドームフラットの visit 番号については前述)

[user@example ~]$ reduceFlat.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=calib_tutorial \
    --id visit=907002..907010:2 \
    --detrendId calibVersion=calib_tutorial \
    --config isr.doBias=False \
        isr.doDark=False \
        isr.doFringe=False \
    --batch-type=smp --cores=4

↑重要: --calib や calibVersion の値は、正規表現 [A-Za-z0-9_+-]+ に マッチするものでなければいけません。例えば、ピリオドを含んではいけません。

↑ "--cores=4" は「4スレッドで」の意味です。スレッドの数は、搭載メモリを 4GB で割った値 (と搭載CPUコア数とのうち小さい方) を目安にしてください。

--rerun=... には rerun 名(任意)を指定します。 データ解析は、Aを処理してBができ、Bを処理してCができ、…と 順々に進んでいきます。この一連の流れ A→B→C→… をrerun と呼びます。

その他のオプションの意味については、コマンドラインオプションを ご覧ください。

平均化されたフラットデータは --calib= に指定したディレクトリの下にできています。具体的には

[user@example ~]$ ls $HOME/HSC/CALIB_tutorial/FLAT/2014-02-03/HSC-G/calib_tutorial
FLAT-000.fits  FLAT-023.fits  FLAT-046.fits  FLAT-069.fits  FLAT-092.fits
FLAT-001.fits  FLAT-024.fits  FLAT-047.fits  FLAT-070.fits  FLAT-093.fits
FLAT-002.fits  FLAT-025.fits  FLAT-048.fits  FLAT-071.fits  FLAT-094.fits
    :              :              :              :              :

本来なら FITS 画像を一つ一つ閲覧すべきでしょうが、 ここでは CCD 112 枚分そろっているかだけ最低限見ておきます

[user@example ~]$ ls $HOME/HSC/CALIB_tutorial/FLAT/2014-02-03/HSC-G/calib_tutorial/FLAT-*.fits |
    wc -l
112

また、$HOME/HSC/rerun/(rerun名前) には、この最終出力以外の、中間ファイルが書きだされています。 興味があればご覧ください。

データができていたら、その情報をレジストリに登録します。

[user@example ~]$ genCalibRegistry.py \
    --root=$HOME/HSC/CALIB_tutorial \
    --camera=HSC \
    --create

↑この作業を忘れると、以降の工程は何一つうまく行きませんのでご注意ください。

フラット作成はこれで終了です。ここではバイアス・ダーク・フリンジを無視しましたが、 同じようなコマンドで作れます(reduceBias.py, reduceDark.py, reduceFringe.py)。

CCDごとの解析

まずアストロメトリのために、astrometry_net_data を指定します。

[user@example ~]$ setup astrometry_net_data sdss-dr9-fink-v5b

次に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 \
    --batch-type=smp --cores=4
出力ファイル一覧

↑ もしこのサイトで配布されているバイナリをご使用で、pthread_create が うんちゃらかんちゃらというエラーメッセージがでて動作しない場合は、 OMP_NUM_THREADS=1 reduceFrames.py $HOME/HSC... のようにして起動してみてください。

Flat を作るときには --rerun=calib_tutorial でした。今回は --rerun=tutorial を指定しています。 較正用データを作るための rerun (ひと繋がりの解析の流れ) を断ち切って、データ処理のための新しい rerun を始めよう、ということです。

その他のオプションの意味については、コマンドラインオプションを ご覧ください。

終了したとき、$HOME/HSC/rerun/(rerun名) の下にファイルが出力されています。 画像は↓

[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 枚だけなのです。ですからこれで問題ありません。

スタック

スカイマップを作る

モザイクの前に "tract" を作ります。 天球をさいの目にぶつ切りにした各々の領域を tract と言います。 一つの tract は一つの平面に投影されます。 (tract が違えば WCS も違うということです)

広い領域のサーベイではこれで良いのですが、 一視野くらいの領域をディザーして撮っただけのデータから モザイクを作るには、単に天球をぶつ切りにした tract を使うのは不便です。

そこで、モザイクしたい手持ちの画像を全て含むような tract を 特別に定義します。それが Discrete Sky Map です。

[user@example ~]$ makeDiscreteSkyMap.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id visit=905496..905504:2
出力ファイル一覧

↑これで最後に

2015-04-08T06:32:05: makeDiscreteSkyMap: tract 0 has corners (134.913, 23.031), (132.683, 23.031), (132.665, 25.084), (134.931, 25.084) (RA, Dec deg) and 11 x 11 patches

などと表示されればOKです。トラクト 0 番が定義されました。

モザイクする

トラクト 0 番をモザイクします。 モザイクは座標を決めるだけで実際に画像を作りはしません。

[user@example ~]$ mosaic.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id visit=905496..905504:2 ccd=0..103 tract=0
出力ファイル一覧

画像を較正しなおす

モザイクによって、各CCD画像の等級ゼロ点やWCSなどが決まります。 これを各CCD画像に反映します(較正されたCCD画像が不要ならスキップ可能です)。

[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 スレッドで
出力ファイル一覧

↑この際、(視野の端のほうにある)「wcsファイルが開けない」という警告がでるかもしれません。 端の方がモザイクに失敗して wcs ができていないということです。 無視してかまいません。

較正されたCCD画像は↓にできます

[user@example ~]$ ls $HOME/HSC/rerun/tutorial/00762/HSC-G/corr/0000/CALEXP-*.fits

カタログを較正しなおす

同様に、各カタログにモザイクの結果を反映します(較正されたカタログが不要ならスキップ可能です)。

[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 ~]$ ls $HOME/HSC/rerun/tutorial/00762/HSC-G/output/0000/CALSRC-*.fits

スタックする

モザイクが終わった時点では、 まだ実際にスタックされた画像はできていませんので、 ここでスタックして Coadd画像を作ります。

[user@example ~]$ stack.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id filter=HSC-G tract=0 \
    --selectId visit=905496..905504:2 ccd=0..103 \
    --batch-type=smp --cores=4
出力ファイル一覧

Coadd画像は↓にできます。

[user@example ~]$ cd $HOME/HSC/rerun/tutorial/deepCoadd/HSC-G/0
[user@example 0]$ ls */calexp-*.fits
0,4/calexp-HSC-G-0-0,4.fits    5,1/calexp-HSC-G-0-5,1.fits
0,5/calexp-HSC-G-0-0,5.fits    5,10/calexp-HSC-G-0-5,10.fits
0,6/calexp-HSC-G-0-0,6.fits    5,2/calexp-HSC-G-0-5,2.fits
            :                              :

トラクト0番を一枚スタックしただけなのに、こんなにたくさんファイルができています。 これはなぜかというと、トラクト一枚が広すぎてメモリに載せるのも一苦労ですので 細かくパッチに分割してあるのです。一つの大きな画像を分割してあるだけですので、 例えば ds9 を使えば、下のようにして 全画像を一度に繋げて表示することもできます(ただし16GBを越えるメモリを消費します)。 

[user@example 0]$ ds9 -mosaic wcs */calexp-*.fits

トラクト全体はこんなふうに見えるはずです: stacked image (8MB; 102倍に縮小)

マルチバンド解析

チュートリアルデータは 1色のみなのですが、これはチュートリアルですので 無理やりマルチバンド解析に進みます。 各々のバンドで作った天体カタログを一つにマージし、さらにそれを測りなおします。

[user@example ~]$ multiBand.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --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

スタック画像をレファレンスにして各CCDを測りなおす

スタック画像の上で検出された天体の位置 (正確には複数バンドのカタログをマージしたもの) を、スタックする前の各々の CCD 画像に立ち戻って測定します。

[user@example ~]$ forcedCcd.py $HOME/HSC \
    --calib=$HOME/HSC/CALIB_tutorial \
    --rerun=tutorial \
    --id visit=905496..905504:2 ccd=0..103 tract=0 \
    --batch-type=smp --cores=4
出力ファイル一覧

結果のカタログは↓に書きだされます

[user@example ~]$ ls $HOME/HSC/rerun/tutorial/00762/HSC-G/tract0/FORCEDSRC-*.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 のコマンド群は、どれもだいたい共通のコマンドラインを持ちます。 reduceFrames.py を例にとりますと:

[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 \
    --batch-type=smp --cores=4

まず解説しなければいけないのが最後の部分です。 reduce** という名のコマンドは、本来、ジョブキュー (Torque PBSなど) にジョブを投げるコマンドです。 しかしそういうシステムを持っていない場合は、--batch-type=smp オプションによって ローカルのコンピュータ上で処理を実行できます。 その際のスレッド数を --cores= に指定します。

[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 \
    --batch-type=smp --cores=4

次に--configオプションです。ここで解析時のコンフィグを変えられます。 指定できるパラメータとデフォルト値を知りたい場合は、コマンドラインの後ろに --show configと付記して起動してください。

コンフィグを変更するとプログラムは動作しなくなります。

コンフィグを変更するとプログラムは動作しなくなります。

大事なことなので二回言いました。 同じ rerun の中に違う解析パラメータでできたデータが混在することは好ましくないからです。

もしプログラムが動作しなくなった場合は、 --clobber-config をコマンドラインの後ろに付記してください。 キャッシュされている過去のコンフィグを上書きします。

[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 \
    --batch-type=smp --cores=4

最後に--idについて。入力データを検索するための条件をここに入れます。 複数指定すると AND で連結されます。たとえば--id field="DEN_C" filter="HSC-G" と指定すれば、DEN_Cフィールドかつ gバンドなデータのすべてが入力になります。 "field" や "filter" などは レジストリの中のテーブルのコラム名です。

なお、visit で直接指定せずに field などで指定する方法は直感的である一方で、 仮にFITSヘッダが間違っていたりするとどのデータが入力されどのデータがされなかったのか わかりにくくなります。ご注意ください。