Pipeline tools¶
HSC pipeline には python をベースにした Pipeline tool が同梱されています。この Pipeline tool を使うと画像データやカタログデータを確認することができます。データを検索したりロードする時には butler, dataRef が、画像を処理するには Exposures, MaskedImages, Images、カタログ内の天体の座標、フラックス、サイズ等の情報を一覧にして取り扱うには SourceCatalogs という tool が使えます。このページでは Pipeline tool の使い方を簡単に紹介します。
butler¶
butler は様々なタイプ（画像、カタログ）のデータを検索したりロードしたりするための Pipeline tool です。HSC pipeline で生成されたデータをbutler で読み込む際には get() と dataId を使います。get() では読み込みたいデータの種類 （ target ）と、データの [visit, ccd]、coadd データの場合は [tract, patch] を指定します。
# butler を呼び出す python モジュール
import lsst.daf.persistence as dafPersist
# 呼び出したいデータのある rerun ディレクトリを指定
dataDir = "~/hsc/rerun/dith_16h_test"
# butler の呼び出し
butler = dafPersist.Butler(dataDir)
# 読み込みたいデータを指定し、butler に検索させる
#
# 今回は CORR-0902798-59.fits 画像
# 'calexp' は CORR-*.fits の意味
dataId = {'visit': 902798, 'ccd': 59}
exp = butler.get('calexp', dataId)
例では CORR-.fits を読み込むために target = ‘calexp’ を設定しています。 これ以外にも様々な target がありますので、以下に一覧を表示します。
| Target | データ種類 | データタイプ |
|---|---|---|
| bias | Bias データ | ExposureF |
| dark | Dark データ | ExposureF |
| flat | Flat データ | ExposureF |
| fringe | Fringe データ | ExposureF |
| postISRCCD | 一次処理済データ（デフォルトでは生成されない） | ExposureF |
| calexp | sky 引き済一次処理済データ | ExposureF |
| psf | 解析で使用される PSF | psf |
| src | 一次処理済データから作られる天体カタログ | SourceCatalog |
| wcs, frc | mosaic.py で用いられた天体カタログ | ExposureI |
| Target | データ種類 | データタイプ |
|---|---|---|
| deepCoadd_calexp | coadd データ | |
| deepCoadd_psf | coadd 画像の PSF | |
| deepCoadd_src | coadd 画像から作れたカタログファイル |
butler 読み込んだデータの種類は dataRef で調べることができます。そこで、先ほど butler で検索したデータの情報を調べてみましょう。
# butler を使い HSC pipeline のデータを検索できる python モジュール
import hsc.pipe.base.butler as hscButler
# dataRef でデータ情報を読み込む
dataRef = hscButler.getDataRef(butler, dataId)
ref_exp = dataRef.get('calexp')
ref_exp
# 出力：
# <lsst.afw.image.imageLib.ExposureF; proxy of <Swig Object of type 'boost::shared_ptr< lsst::afw::image::Exposure< float,lsst::afw::image::MaskPixel,lsst::afw::image::VariancePixel > > *' at 0x7f24376a5d50> >
# 同じ dataId を持つカタログファイルの情報を読み込む
ref_cat = dataRef.get('src')
# 出力：
# <lsst.afw.table.tableLib.SourceCatalog; proxy of <Swig Object of type 'lsst::afw::table::SortedCatalogT< lsst::afw::table::SourceRecord > *' at 0x7f24376a5f90> >
# 同じ dataId データの解析で使用された PSF の情報を読み込む
ref_psf = dataRef.get('psf')
# 出力：
# <lsst.afw.detection.detectionLib.Psf; proxy of <Swig Object of type 'boost::shared_ptr< lsst::afw::detection::Psf > *' at 0x7f24376a5e70> >
Exposures, MaskedImage, Images¶
次に、画像を表示する 3 つの Pipeline tool についてご紹介します。これら 3 つの tool のうち最もシンプルなのは Images です。Images では 2 次元画像を読み込みます。一方、 MaskedImage では CORR-.fits や [patch].fits に含まれる 3 つの画像データを読み込み （天体画像、マスク画像、variance 画像）、Exposures では MaskedImage が読み込む画像とそれに付随するヘッダー情報を読み込みます。
以下の例では、get() を使い、MaskedImage を呼び出し、sky 引き済一次処理済データを matplotlib で表示・保存しています。
# python モジュールの呼び出し
# python ベースの plotter を呼び出している
import matplotlib.pyplot as pyplot
import numpy
import argparse
# sky 引き済一次処理済データを読み込み
mimg = exp.getMaskedImage()
img = mimg.getImage()
# 画像データを array 形式に変換し、画像を test.png として保存
nimg = img.getArray()
pyplot.imshow(numpy.arcsinh(nimg), cmap='gray')
pyplot.gcf().savefig("test.png")
図1 test.png
MaskedImage により読み込んだ画像は、ds9 でも開くことができます。開き方は以下の通りです。
# ds9 や画像表示用のモジュールを呼び出す
import lsst.afw.display.ds9 as ds9
import lsst.afw.image as afwImage
# ds9 で画像を表示
ds9.mtv(mimg)
図2 ds9 で表示した sky 引き済一次処理済データ。左から、マスク画像、マスク画像を重ねた天体画像
図3 ds9 でマスク画像のパラメーターを調整する方法
図2 では sky 引き済一次処理済データを ds9 で表示した際の結果を表示しています。ds9 で画像を開くと、 上からマスク画像、天体画像、variance 画像の順で、全て 1 つの画像ウィンドウに重ねて表示されます。 図2左 のようにマスク画像を半透明にし天体画像を重ねて表示させたい場合は、ds9 のツールバーから [Analysis] > [Mask Parameters] を開き（図3）、Transparency （図2右） を小さくすればマスク画像が半透明になります。
次にマスク画像について少し詳しく説明します。マスクとして flag がつくピクセルの多くは、 宇宙線やサチッた星によるオーバーフローなどの問題があるピクセルです。それ以外にも、検出された天体の footprint の位置にも flag がつけられる場合があります（表3）。ds9 でマスク画像を表示すると、 こうしたマスクの種類毎に異なった色がつけられています。表3 に対応表を示しています。 例えば、マスク画像（図2左）を見ると、表3 にあるように、検出した天体は青、 サチっている星のオーバーフローは緑、バッドピクセル（だと思われるピクセル）は赤が着色されています。
| flag のラベル | 意味 | ds9 上の色 |
|---|---|---|
| BAD | バッドピクセル（HSC カメラ側の問題） | 赤 |
| CR | 宇宙線 | マゼンタ |
| CROSSTALK | ピクセルのクロストーク | |
| EDGE | 端にある CCD | 黄色 |
| INTERPOLATED | 周囲のピクセルから補間されたカウント値になっているピクセル | 緑 |
| INTRP | INTERPOLATED と同義 | |
| SATURATED | サチったピクセルからのオーバーフロー | 緑 |
| SAT | SATURATED と同義 | |
| SUSPECT | ほぼ サチっているピクセル（非線形性がうまく補正されていないため生じている） | 黄色 |
| UNMASKEDNAN | ISR 画像でピクセルの値が NaN になっているピクセル | |
| DETECTED | 検出された天体の footprint の一部 | 青 |
| DETECTED_NEGATIVE | 天体の footprint の一部で値が 負 になっているピクセル | シアン |
| CLIPPED | （coadd データのみ）coadd 処理においてクリップ (coadd に使われていない) されたピクセル | |
| NO_DATA | （coadd データのみ）coadd 処理において入力データがないピクセル |
SourceCatalog¶
HSC pipeline で生成された天体カタログは SourceCatalog と butler を使って検索・参照できます。 天体カタログには、各列に天体が、各行には天体の測定結果（フラックス、位置等）が格納されています。 各行に含まれる測定結果のうち、変数が float や int 型のものは sources.get(” ”) を使って抽出できます。また、全ての天体に対して同じ測定量を調べたい場合は、 ループ処理を使ってアクセスすることも可能です。カタログファイルの全測定情報については 一次処理済天体データのカタログファイル（src）中身 や coadd 画像による天体カタログ（deepCoadd_src）の中身 をご覧ください。
以下では、カタログ内にある天体の PSF フラックスと、そのフラックスを等級に変換して表示する場合の 2 例を示します。
# python モジュールを呼び出す
import numpy
""" PSF フラックス を調べる """
# butler で SourceCatalog を読み込む
sour = butler.get("src", dataId)
# カタログ内の天体数を取得
n = len(sources)
n
# > 1922
# PSF フラックスの値を取得
psfflux = sources.get("flux.psf")
# 各天体の PSF フラックスと 'extendedness' の値を SourceRecord から抽出
for i, src in enumerate(sour):
print i, psfflux[i], src.get("classification.extendedness")
#（結果: 左から、天体番号[0~1922], PSF フラックス, extendedness）
# 0 nan 1.0
# 1 nan 1.0
# 2 nan 1.0
# 3 nan 1.0
# :
# :
# 1918 600.487226674 1.0
# 1919 618.323853071 1.0
# 1920 578.070700843 1.0
# 1921 nan 1.0
""" PSF フラックス を等級に変換する """
# CCD の原点情報を取得する
metadata = butler.get("calexp_md", dataId)
zeropoint = 2.5 * numpy.log10(metadata.get("FLUXMAG0"))
zeropoint
# > 30.595966894420105
# PSF フラックスを等級に変換する
psfmag = zeropoint - 2.5 * numpy.log10(psfflux)
# 各天体の PSF フラックスから求めた等級の値を SourceRecord から抽出
for i, src in enumerate(sour):
print i, psfmag[i]
#（結果: 左から、天体番号[0~1922], PSF フラックス, extendedness）
# 0 nan
# 1 nan
# 2 nan
# 3 nan
# :
# :
# 1918 23.6497074602
# 1919 23.6179268929
# 1920 23.6910144995
# 1921 nan
このように得られた天体の測定結果を使って matplotlib で図を作ることが可能です。例えば、PSF フラックスから得られた等級の頻度分布を調べてみます。
# python モジュールを呼び出す
from matplotlib.pyplot import *
from matplotlib.figure import *
# PSF フラックスから得た等級の頻度分布を作成し、表示
hist(psfmag, bins=40, range=(10,30))
xlabel('magnitude(PSF)')
show()
図4 PSF フラックスから得た等級の頻度分布
次に、カタログ内の座標情報について調べる方法を紹介します。HSC pipeline で生成されたカタログでは、 天体の座標情報は coord という変数の中に格納されています。この座標情報は Angle と呼ばれる特殊なデータ形式をしており、その中には RA, DEC, 座標変換の情報が含まれています。coord には ICRS, FK5, Galactic, Ecliptic（ただし ICRS と FK5 は同じもの）の座標系が含まれています。 また、座標の単位としては degrees, radians, arcminutes, arcseconds という様々なフォーマットが利用可能です。
では、カタログの天体の様々な座標系での情報を調べてみることにします。
# カタログ内の全ての天体について座標情報を取得（ベースは ICRS 座標系）
for src in sour[0:n]:
icrs = src.get('coord')
galactic = icrs.toGalactic()
fk5 = icrs.toFk5()
ecliptic = icrs.toEcliptic()
# ICRS 系での RA, DEC [deg] を取得
ra, dec = icrs.getRa(), icrs.getDec()
# 他座標系での座標情報を取得
l, b = galactic.getL(), galactic.getB()
ra2, dec2 = fk5.getRa(), fk5.getDec()
lamb, beta = ecliptic.getLambda(), ecliptic.getBeta()
lon, lat = icrs.getLongitude(), icrs.getLatitude()
sid = src.getId()
# 座標情報をターミナルに表示
print "ID: ", sid
print " ICRS RA/Dec (deg)", ra.asDegrees(), dec.asDegrees()
print " FK5 RA/Dec (rad)", ra2.asRadians(), dec2.asRadians()
print " Galactic l/b (arcmin)", l.asArcminutes(), b.asArcminutes()
print " Ecliptic lamb/beta (arcsec)", lamb.asArcseconds(), beta.asArcseconds()
print " Generic Long/Lat (str)", lon, lat
""" print の結果
ID: 775497830381912065
ICRS RA/Dec (deg) 237.77835192 10.1021786534
FK5 RA/Dec (rad) 4.15001513096 0.176316279126
Galactic l/b (arcmin) 1189.90186102 2667.63454255
Ecliptic lamb/beta (arcsec) 838482.263154 106153.9408
Generic Long/Lat (str) 4.15002 rad 0.176316 rad
ID: 775497830381912066
ICRS RA/Dec (deg) 237.778311856 10.0840572239
FK5 RA/Dec (rad) 4.15001443172 0.176000000516
Galactic l/b (arcmin) 1188.5591858 2667.12102399
Ecliptic lamb/beta (arcsec) 838500.362538 106090.634831
Generic Long/Lat (str) 4.15001 rad 0.176 rad
:
:
ID: 775497830381913985
ICRS RA/Dec (deg) 237.585304192 10.1094545794
FK5 RA/Dec (rad) 4.14664581249 0.176443267991
Galactic l/b (arcmin) 1182.84956225 2677.87959672
Ecliptic lamb/beta (arcsec) 837712.866183 106012.221798
Generic Long/Lat (str) 4.14665 rad 0.176443 rad
ID: 775497830381913986
ICRS RA/Dec (deg) 237.583434907 10.1064926507
FK5 RA/Dec (rad) 4.14661318733 0.176391572584
Galactic l/b (arcmin) 1182.55610941 2677.89237965
Ecliptic lamb/beta (arcsec) 837708.488308 106000.261247
Generic Long/Lat (str) 4.14661 rad 0.176392 rad
"""