本記事では、天文データ解析でよく使う list や numpy array の操作について思いついたものを記述します。思いつき次第追記していきます。

途中、例として国立天文台の FUGIN プロジェクトで得られた野辺山45m電波望遠鏡の CO 輝線のアーカイブデータを用います。データは
http://jvo.nao.ac.jp/portal/nobeyama/
から
FGN_03100+0000_2x2_12CO_v1.00_cube.fits
をダウンロードしました (重いです)。

まず必要なものを import します。

import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt

list の生成

a = [0, 7, 2, 8] # 直打ち
print(a)
#[0, 7, 2, 8]

a = [] # 空の list を用意
for i in range(1, 6): # for文で一つずつ入れていく
    a.append(i)
print(a)
#[1, 2, 3, 4, 5]

a = []
for i in range(5): # range(0, 5) と等価
    a.append(i)
print(a)
#[1, 2, 3, 4, 5]

a = [i for i in range(5)] # リスト内包表記 (上と等価)
print(a)
#[1, 2, 3, 4, 5]

a = [0, 7, 2, 8] * 3 # 繰り返し
print(a)
#[0, 7, 2, 8, 0, 7, 2, 8, 0, 7, 2, 8]

a = [] * 3
print(a)
# []

a = [[]] * 3
print(a)
# [[], [], []]

a = [{}, {}] # list の要素は辞書などでも ok
print(a)
# [{}, {}]

a = np.array([0, 7, 2, 8])
a = list(a) # numpy array 等から list へ変換
print(a)
#[0, 7, 2, 8]

numpy array の生成

a = np.array([0, 7, 2, 8]) # list を numpy array に変換
print(a)
#[0 7 2 8]

a = np.array([0, 7, 2, 8], dtype="float32") # 要素の型を指定
print(a)
#[0. 7. 2. 8.]
a = a.astype("int")
print(a)
#[0 7 2 8]

a = np.arange(4) # np.array([i for i in range(4)]) と等価
print(a)
#[0 1 2 3]

a = np.linspace(0, 2, num=5, endpoint=True) # linear で生成
print(a)
#[0.  0.5 1.  1.5 2. ]

a = np.logspace(0, 2, num=5, endpoint=True) # log で生成
print(a)
#[  1.           3.16227766  10.          31.6227766  100.        ]

a, b = np.zeros(4), np.ones(4) # 0 や 1 の array を生成
print(a, b)
#[0. 0. 0. 0.] [1. 1. 1. 1.]
a = np.zeros((2, 2)) # 2D 
print(a)
#[[0. 0.]
# [0. 0.]]
b = np.ones(4)
a = np.zeros_like(b) # 同じ shape の zeros array を返す
print(a)
# [0. 0. 0. 0.]

list の操作

a = [0, 7, 2, 8]

len(a) # 長さを取得
max(a), min(a) # 最大値, 最小値 (NaN は無視される)
[i for i in a if i%2==0] # 条件に合うもののみにする (この場合偶数)
a[::2] # n個ごと (n-1個とばし) にとってくる (この場合1個とばし)
a[::-1] # 逆順にする

numpy array の操作

a = np.array([0, 7, 2, 8])

len(a) # (0軸目の) 長さを取得
a.shape # shape を取得
a.reshape(2, 2) # shape を変更
np.max(a), np.min(a) # 最大値, 最小値 (NaN が含まれていると NaN が返ってくる)
np.nanmax(a), np.nanmin(a) # 最大値, 最小値 (NaN は無視される)
np.mean(a), np.nanmean(a), np.median(a), np.nanmedian(a) # 平均値, 中央値
np.argmax(a), np.argmin(a) # 最大値, 最小値の index が返ってくる (複数ある場合は最初の index)
                           # (また、多次元 array の場合、1列にした時の index が返されるので以下のようにする)
a = np.array([[0, 7, 2, 8], [9, 3, 5, 1]])
print(np.unravel_index(np.argmax(a), a.shape))
#(1, 0)

b = np.arange(24).reshape(2, 3, 4)
print(b)
#[[[ 0  1  2  3]
#  [ 4  5  6  7]
#  [ 8  9 10 11]]
# [[12 13 14 15]
#  [16 17 18 19]
#  [20 21 22 23]]]
a = np.swapaxes(b, 0, 1) # 軸を入れ替える
print(a)
#[[[ 0  1  2  3]
#  [12 13 14 15]]
# [[ 4  5  6  7]
#  [16 17 18 19]]
# [[ 8  9 10 11]
#  [20 21 22 23]]]
a = b.T # 転置
print(a)
#[[[ 0 12]
#  [ 4 16]
#  [ 8 20]]
# [[ 1 13]
#  [ 5 17]
#  [ 9 21]]
# [[ 2 14]
#  [ 6 18]
#  [10 22]]
# [[ 3 15]
#  [ 7 19]
#  [11 23]]]

b = np.array([0, 7, 2, 8])
a = b * 3 # 要素に掛け算
print(a)
#[ 0 21  6 24]
a = b * b
print(a)
#[ 0 49  4 64]

fits 画像への適用例

基本的な map

import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt

hdu = fits.open("~/your/fits/dir/FGN_03100+0000_2x2_12CO_v1.00_cube.fits")[0]
data = hdu.data
print(data.shape)
#(462, 848, 848) # V, Y, X

channel map

ch_map = data[300] # 300チャンネル目をスライス
plt.imshow(ch_map, vmin=0, vmax=20, origin="lower left")
plt.show()

integrated intensity map

integ_map = np.nansum(data[280:320], axis=0)*hdu.header["CDELT3"]/1000.0 # 280-319チャンネルを積分
plt.imshow(integ_map, vmin=0, vmax=300, origin="lower left")
plt.show()

参考

peak intensity map

peak_map = np.nanmax(data, axis=0) # 全チャンネルでの最大値
plt.imshow(peak_map, vmin=0, vmax=30, origin="lower left")
plt.show()

position–velocity diagram

integ_map = np.nansum(data, axis=1)*np.abs(hdu.header["CDELT2"]) # Y軸方向に積分
plt.imshow(integ_map, vmin=0, vmax=8, origin="lower left")
plt.show()

参考

map への mask

cut = 1.5
ch_map = data[300] # 300チャンネル目をスライス
ch_map[ch_map<cut] = np.nan # cut 未満のものは NaN にする (「ch_map<cut」 は True (cut未満) と False で構成される array)
plt.imshow(ch_map, vmin=0, vmax=20, origin="lower left")
plt.show()

# 逆バージョン
ch_map = data[300] # 300チャンネル目をスライス
ch_map[ch_map!=ch_map] = 0.0 # NaN が入っているところに 0.0 を入れる (a!=a は NaN の時だけ Falseが返される)

cut = 1.5
integ_map = data[280:320]# 280-319チャンネルをとってくる
integ_map[integ_map<cut] = np.nan # cut 未満のものは NaN にする (厳密には、-cut<T<cut を NaN にした方が良い)
integ_map = np.nansum(integ_map, axis=0)*hdu.header["CDELT3"]/1000.0 # 積分
plt.imshow(integ_map, vmin=0, vmax=300, origin="lower left")
plt.show()

map 全体に対する計算

例として、励起温度 Tex を求めます (12CO(1-0)が光学的に厚いと仮定)。

cut = 1.5
ch_map = data[300]
ch_map[ch_map<cut] = np.nan
Tex_map = 5.53/(np.log(1.0 + (5.53/(ch_map + 0.84)))) # ちゃんとブロードキャストしてくれるのでこのように大抵一行でかける
plt.imshow(Tex_map, vmin=0, vmax=30, origin="lower left", cmap="jet")
plt.colorbar()
plt.show()

以上です。

リンク
目次