- 天文データの解析でよく使う python のパッケージ/モジュールを紹介
- 興味の都合で、3次元 (空間:x, 空間:y, 視線速度:v) データの取り扱いに偏っていると思う
astropy
- web : http://www.astropy.org/
- document : http://docs.astropy.org/en/stable/
- github : https://github.com/astropy/astropy
天文データ解析パッケージ。元々は個別のパッケージとして開発されていたが、2012年くらいに統合された。データの取り扱いから単位系、座標系の操作まで天文解析固有の処理は概ねカバーされている印象 (パッケージが大きすぎてとても全部は知らない)。
astropy でよく使うモジュール
astropy.io.fits
FITS ファイルを読み書きするモジュール。pyfits から完全移行された (今後、pyfits の使用は不推奨)。FITS ファイルは天文データをなんでもかんでも保存できることを目的に作られた、かなり自由度の高いファイルフォーマット。通常は2次元もしくは3次元の画像データが座標情報付きで保存されるが、任意のデータ構造をテーブルとして保存することも可能。FITSファイルについての詳細は「FITSユーザーズガイド」が詳しい。
ipython
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
hdu = fits[0]
header = hdu.header
data = hdu.data
# ヘッダ情報を確認
print(repr(header))
naxis = header['NAXIS'] # 辞書のように扱える
naxis1 = header.get('NAXIS1')
# データ構造を確認
type(data)
>>> numpy.ndarray
data.ndim
>>> 2
data.shape
>>> (480, 720)
astropy.wcs
WCS (World Coordinate System) とは、FITS に記録されたデータが、天球上の座標にどのように対応するかを記述する仕様。詳細は、http://fits.gsfc.nasa.gov/fits_wcs.html。データの基準インデックス(CRPIX)、基準インデックスの物理座標値(CRVAL)、データグリッドの物理座標対応値(CDELT)、天球面への投影方法(CTYPE)などの情報をもとにして、インデックス値を物理座標値に変換する。
ipython
import astropy.wcs
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
wcs = astropy.wcs.WCS(fits[0].header)
print(wcs)
>>> ... # 情報が出力される
wcs.wcs_world2pix([[210, -19]], 0) # (l=210, b=-19) の index を取得
>>> array([[ 359.00009815, 120.00221429]])
wcs.wcs_pix2world([[150, 60]], 0) # (x=150, y=60) の物理座標を取得
>>> array([[ 213.48334194, -20.0000389 ]])
astropy.coordinates
お空の座標系 (赤経赤緯座標、天の川銀河座標など) を変換したり計算できるモジュール。
ipython
import numpy
import astropy.coordinates
import astropy.wcs
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
header = fits[0].header
data = fits[0].data
maxind = numpy.unravel_index(numpy.nanargmax(data), data.shape)
wcs = astropy.wcs.WCS(fits[0].header)
maxworld = wcs.wcs_pix2world([maxind[::-1]], 0)
# wcs_pix2world の要求するインデックスの軸の順は x, y(, z)
# fits[0].data の軸の順は、(z,) y, x
# この 2 つの向きが逆なので、maxind[::-1] で反転させている
galx, galy = maxworld[0]
# >>> (208.99999963578261, -19.383371004831144)
# 銀河座標で作成
coord = astropy.coordinates.SkyCoord(galx, galy, frame='galactic', unit='deg')
# >>> <SkyCoord (Galactic): (l, b) in deg
# (208.99999964, -19.383371)>
# 星座を確認
coord.get_constellation()
>>> u'Orion'
# 赤経赤緯 (epoch=J2000) で出力
coord.fk5
>>> <SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
(83.81461251, -5.38035224)>
# 赤経赤緯をフォーマットして出力
coord.fk5.to_string('hmsdms')
>>> u'05h35m15.507s -05d22m49.268s'
astropy.units
天文関連で使う単位を取り扱うモジュール。
ipython
import astropy.units
# 1: 5km/s で動く分子雲コアが 5pc 移動するのにかかる時間を計算してみる
v = 5 * astropy.units.kilometer / astropy.units.second
# >>> <Quantity 5.0 km / s>
d = 5 * astropy.units.pc
# >>> <Quantity 5.0 pc>
(d/v).to('megayear')
>>> <Quantity 0.9777922216731284 Myr>
# 2: 10^-6 Mo/yr の質量降着率を、kg/sec にしてみる
(1e-6 * astropy.units.solMass / astropy.units.year).to('kg/second')
>>> <Quantity 6.303077547088497e+16 kg / s>
APLpy
- web : https://aplpy.github.io/
- document : http://aplpy.readthedocs.org/en/stable/
- github : https://github.com/aplpy/aplpy
FITS データの可視化パッケージ。論文出版に使えるクオリティの画像を出力できる。
使用例
ipython
import matplotlib.pyplot
import astropy.io.fits
import aplpy
# FITS ファイルを web からダウンロードする場合
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
fig = aplpy.FITSFigure(fits)
# FITS ファイルが local に保存されている場合
# fig = aplpy.FITSFigure('Orion.CO1221.Osaka.beam204.mom0.fits.gz')
fig.show_colorscale()
fig.add_colorbar()
# fig.save('Orion.CO1221.Osaka.beam204.mom0.png') # 画像を保存する場合
matplotlib.pyplot.show()
numpy
- web : http://www.numpy.org/
数値計算をする上での基盤となるデータ構造、操作を提供するモジュール。
numpy でよく使う関数
numpy.sum, .nansum
numpy.ndarray もしくは list の、合計値を計算する。
axis 引数を指定すれば、積算する軸を指定できる。
ipython
import numpy
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
# 未観測の領域 (voxel) には nan が入っていることが多いので、nansum を使う
numpy.nansum(fits[0].data)
>>> 1630958.2 # 単位は K km/s pix^2
# 3D データ (cube) から 2D データ (積分強度) を出すときにも使う
# !!!! ダウンロード容量注意 (1.8GB) !!!!
# fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.cube.fits.gz')
fits = astropy.io.fits.open('Orion.CO1221.Osaka.beam204.cube.fits.gz')
fits[0].data.shape
>>> (2521, 480, 720) # z 軸 (速度), y 軸 (横), x 軸 (縦)
mom0 = numpy.nansum(fits[0].data, axis=0) # 0 軸を合計する
mom0.shape
>>> (480, 720)
numpy.min, .max, .nanmin, .nanmax
numpy.ndarray もしくは list の、最小、最大を取得する。
ipython
import numpy
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
numpy.nanmin(fits[0].data)
>>> -0.61183316 # K km/s
numpy.nanmax(fits[0].data)
>>> 431.83881 # K km/s
numpy.argmin, .argmax, .nanargmin, .nanargmax
numpy.ndarray もしくは list の、最小、最大値を保有する index を取得する。
ipython
import numpy
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
data = fits[0].data
maxind = numpy.nanargmax(data)
# >>> 70259
maxcoord = numpy.unravel_index(maxind, data.shape)
# >>> (97, 419)
data[maxcoord] == numpy.nanmax(data)
>>> True
numpy.meshgrid
1次元の array を二つ使って、2次元の array を作る。
座標情報のマップが欲しいときに使う。
ipython
import numpy
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
def get_axis(header, axis=1):
crval = header.get('CRVAL%1d'%(axis))
crpix = header.get('CRPIX%1d'%(axis))
cdelt = header.get('CDELT%1d'%(axis))
num = header.get('NAXIS%1d'%(axis))
return crval + (numpy.arange(num) - crpix + 1) * cdelt
x1d = get_axis(fits[0].header, 1)
y1d = get_axis(fits[0].header, 2)
X2d, Y2d = numpy.meshgrid(x1d, y1d)
X2d.shape
>>> (480, 720)
# 3d data の場合 (meshgrid が使えないので以下のようにする)
# !!!! ダウンロード容量注意 (1.8GB) !!!!
# fits_cube = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.cube.fits.gz')
fits_cube = astropy.io.fits.open('Orion.CO1221.Osaka.beam204.cube.fits.gz')
x1d_ = get_axis(fits_cube[0].header, 1)
y1d_ = get_axis(fits_cube[0].header, 2)
v1d_ = get_axis(fits_cube[0].header, 3)
X3d = x1d_[None, None, :]*1 + y1d_[None, :, None]*0 + v1d_[:, None, None]*0
Y3d = x1d_[None, None, :]*0 + y1d_[None, :, None]*1 + v1d_[:, None, None]*0
V3d = x1d_[None, None, :]*0 + y1d_[None, :, None]*0 + v1d_[:, None, None]*1
X3d.shape
>>> (2521, 480, 720)
# None の代わりに numpy.newaxis も使える (そちらのほうが推奨かも？)
numpy.polyfit, .polyval
多項式でのフィッティングと、フィッティングパラメータから曲線を復元する関数。
scipy
scipy でよく使う関数
scipy.signal.convolve
scipy.ndimage.map_coordinates
scipy.interpolate.interp1d
matplotlib
...(加筆予定)...