PythonでHEALPixを描画する

Introduction

Pythonでお手軽にHEALPixを描画するためのツールである、healpyの使い方をご説明します。

インストール

pip install

healpyはPyPIに登録されているため、pipで簡単にインストールすることができます。

pip install healpy

とりあえず描画

スクリプト

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

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をお使いください。