概要

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

install

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

pip install astropy-healpix

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

スクリプト

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

#!/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()

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

スクリプト詳細

fits.open

以下の行でPlanckの全天データmapを読み込んでいます。

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

[1].header

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

「astropy.ioでfitsファイルの情報を取得する」を参考に、N_sideやOrderingなどのfitsファイルの情報を取得します。ここでTTYPE1を選択しているので、I_Stokesのデータをプロットしていることになります。

HEALPix()

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

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

SkyCoord()

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

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

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

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

interpolate_bilinear_skycoord()

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の強みかなという気がします。