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