HSC 解析講習会用チュートリアル

このページは 2016年 7 月 20, 21 日に行われた HSC 画像解析講習会用のチュートリアルです。実習は Pipeline を用いた解析と、最終的なデータで疑似カラー画像作成と色等級図作成を行います。解析に使うコマンドの中には今回の講習会に特化したコマンドオプションもありますので、注意してください。
HSC pipeline マニュアルHSC pipeline バイナリ配布所チュートリアル にも詳細がありますので、参照下さい。 また、上級者には こちら も参考になると思います。Pipeline についてかなり細かく書かれていますが、古いバージョンの情報もありますので、参照する際にはご注意下さい。
今回使用するデータは 3 バンド (i, z, Y)で、中心 6 CCD を 5 ショット分です。ドームフラットについても同様です。
また、HSC pipeline のバージョンは最新の 4.0.2 になります。

Pipeline のインストールとアストロメトリ用データの設置

Pipeline の動作環境を確認の上、下記ページを参考にインストールを行って下さい。

Rawdata の準備

まず、rawdata をダウンロードしましょう。

 # 個人の解析ディレクトリを作成します。
 cd /data/users
 mkdir <User Name>
 # 自分の好きなように作って下さい。

 # Rawdata をローカルディスクにダウンロード
 wget http://hsc.mtk.nao.ac.jp/HSC_Training_download/hdfn-6CCD.tar.gz -P /data/users/<User Name>

 # Rawdata を展開
cd /data/users/<User Name>
 tar xvf hdfn-6CCD.tar.gz

 # Rawdata の確認
 cd hdfn-6CCD/
 ls | wc -l
 # 結果が 180 であれば問題ありません。

 # どのような画像が撮れているのか確認してみましょう。
 ds9 xxxxx.fits
 # どれでもいいので、自分の好きなファイルを選んでください。

続いて、解析用ディレクトリを準備し、レジストリを作成します。

# HSC pipeline セットアップコマンドの登録
source /data/users/opt/hscpipe/4.0.2/bashrc

 # HSC pipeline を呼び出すためのセットアップコマンド (ログインの度に必ず行ってください)
 setup-hscpipe

 # 解析用ディレクトリ(リポジトリの親ディレクトリ)の作成
 mkdir /data/users/<User Name>/HSC

 # _mapper ファイルの作成 (フォルダを初期化するおまじないのようなもの)
 echo 'lsst.obs.hsc.HscMapper' > /data/users/<User Name>/HSC/_mapper

/data/users/<User Name>/HSC の下を確認すると、_mapper という名前のファイルが作られています。これで解析用ディレクトリの準備は完了です。

次に、rawdata を解析用ディレクトリに配置し、登録します。

# 生データのあるディレクトリに移動
cd /data/users/<User Name>/hdfn-6CCD

# 解析ディレクトリ下にレジストリを作成
hscIngestImages.py /data/users/<User Name>/HSC  HSC*.fits --mode=link --create
この後、/data/users/<User Name>/HSC の下には、DOMEFLAT、HDFN という名前のフォルダと、registry.sqlite3 というファイルが作られています。
では、レジストリに登録された画像の情報を確認しましょう。 registry.sqlite3 に登録された情報は SQL というデータベース言語を用いて検索することができます (SQL の詳細についてはここでは触れません)。
cd /data/users/<User Name>/HSC
sqlite3 /data/users/<User Name>/HSC/registry.sqlite3

sqlite> .header on
sqlite> SELECT visit, filter, field, count(visit)
        FROM raw
        GROUP BY visit, filter, field;

# 以下のように、カラムの右端 (= CCD 数) が全て 6 になっていればOK.
visit|filter|field|count(visit)
212|HSC-Y|DOMEFLAT|6
214|HSC-Y|DOMEFLAT|6
216|HSC-Y|DOMEFLAT|6
218|HSC-Y|DOMEFLAT|6
220|HSC-Y|DOMEFLAT|6
710|HSC-Y|HDFN|6
712|HSC-Y|HDFN|6
714|HSC-Y|HDFN|6
716|HSC-Y|HDFN|6
718|HSC-Y|HDFN|6
950|HSC-I|HDFN|6
952|HSC-I|HDFN|6
954|HSC-I|HDFN|6
...

# 終了
sqlite> .q

以上で rawdata の準備は終了です。

較正用データの作成

フラットデータの作成

ドームフラットからフラットデータを作成します。これは全てのバンドで実行します。

reduceFlat.py /data/users/<User Name>/HSC \
        --rerun=calib_tutorial \
        --id field=DOMEFLAT filter=HSC-Y \ # これは Y バンドの場合です。filter ごとに計 3 回走らせます。
        --detrendId calibVersion=calib_tutorial \
        --config isr.doBias=False \
                isr.doDark=False \
        --batch-type=smp --cores=6

# rerun とは HSC pipeline で使われる解析プロセスの概念です。ある共通の解析パラメータを用いて画像を生成するまでの一つの解析プロセスを rerun と呼んでいます。
# --id は入力データを検索する条件を入れます。field と filter を指定していますが、visit 番号を指定することもできます (visit 212 から 220 まで二つ飛び;visit=212..220:2)。
# ただし、 field などで指定する場合は、 FITS ヘッダが間違っていたときに、どのデータが入力されたのかわかりにくくなるのでご注意下さい。
# 今回はバイアス、ダークは使用しないので、-- config isr.doBias=False、isr.doDark=False というオプションをつけています。
# --batch-type=smp --cores=6 は「ローカルコンピュータ上で処理を実行、6 スレッドで」という意味です。
# このあたりのパラメータについては後ほど解説します。

上のYバンドの例では、 /data/users/<User Name>/[dateObs]/[filter]/[calibVersion](今回は calibVersion = rerun)の下にフラットデータ (FLAT-[ccd].fits) ができます。ds9 で好きな画像ファイルを選んで、できた画像を確認してみましょう。 続いて、filter=HSC-Y を HSC-I, HSC-Z として同様に走らせましょう。

そして、出来上がったフラットデータをレジストリに登録します。これを忘れると、以降の解析がすべて失敗します。

# レジストリに登録
genCalibRegistry.py \
        --root=/data/users/<User Name>/HSC/CALIB \
        --camera=HSC \
        --create

これでフラットデータの作成は終了です。

フリンジデータの作成

フリンジデータは天体画像のスタックから生成されます。今回のデータのうち、Yバンドについてはフリンジデータから画像の干渉縞を除去する必要があります。 Pipeline ではフリンジデータ作成は Yバンドと NB0921 のみ対応しています。

# フリンジデータ(Y-band のみ)
reduceFringe.py /data/users/<User Name>/HSC \
        --rerun=calib_tutorial \
        --id field=HDFN filter=HSC-Y \
        --detrendId calibVersion=calib_tutorial \
        --config isr.doBias=False \
                isr.doDark=False \
        --batch-type=smp --cores=6

# 天体画像を使うため、field=HDFN を指定。

完了すると、 /data/users/<User Name>/HSC/CALIB/FRINGE/[dateObs]/[filter]/[calibVersion] の下にフリンジデータ (FRINGE-[ccd].fits) ができます。フリンジデータができていることを確認したら、フラットデータと同様にレジストリに登録しましょう。コマンドはフラットデータの時と同じです。どんな画像ができたのか気になる人は ds9 で確認してみましょう。

# レジストリに登録
genCalibRegistry.py \
        --root=/data/users/<User Name>/HSC/CALIB \
        --camera=HSC \
        --create

以上でフリンジデータの作成は終了です。

CCD ごとの解析

ここでは、天体データの 1 次処理後に、次の解析プロセスである天体データの位置合わせとフラックススケール決めのためにラフな等級原点と座標決めも行われます。そこでまずアストロメトリ用のカタログファイルを設定します。

# ログインの度に必ず行ってください。
setup astrometry_net_data sdss-dr9-fink-v5b

次に CCD 毎の解析プログラムを走らせます。まず、天体データの一次処理 (フリンジ除去、フラット割り) が行われます。そして、sky を引き、天体検出を行います。検出された天体はカタログファイルとマッチングされ、等級原点と座標決めが行われます (reduceFrames.py ではこの等級原点と座標決めが2回走ります)。最終的に一次処理済み画像と sky 画像などの画像データと、カタログデータが生成されます。

reduceFrames.py /data/users/<User Name>/HSC \
        --rerun=tutorial \
        --id field=HDFN \
        --config processCcd.isr.doBias=False \
                processCcd.isr.doDark=False \
                doSolveTansip=False \
        --batch-type=smp --cores=6

# --rerun=tutorial を指定しています。較正用データを作るための rerun をひとまず終了させて、データ処理のための新しい rerun を始める、ということです。
# --config doSolveTansip=False : これは視野全体の WCS を決定するプロセスで、デフォルトは True。今回は CCD 6 枚でしか座標決めを行っておらず、うまく解けない可能性があるのでスキップします。
# 以降で mosaic を行う場合はスキップしても大きな影響はありませんが、通常のすべての CCD を解析する場合は、精度を上げるためにもスキップしないほうが良いです。
reduceFrames を走らせている間に、コマンドのパラメータについて説明します。

こちらのページ -> パラメータの説明

reduceFrames が終了したら、まず画像データを確認していきましょう。
画像データの出力ディレクトリは /HSC/rerun/[rerun]/[pointing]/[filter]/corr、サムネイル画像の出力ディレクトリは /hsc/rerun/[rerun]/[pointing]/[filter]/thumbs です。ここで、[rerun] は –rerun=tutorial で指定した rerun 名、[pointing] は観測時のポインティングで、Y バンドでは 00815、i, z バンドでは 00816 になります。
つまり、Y バンドでは画像データは /data/users/<User Name>/HSC/rerun/tutorial/00815/HSC-Y/corr の下に生成されることになります。
  • 1 次処理済天体画像(sky も引き済)   :CORR-[visit]-[ccd].fits
  • sky 引きで用いられた sky 画像      :BKGD-[visit]-[ccd].fits
  • overscan 引きの結果サムネイル     :oss-[visit]-[ccd].png
  • Flat データで補正した天体画像サムネイル :flattened-[visit]-[ccd].png

次に、ラフな等級原点と座標決めのために用いられたカタログデータを紹介します。生成された全てのカタログデータは /data/users/<User Name>/HSC/rerun/[rerun]/[pointing]/[filter]/output の下にあります。

  • 1 度目の等級原点、座標決めのために天体画像から検出された比較的明るい天体のカタログ : ICSRC-[visit]-[ccd].fits
  • 1 度目の等級原点、座標決めに用いられた天体のマッチリスト (検出された天体のうち、アストロメトリ用カタログの中に一致するものがあった天体) : MATCH-[visit]-[ccd].fits
  • MATCH-[visit]-[ccd].fits をカラムに展開したもの : ML-[visit]-[ccd].fits
  • 2 度目の等級原点、座標決めのために天体画像から検出された天体のカタログ : SRC-[visit]-[ccd].fits
  • 2 度目の等級原点、座標決めに用いられた天体のマッチリスト : SRCMATCH-[visit]-[ccd].fits
  • SRCMATCH-[visit]-[ccd].fits をカラムに展開したもの : SRCML-[visit]-[ccd].fits

このリストの中にある MATCH-ML- の中の情報は基本的には同じです。しかし、MATCH- では HSC pipeline の中でのみ参照できる ID 番号が割り振られているのに対し、ML- は HSC pipeline 以外でも参照できる仕様になっています。 これらカタログデータは .fits の形式ですが画像データではありませんので、ds9 では中身を見ることはできません。 中身を確認するには。。。

—> こちら

以上で CCD ごとの解析は終了です。

位置合せとフラックススケール決め (モザイク)

tract 定義

全ての天体データを使い、 tract を定義します。天球面をぶつ切りにした領域を tract と呼び、一つの tract は一つの平面に投影されます。 1 視野くらいのデータの場合は、その画像を全て含むような tract を定義しなおします。

makeDiscreteSkyMap.py /data/users/<User Name>/HSC \
        --rerun=tutorial \
        --id field=HDFN

tract 0 has corners (189.998, 61.807), (188.417, 61.807), (188.398, 62.553), (190.018, 62.553) (RA, Dec deg) and 4 x 4 patches となれば成功です。ここで、patch とは tract をさらに細かく分けたものです。デフォルトの設定では 1 tract 中に 100 patch あり、 1 patch が 4200 × 4200 ピクセル四方の正方形で定義されています。これで tract 0 番が定義されました。

モザイク

モザイクでは、各 [visit, ccd] の等級原点と座標決めの情報に複数の visit の情報を加え、より精密に各 [visit, ccd] の位置合わせとフラックススケールを決定します。

mosaic.py /data/users/<User Name>/HSC \
        --rerun=tutorial \
        --id field=HDFN filter=HSC-Y tract=0 # filter を書き換えて計 3 回走らせます。
モザイクではフラックススケールファイルと WCS 座標ファイルが生成され、画像については特に生成されません。i、zバンドについても filter=HSC-I, HSC-Z として実行します。
では、生成されたファイルを確認していきましょう。/data/users/<User Name>/HSC/rerun/[rerun]/[pointing]/[filter]/corr/[tract] の下に fcr-[visit]-[ccd].fitswcs-[visit]-[ccd].fits というファイルができているはずです。これらがそれぞれフラックススケールファイルと WCS 座標ファイルです。[tract] は tract 番号で、今回は 0000 になります。

天体データの重ね合わせ (スタック)

モザイクによって生成されたデータを用いて、天球面座標を平面座標に投影したデータを生成し(warp)、観測された全ショットの積分を行います(coadd)。

stack.py /data/users/<User Name>/HSC \
        --rerun=tutorial \
        --id filter=HSC-Y tract=0 \ # filter はここと
        --selectId field=HDFN filter=HSC-Y \ # ここにあるので注意!!
        --batch-type=smp --cores=6
filter を HSC-I, HSC-Z として同様に走らせて下さい。上記のコマンドでは、filter 情報は二箇所に入りますので両方書き換えて下さい。
Coadd 画像は /data/users/<User Name>/HSC/rerun/[rerun]/deepCoadd/[filter]/[tract]/[patch] の下の calexp-[filter]-[tract]-[patch].fits です。また、warp 画像 (warp-[filter]-[tract]-[patch]-[visit].fits) も同じフォルダに生成されます。
さらに、/data/users/<User Name>/HSC/rerun/[rerun]/deepCoadd-results/[filter]/[tract]/[patch] には検出された天体のカタログ (det-[filter]-[tract]-[patch].fits) と sky 引きのパターン (bkgd-[filter]-[tract]-[patch].fits) ができます。この det- カタログが次のステップのマルチバンド解析で使われることになります。
スタック後の画像を一度に表示させたい場合は、Iバンドの場合であれば、
cd /data/users/<User Name>/HSC/rerun/tutorial/deepCoadd/HSC-I/0
ds9 -mosaic wcs */calexp-*.fits

で確認することが可能です (ただし、全 104 CCD でこれをやるとかなりのメモリを消費しますのでご注意ください)。

今回は中心の 6 CCD のみを使用しましたが、全ての CCD を使用した場合はこのようになります

–> こちら

以上でスタックが完了です。

マルチバンド解析

スタックによって生成されたそれぞれのフィルタの天体カタログをマージして、マルチバンドで測光しなおし、新たにカタログを生成します。

multiBand.py /data/users/<User Name>/HSC \
        --rerun=tutorial \
        --id tract=0 filter=HSC-Z^HSC-I^HSC-Y \ # 使うフィルタを全選択。
        --batch-type=smp --cores=6
# --id ではマルチバンド解析で使用する全てのデータを指定します。
まず最初に各 filter で検出された天体カタログを集めて天体の位置情報のマージカタログを作成します(mergeDet)。この時、デブレンド(2 つの天体が近接していて重なって見える時にその 2 つの天体を分離するプロセス)も行われます。 mergeDet カタログをもとに各 filter の stack 画像に戻って天体の測光を行います(meas)。 その後、各 filter の meas カタログで得られた天体の位置情報を再度集め直して、全ての filter で共通の座標情報を持つ天体カタログを作成します(ref)。 この ref カタログの位置情報をもとに天体の測光を filter 毎に改めて行い、forced_src- カタログが生成されます。
生成されるファイルを以下にまとめます。

/data/users/<User Name>/HSC/rerun/[rerun]/deepCoadd-results/merged/[tract]/[patch] ;

  • 天体の位置情報のマージカタログ : mergeDet-[tract]-[patch].fits
  • meas カタログをもとに生成される天体の位置情報カタログ : ref-[tract]-[patch].fits

/data/users/<User Name>/HSC/rerun/[rerun]/deepCoadd-results/[filter]/[tract]/[patch] ;

  • mergeDet をもとに測光した天体カタログ : meas-[filter]-[tract]-[patch].fit
  • 最終天体カタログ : forced_src-[filter]-[tract]-[patch].fits

スタック画像をレファレンスにした各 CCD の測光

ここでは、スタック後の天体カタログ (複数バンドをマージしたもの) の位置情報を使い、 各 [visit, ccd] データで天体を検出、測光し、カタログを生成し直すという作業が行われます。

forcedCcd.py /data/users/<User Name>/HSC \
        --rerun=tutorial \
        --id field=HDFN tract=0 \
        --batch-type=smp --cores=6

これを実行すると /data/users/<User Name>/HSC/rerun/[rerun]/[pointing]/[filter]/tract[tract] の下に FORCEDSRC-[visit]-[ccd].fits というカタログができます。

各 CCD 画像, カタログを較正しなおす

最後にモザイクによって決められた等級原点と WCS を、各 CCD 画像、カタログに反映させましょう。 WCS とゼロ点の較正された画像/カタログが不要な場合は行う必要はありませんが、講習会なので一通りやっていきます。

# 画像を較正しなおす
calibrateExposure.py /data/users/<User Name>/HSC \
        --rerun=tutorial \
        --id field=HDFN tract=0 \
        -j 6
# -j : スレッド数を指定。-j 6 とすると、6 スレッドで、という意味になる。


# カタログを較正しなおす
calibrateCatalog.py /data/users/<User Name>/HSC \
        --rerun=tutorial \
        --id field=HDFN tract=0 \
        --config doApplyCalib=False \
        -j 6
# --config doApplyCalib=False : カタログのフラックスを等級に較正しない場合は False にする
calibrateExposure.py が終了すると、/data/users/<User Name>/HSC/rerun/[rerun]/[pointing]/[filter]/corr/[tract] 以下に CALEXP-[visit]-[ccd].fits というファイルが生成されます。
また、calibrateCatalog.py で較正されたカタログは、/data/users/<User Name>/HSC/rerun/[rerun]/[pointing]/[filter]/output/[tract] 以下の CALSRC-[visit]-[ccd].fits です。

以上で全ての解析が終了です。

疑似カラー画像作成

それでは、解析で出来上がったデータを元に、疑似カラー画像を作成してみましょう。

  1. 必要なスクリプトをコピーしてきます。

    cd /data/users/<User Name>/
    wget http://hsc.mtk.nao.ac.jp/HSC_Training_download/script/stitchPatches.py
    # スクリプトの中身は各自で確認下さい。
    
    # 日本語が文字化けしている場合は、
    # メニューから「Terminal/Set Character Encoding/Unicord(UTF-8)」を選び、export LC_ALL=ja_JP.UTF-8 をしてください。
    
  2. パッチに分かれていた画像を1つの大きな画像にします。

    python stitchPatches.py -o HSC-Y.fits -x 5000 11000 -y 5000 11000 HSC/rerun/tutorial/deepCoadd/HSC-Y/0/*/calexp-*.fits
    python stitchPatches.py -o HSC-Z.fits -x 5000 11000 -y 5000 11000 HSC/rerun/tutorial/deepCoadd/HSC-Z/0/*/calexp-*.fits
    python stitchPatches.py -o HSC-I.fits -x 5000 11000 -y 5000 11000 HSC/rerun/tutorial/deepCoadd/HSC-I/0/*/calexp-*.fits
    
    # または、
    
    for f in Y Z I ; do python stitchPatches.py -o HSC-$f.fits -x 5000 11000 -y 5000 11000 HSC/rerun/tutorial/deepCoadd/HSC-$f/0/*/calexp-*.fits ; done
    
  3. ds9のカラーモードで表示します。

    ds9 -rgb -red HSC-Y.fits -green HSC-Z.fits -blue HSC-I.fits
    
  4. 色を調整します。

  • 各フィルタごとに次の操作をして色を調整します
    1. メニューから Frame/RGB... を選ぶ。
      • 「R」「G」「B」からチャンネルを選ぶ。
    2. メニューから Scale/Scale Parameters... を選ぶ。
      • スケールの最小値、最大値を調整する。
    3. メニューから Color/Colormap Parameters... を選ぶ。
      • Contrast, Biasを調整する。
  • (おすすめの調整の仕方)
    • フィルタごとに
      1. メニューから Scale/Log を選ぶ。

      2. メニューから Scale/Scale Parameters を選ぶ。
        • 最小値を0、最大値を100くらいにする。
  1. png形式で保存します。
    1. メニューから File/Export/png を選ぶ。
    2. HDFN.png の名前で保存する。
  2. 保存した画像を確認します。

    eog HDFN.png
    
    # または
    
    firefox HDFN.png
    

色等級図作成

最後に、カタログ情報から色等級図を描いてみましょう。
ここでは、butlerと呼ばれている、様々なタイプ(画像、カタログ)のデータを検索したりロードしたりするためのPipeline toolを使って色等級図の作成を行います。butlerの詳細について知りたい方は こちらのページ を参照してください。
  1. 必要なスクリプトをコピーします。

    cd /data/users/<User Name>
    wget http://hsc.mtk.nao.ac.jp/HSC_Training_download/script/catalog_all_{I,Y,Z}.py
    wget http://hsc.mtk.nao.ac.jp/HSC_Training_download/script/catalog.sh
    wget http://hsc.mtk.nao.ac.jp/HSC_Training_download/script/fig_izy.py
    
  2. 各バンドで作成したカタログから、butlerを用いて天体のIDとPSF fluxを取得し、一つのファイルにまとめます。

    python catalog_all_I.py /data/users/<User Name>/HSC/rerun/tutorial 0 HSC-I > HSC-I.txt
    python catalog_all_Z.py /data/users/<User Name>/HSC/rerun/tutorial 0 HSC-Z > HSC-Z.txt
    python catalog_all_Y.py /data/users/<User Name>/HSC/rerun/tutorial 0 HSC-Y > HSC-Y.txt
    sh catalog.sh
    
  3. 一つにまとめたファイルを使って、色等級図を表示します。

    python fig_izy.py
    
  4. 表示された図を保存します。

    1. save the figureのボタンを押す。
    2. 保存する名前を決めてsaveボタンを押す。
  5. 保存した図を確認します。

    eog hoge.png
    
    # または
    
    firefox hoge.png
    

スクリプトの中身は各自で確認してください。

以上で実習の内容はすべて終了です。お疲れ様でした。