Introduction
Pythonでお手軽にHEALPixを描画するためのツールである、healpyの使い方をご説明します。
インストール
pip install
healpyはPyPIに登録されているため、pipで簡単にインストールすることができます。
pip install healpy
とりあえず描画
スクリプト
以下のスクリプトを適当なエディタでコピペし、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を使用しています。
nside = 32
でpixel数を指定します。実際に使われるpixel数はN_pixel = 12 * nside ** 2です。
.fitsファイルデータを描画してみる
Planckのデータを用いる
適当なデータをPlanckのページからダウンロードしてきます。今回はLFI_SkyMap_030_1024_R2.01_full.fitsを使用しました。
.fitsファイルデータの描画
先ほどダウンロードした.fitsファイル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()
すると以下のようなウィンドウが出てきます。
銀河座標から赤道座標への変換
もし銀河座標でなく赤道座標を用いて描画したい場合は、以下のようにcoordオプションを付けます。
# make healpix map
hp.mollview(planck_map, coord=['G','E'], title='Example of Planck',
unit='mK', norm='hist', min=-1, max=1)
すると結果は以下のようになります。
カラーを変更する
カラーマップを変更したい場合は、以下のようにcmapオプションでカラーマップを指定します。
# make healpix map
hp.mollview(planck_map, cmap='jet', coord=['G','E'],
title='Example of Planck', unit='mK', norm='hist', min=-1, max=1)
すると以下のようにカラーマップが変化します。好みのカラーマップで描画を楽しみましょう。
違うデータセットを描画する
このPlanckのデータはI_Stokes, Q_Stokes, U_Stokesなどの違うデータが一緒くたになった.fitsファイルデータです。デフォルトで何も設定しなければ、最初のデータセットであるI_Stokesの結果が描画されます。違うデータセットを描画したい場合には、以下のように設定します。以下はQ_Stokesの描画例です。
# read Planck data
planck_map = hp.read_map('./LFI_SkyMap_030_1024_R2.01_full.fits', ('Q_Stokes', ))
結果は以下のようになります。
どのデータがひとまとめになっているかは.fitsファイルをダウンロードしてきたPlanckのページを参照してください。
球面調和関数で展開する
healpyは観測結果を球面調和関数で展開することができます。以下のスクリプトを実行してみてください。
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()
すると以下のような結果が出てきます。
結果を描画するだけでなく解析までできてしまう、優れたモジュールです。
.csvファイルのデータをHEALPix上に描画する
healpyは.fitsファイルデータ以外のデータも描画することができます。
.csvファイルデータの描画
データはKaggleのPLAsTiCCデータを使用しました。HEALPix上に、各データ点の赤方偏移をカラーでプロットします。
スクリプト
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を利用しています。結果は以下のようになります。
hp.projscatterのs=5でプロットのサイズを、c=df_gal_lb.hostgal_speczでカラープロットを行うデータを指定しています。
宇宙物理でHEALPixを使った描画が必要になりましたら、ぜひPythonのhealpyをお使いください。