Astro : 解析によく使う python モジュール/関数

• 天文データの解析でよく使う 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

...(加筆予定)...

データについて