はじめに
データ解析や観測結果の処理において、領域をスライスすることは一般的な作業ですが、そのスライスに対応するようにWCS（World Coordinate System）座標を取得したい場合があります。
本記事では、PythonのライブラリであるAstropyを使用して、FITSファイルからデータをスライスし、新しいFITSファイルに保存する方法について説明します。
方法
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from datetime import datetime
# FITSファイルのパス
infile = 'infile.fits'
# FITSファイルを開く
hdul = fits.open(infile)
# データとWCS情報を取得
image = hdul[0].data
wcs = WCS(hdul[0].header)
# データをスライスするための範囲を定義
view = (slice(100, 151, 1), slice(10, 61, 1))
# データをスライス
sliced_image = image[view]
# WCSをスライスに合わせて更新
sliced_wcs = wcs[view]
# ヒストリー情報を更新
history = "Data sliced using NumPy. Original data shape: {}. Sliced shape: {}. Date: {}".format(image.shape, sliced_image.shape, datetime.now().isoformat())
# 元のヘッダー情報をコピー
original_header = hdul[0].header.copy()
# 新しいFITSファイルにデータとWCS情報を書き込む
outfile = 'outfile.fits'
sliced_hdul = fits.PrimaryHDU(data=sliced_image, header=original_header)
sliced_hdul.header.extend(sliced_wcs.to_header(), update=True) # WCS情報をヘッダーに追加
sliced_hdul.header.add_history(history) # ヒストリー情報を追加
sliced_hdul.writeto(outfile, overwrite=True)
コードの説明
- 必要なライブラリをインポート
- FITSファイルを開いて、データとWCS情報を取得
- データをスライスするための範囲を定義し、その範囲でデータをスライスします。WCS情報もスライスに合わせて更新
（この例では、スライスの範囲を[100:151, 10:61:1]と指定していますが、sliced_image = image[100:151:1, 10:61:1]やsliced_wcs = wcs[100:151:1, 10:61:1]と書いても同じ結果が得られます。ただ、同じスライス内容を2度書くことは誤解を招く恐れや可読性を下げるため、一つのスライスオブジェクトを使用して書いています。スライスオブジェクトについての説明は、Pythonのスライスはコロン":"だけじゃない？slice()の活用法をご覧ください。） - ヒストリー情報を更新し、新しいFITSファイルにデータとWCS情報を書き込む
- FITSファイルからデータをスライスし、新しいFITSファイルに保存
補足: wcs.slice()を使った書き方
wcs.slice()メソッドを使用してWCS情報を更新する方法もあります。
# WCSをスライスに合わせて更新
# sliced_wcs = wcs[view] # numpy_order=Trueの場合、通常のスライスと同じ意味
sliced_wcs = wcs.slice(view, numpy_order=True)
numpy_order=Falseに設定すると、WCSはFITS順序でスライスされます。つまり、最初のスライスは、最後のnumpyインデックスに適用されますが、最初のWCS軸に適用されます。
詳しくはAstropyの公式ドキュメントをご覧ください。
まとめ
本記事では、Astropyを使用してFITSファイルからデータとWCS座標をスライスし、新しいFITSファイルに保存する方法を説明しました。スライス操作を使用することで、任意の領域のデータを抽出し、解析することが容易になります。また、WCS情報も適切に更新することで、空間座標系の情報を保持したままデータを操作できます。