PythonでHEALPixを描画する

# Introduction

Pythonでお手軽に[HEALPix](https://healpix.jpl.nasa.gov)を描画するためのツールである、[healpy](https://healpy.readthedocs.io/en/latest/)の使い方をご説明します。

# インストール

## pip install

healpyは[PyPI](https://pypi.org)に登録されているため、pipで簡単にインストールすることができます。

```
pip install healpy
```

# とりあえず描画

## スクリプト

以下のスクリプトを適当なエディタでコピペし、sample.pyとして保存してください。

```python:sample.py
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt 
 
 
# define nside
nside = 32 

# make healpix map
m = np.arange(hp.nside2npix(nside))
hp.mollview(m, title='example')
 
# Draw
plt.show()
```

これを

```
$ python sample.py
```

として実行します。すると以下のようなウィンドウが現れます。描画にはmatplotlibを使用しています。

![スクリーンショット 2019-03-11 18.36.18.png](https://qiita-image-store.s3.amazonaws.com/0/335400/932413cf-9dd1-9481-eb48-c1fd91198503.png)

```python
nside = 32
```

でpixel数を指定します。実際に使われるpixel数は```N_pixel = 12 * nside ** 2```です。

# .fitsファイルデータを描画してみる

## Planckのデータを用いる

適当なデータを[Planckのページ](https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/previews/LFI_SkyMap_030_1024_R2.01_full/index.html)からダウンロードしてきます。今回はLFI_SkyMap_030_1024_R2.01_full.fitsを使用しました。

## .fitsファイルデータの描画

先ほどダウンロードした.fitsファイルPython実行ディレクトリに移動させ、以下のスクリプトを実行します。

```python

import healpy as hp
import matplotlib.pyplot as plt 
 

# read Planck data
planck_map = hp.read_map('./LFI_SkyMap_030_1024_R2.01_full.fits')
 
# make healpix map
hp.mollview(planck_map, title='Example of Planck', 
            unit='mK', norm='hist', min=-1, max=1)
 
# Draw
plt.show()

```

すると以下のようなウィンドウが出てきます。

![スクリーンショット 2019-03-12 10.48.10.png](https://qiita-image-store.s3.amazonaws.com/0/335400/0b4eead5-1c2e-1a0c-a6bb-191266dd4e74.png)

## 銀河座標から赤道座標への変換

もし銀河座標でなく赤道座標を用いて描画したい場合は、以下のように`coord`オプションを付けます。

```python
# make healpix map
hp.mollview(planck_map, coord=['G','E'], title='Example of Planck', 
            unit='mK', norm='hist', min=-1, max=1)
```

すると結果は以下のようになります。

![スクリーンショット 2019-03-12 10.52.24.png](https://qiita-image-store.s3.amazonaws.com/0/335400/9fea28ff-c80d-e78e-a6bf-58ae033ba80f.png)

## カラーを変更する

カラーマップを変更したい場合は、以下のように`cmap`オプションでカラーマップを指定します。

```python
# make healpix map
hp.mollview(planck_map, cmap='jet', coord=['G','E'], 
            title='Example of Planck', unit='mK', norm='hist', min=-1, max=1)
```

すると以下のようにカラーマップが変化します。好みのカラーマップで描画を楽しみましょう。

![スクリーンショット 2019-03-12 12.26.00.png](https://qiita-image-store.s3.amazonaws.com/0/335400/4acf8249-0c27-3fe3-8f56-03a32772d5b7.png)


## 違うデータセットを描画する

このPlanckのデータはI_Stokes, Q_Stokes, U_Stokesなどの違うデータが一緒くたになった.fitsファイルデータです。デフォルトで何も設定しなければ、最初のデータセットであるI_Stokesの結果が描画されます。違うデータセットを描画したい場合には、以下のように設定します。以下はQ_Stokesの描画例です。

```python
# read Planck data
planck_map = hp.read_map('./LFI_SkyMap_030_1024_R2.01_full.fits', ('Q_Stokes', ))
```

結果は以下のようになります。

![スクリーンショット 2019-03-12 18.40.17.png](https://qiita-image-store.s3.amazonaws.com/0/335400/0a673234-76af-e2fd-daee-a4b81e7cdc07.png)

どのデータがひとまとめになっているかは.fitsファイルをダウンロードしてきたPlanckのページを参照してください。

# 球面調和関数で展開する

healpyは観測結果を球面調和関数で展開することができます。以下のスクリプトを実行してみてください。

```python

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt 
 

# read Planck data
planck_map = hp.read_map('./LFI_SkyMap_030_1024_R2.01_full.fits')
 
# Spherical harmonic transforms
LMAX = 1024
cl = hp.anafast(planck_map, lmax=LMAX)

# set ell list
ell = np.arange(len(cl))

# plot spectrum
plt.plot(ell, ell*(ell+1)*cl)
plt.xlabel('ell')
plt.ylabel('ell(ell+1)cl')
plt.grid()
plt.show()

```

すると以下のような結果が出てきます。

![スクリーンショット 2019-03-12 11.08.33.png](https://qiita-image-store.s3.amazonaws.com/0/335400/75cf3618-7dae-bb70-5098-92bb4ffc7764.png)

結果を描画するだけでなく解析までできてしまう、優れたモジュールです。

# .csvファイルのデータをHEALPix上に描画する

healpyは.fitsファイルデータ以外のデータも描画することができます。

## .csvファイルデータの描画

データは[KaggleのPLAsTiCCデータ](https://www.kaggle.com/c/PLAsTiCC-2018)を使用しました。HEALPix上に、各データ点の赤方偏移をカラーでプロットします。

## スクリプト

```python

import healpy as hp
import pandas as pd
import matplotlib.pyplot as plt


# read data
df_gal_lb = pd.read_csv('./training_set_metadata.csv')

# set grid lines
hp.graticule()

# draw scatter plot in HEALPix
obj = hp.projscatter(df_gal_lb.gal_l, df_gal_lb.gal_b, s=5,
                     c=df_gal_lb.hostgal_specz, cmap='jet', lonlat=True, coord='G')

# draw colorbar
cbar = plt.colorbar(obj, orientation='horizontal')
cbar.set_label('red shift', size=14)

# show image
plt.show()
```

データの読み込みにpandasを利用しています。結果は以下のようになります。

![スクリーンショット 2019-03-12 12.35.29.png](https://qiita-image-store.s3.amazonaws.com/0/335400/2dde7059-ad1a-012d-3da0-7fcb083f8a76.png)

`hp.projscatter`の`s=5`でプロットのサイズを、`c=df_gal_lb.hostgal_specz`でカラープロットを行うデータを指定しています。

宇宙物理でHEALPixを使った描画が必要になりましたら、ぜひPythonのhealpyをお使いください。