はじめに
Pythonとmatplotlibを使って、チャンドラ（Chandra）X線観測衛星のFITS画像に光軸（観測時のRA / DEC / ROLL）を可視化したい——
そう思って調べても、意外とシンプルな方法がまとまっていなかったりします。
ヘッダーを見たり、DS9に手入力したりと、地味に手間がかかりますよね。
でも、この光軸の位置って実はとても重要で、解析の精度に直結します。たとえば複数の観測をmergeするときに位置がずれていると、PSF（Point Spread Function）の形が崩れ、構造がぼやけてしまうこともあります。実際に可視化してみると、PSFの広がりやズレの影響が見えてくることもあります。
本記事では、Python・matplotlibを用いて、超新星残骸 Cassiopeia A（以下、Cas A）のFITS画像を例に、光軸の位置を画像上にかっこよく可視化する方法をご紹介します。
RA / DECはhms / dms形式でタイトルに表示し、DS9風のカラーマップ（ds9_b）も再現可能。論文やスライド資料にもそのまま使える、シンプルかつ見栄えのいい可視化手法です。
実装コードはGoogle Colab こちら からも閲覧できます。
コードデモ
ここでは、Cas Aを例に紹介しますが、他のChandra衛星の観測画像でも活用可能です。
デモに使用した画像はこちらから入手できます。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colormaps
from matplotlib.colors import LogNorm, LinearSegmentedColormap
from scipy.ndimage import gaussian_filter
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import astropy.units as u
# --- DS9 B colormap 登録（天体画像用） ---
ds9_b_cdict = {
'red': [(0.0, 0.0, 0.0), (.25, 0.0, 0.0), (.5, 1.0, 1.0), (1.0, 1.0, 1.0)],
'green': [(0.0, 0.0, 0.0), (.5, 0.0, 0.0), (.75, 1.0, 1.0), (1.0, 1.0, 1.0)],
'blue': [(0.0, 0.0, 0.0), (.25, 1.0, 1.0), (.5, 0.0, 0.0), (.75, 0.0, 0.0), (1.0, 1.0, 1.0)]
}
ds9_b_cmap = LinearSegmentedColormap('ds9_b', ds9_b_cdict)
colormaps.register(name='ds9_b', cmap=ds9_b_cmap)
# --- FITSファイル読み込み ---
filename = 'casA_4636_broad_flux.img'
with fits.open(filename) as hdul:
data = hdul[0].data
header = hdul[0].header
# --- WCSおよびpointing情報取得 ---
wcs = WCS(header)
ra_pnt = header['RA_PNT']
dec_pnt = header['DEC_PNT']
roll_pnt = header['ROLL_PNT']
sky = SkyCoord(ra=ra_pnt * u.deg, dec=dec_pnt * u.deg)
ra_str = sky.ra.to_string(unit=u.hour, sep=':', precision=2, pad=True)
dec_str = sky.dec.to_string(unit=u.deg, sep=':', precision=2, alwayssign=True, pad=True)
# --- 有効なピクセル値からvmin/vmaxを自動設定（2–100%分位） ---
valid_data = data[np.isfinite(data) & (data > 0)] # LogNormなので0以下は除外
vmin, vmax = np.percentile(valid_data, [2, 100])
vmin = max(vmin, 1e-10) # vminが0に近すぎる場合の保険
# --- pointing位置をピクセル座標に変換 ---
x_pnt, y_pnt = wcs.world_to_pixel_values(ra_pnt, dec_pnt)
# --- 表示用データを作成（0以下やNaNはvminに） ---
view_data = np.where((data > 0) & np.isfinite(data), data, vmin)
# --- 画像表示（logスケール + DS9 Bカラーマップ） ---
plt.figure(figsize=(8, 5))
ax = plt.subplot(projection=wcs)
im = ax.imshow(view_data, cmap='ds9_b', norm=LogNorm(vmin=vmin, vmax=vmax), origin='lower')
ax.scatter(x_pnt, y_pnt, s=100, marker='x', color='lime', linewidths=2)
ax.set_xlabel('RA (J2000)')
ax.set_ylabel('Dec (J2000)')
ax.set_title(f'Pointing:\nRA={ra_pnt:.3f} ({ra_str}),\n Dec={dec_pnt:.3f} ({dec_str}),\n Roll={roll_pnt:.1f}')
plt.colorbar(im)
plt.tight_layout()
plt.show()
出力結果
コードのポイント
このコードは、Chandra衛星のFITS画像から観測時のpointing情報（RA / DEC / ROLL）を取得し、画像上にその位置を可視化するものです。
主な処理は以下のとおりです：
- FITSファイルから画像データとヘッダーを読み込み
- ヘッダーに記録された光軸情報（RA/DEC/ROLL）を取得
- WCSを使って、RA/DECを画像上のピクセル座標に変換
- 画像をLogスケールで描画し、光軸位置に×印を重ねて表示
- カラーマップはDS9カラーb（
ds9_b）を再現（詳しくはこちらの記事をご参照ください）
表示には、LogNormによる0除算を避けるため、ゼロやNaNを含む画素を最小表示値 vmin に置き換えた view_dataを使用しています。これにより、データの構造を保ちつつ、視認性の高い表示が可能になります。
おわりに
Chandraは、打ち上げから25年経った今も、最高の空間分解能で私たちに宇宙の姿を届けてくれる、素晴らしいX線衛星です。
それでも、光軸がわずかにずれるだけで若干ボケるし、視野内の位置によっても歪みの程度が違ってきます。
だからこそ、私たちはこうして1ピクセルの違いにも目を向け、見落とさないように工夫を重ねています。
実際、ありがたいことに、Chandraの“ボケ方”は先人たちの研究によってすでにかなり解明されており、PSFとして数値的に記述し、シミュレーションすることも可能です。
こうしたPSFをうまく取り入れることで、観測のズレやボケを補正し、より本来の構造に近づけることもできます。
私自身も、PSFの影響を考慮した上で空間構造を検討する試みとして、以下のような研究を行いました：
とはいえ、PSFの扱いも解析技術も、まだまだ発展途上です。
この拙いコードや視点が、どこかで誰かの気づきへとつながっていたなら、何よりです。
そして最後に。
Chandraの膨大な記憶の房に、ひっそりと埋もれていたひと粒を手に取るように——
忘れられかけたその果実が、またひとつの知として息を吹き返すのなら、
それほど嬉しいことはありません。