astropy_healpix

# 概要

HEALPixを描画するときはhealpyを使う方が良い(参考記事:[「PythonでHEALPixを描画する」](https://qiita.com/github-nakasho/items/6d2e06a3caca1f00f58a))と考えておりましたが、[astropy](https://www.astropy.org/)では現在HEALPixに対応したツールなどが開発中なのだそう。ということで今回はこのツールを使ってデータをプロットすることを目標にした記事を書きます。

# install
PyPIに登録されているので、以下コマンドで簡単にインストールができます。

```
pip install astropy-healpix
```

公式ドキュメントではcondaを使ったインストール方法なども紹介されているので、適宜参照してください。

# スクリプト

今回は以下のスクリプトを使用しました。

```python
#!/usr/bin/env python3

from astropy import units as u
from astropy.coordinates import Galactic
from astropy.coordinates import SkyCoord
from astropy_healpix import HEALPix
from astropy.io import fits 
import numpy as np
import matplotlib.pyplot as plt 


# read Planck data
planck_map = fits.open('./LFI_SkyMap_030_1024_R2.01_full.fits')

# get fits info
nside = planck_map[1].header['NSIDE']
order = planck_map[1].header['ORDERING']
# set map data
ttype = planck_map[1].header['TTYPE1']
tmap = planck_map[1].data[ttype]

hp = HEALPix(nside=nside, order=order, frame=Galactic())

# sample a 360, 180 grid in RA/Dec
ra = np.linspace(180., -180., 360) * u.deg
dec = np.linspace(90., -90., 180) * u.deg
ra_grid, dec_grid = np.meshgrid(ra, dec)

# set up astropy coordinate objects
coords = SkyCoord(ra_grid.ravel(), dec_grid.ravel(), frame='galactic')

# interpolate values
tmap = hp.interpolate_bilinear_skycoord(coords, tmap)
tmap = tmap.reshape((180, 360))

# make a plot of the interpolated temperatures
plt.figure(figsize=(9, 4.5))
im = plt.imshow(tmap, extent=[-180, 180, -90, 90], cmap='jet', aspect='auto', vmin=0, vmax=1e-3)
plt.colorbar(im)
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.show()
```

これを実行すると以下のような図が作成されます。

![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/335400/6402e881-2f36-9a19-7bf8-bb18b8e6a480.png)

## スクリプト詳細

### fits.open

以下の行で[Planckの全天データmap](https://pla.esac.esa.int/)を読み込んでいます。

```python
planck_map = fits.open('./LFI_SkyMap_030_1024_R2.01_full.fits')
```

### [1].header

```python
nside = planck_map[1].header['NSIDE']
order = planck_map[1].header['ORDERING']
ttype = planck_map[1].header['TTYPE1']
```

[「astropy.ioでfitsファイルの情報を取得する」](https://qiita.com/github-nakasho/items/255fb19fb99351dbf92d)を参考に、`N_side`や`Ordering`などのfitsファイルの情報を取得します。ここで`TTYPE1`を選択しているので、I_Stokesのデータをプロットしていることになります。

### HEALPix()

```python
hp = HEALPix(nside=nside, order=order, frame=Galactic())
```

ここでHEALPixを設定しています。`frame=Galactic()`で銀河座標系のデータであることを宣言しています。

### SkyCoord()

```python
coords = SkyCoord(ra_grid.ravel(), dec_grid.ravel(), frame='galactic')
```

ここで、**銀河座標**で`ra_grid`, `dec_grid`の範囲の角度のプロットを設定しています。試しに

```python
coords = SkyCoord(ra_grid.ravel(), dec_grid.ravel(), frame='icrs')
```

にしてみると、以下のように赤道座標でデータが描画されます。

![Figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/335400/982014b3-27cb-437e-2b25-a7be7dc4ed10.png)

### interpolate_bilinear_skycoord()

```python
tmap = hp.interpolate_bilinear_skycoord(coords, tmap)
tmap = tmap.reshape((180, 360))
```

`interpolate_bilinear_skycoord`関数を使って、先ほど設定した`coords`の座標に`ttype`で入れたデータセット`tmap`を変換し、`reshape`関数でデータの大きさを整形しています。

# エラー

## エラー文

astropyのバージョンが古いと、以下のようにinterpolate_bilinear_skycoordでエラーが出てくるようです。

```
Traceback (most recent call last):
  File "example5.py", line 33, in <module>
    tmap = hp.interpolate_bilinear_skycoord(coords, tmap)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy_healpix/high_level.py", line 362, in interpolate_bilinear_skycoord
    return self.interpolate_bilinear_lonlat(lon, lat, values)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy_healpix/high_level.py", line 201, in interpolate_bilinear_lonlat
    return interpolate_bilinear_lonlat(lon, lat, values, order=self.order)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy_healpix/core.py", line 551, in interpolate_bilinear_lonlat
    indices, weights = bilinear_interpolation_weights(lon, lat, nside, order=order)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy_healpix/core.py", line 512, in bilinear_interpolation_weights
    lon = lon.to_value(u.rad)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/astropy/units/quantity.py", line 888, in __getattr__
    self.__class__.__name__, attr))
AttributeError: Longitude instance has no attribute 'to_value'
```

## 解消策

以下のコマンドでastropyをアップデートしたら、エラーが出なくなりました。

```
pip install -U astropy
```

# 結言

宇宙物理の色々な便利ライブラリが使えたりするastropyを使って、fitsを読み込みデータを描画する試みをしました。ただ単にHEALPixを描画するだけであればhealpyに軍配が上がります。しかし、何よりastropyと連携してfitsファイルの読み込み・データの加工ができたりするというのは、astropy_healpixの強みかなという気がします。