天文データ解析入門 その12 (Python: よく使う list や numpy array の操作)

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

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

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

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

# list の生成

```python
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 の生成

```python
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 の操作

```python
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 の操作

```python
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

```python
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

```python
ch_map = data[300] # 300チャンネル目をスライス
plt.imshow(ch_map, vmin=0, vmax=20, origin="lower left")
plt.show()
```
![ダウンロード.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/3e491e80-0a4a-ce8e-90ff-666e689e240d.png)

### integrated intensity map

```python
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()
```
![ダウンロード (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/4c4afe53-7282-c91c-f231-3bac64f354d4.png)
[参考](https://qiita.com/Shinji_Fujita/items/d09f0bef16994bb0ec38)

### peak intensity map

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

![ダウンロード (2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/ffd1b22c-bbe7-34fe-bd71-6b69b8765a64.png)

### position–velocity diagram

```python
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()
```
![ダウンロード (3).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/add4f10f-9c99-56f1-88d9-62f25f5d83e3.png)
[参考](https://qiita.com/Shinji_Fujita/items/ac6351c0749d2f8db1ff)

## map への mask

```python
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()
```
![ダウンロード (4).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/a256344d-ab89-b4c5-b979-f39648e812b3.png)

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

```python
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()
```
![ダウンロード (6).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/bb73d087-5e37-79b1-75e0-ac79d4fd9b24.png)

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

```python
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()
```
![ダウンロード (5).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/117814/390224b5-8744-6e4d-deae-44368b3253b3.png)

以上です。

リンク
[目次](https://qiita.com/Shinji_Fujita/items/ecb869e5b6a9fa468483)