Astropyで天文データ解析 [1] -FITSキューブの基本操作-
はじめに
本記事ではFITSデータの取り扱い方法をAstropy公式ドキュメントのWorking with FITS-cubesを参考にしながら、記事としてまとめる。なお、筆者の備忘録として活用できるように公式ドキュメントのみでは理解しにくかった部分を掘り下げ、コードに関する要点や、天文学の基礎用語の解説など様々な点に配慮した。
AstropyでFITSキューブを操作する
必要なライブラリ
| 今回必要なライブラリ | どんなライブラリ？ |
|---|---|
| astroquery | データを検索したり、DLしたりする |
| spectral_cube | FITSデータキューブを表示する |
| reproject | データの分解能に合わせて座標を投影しなおす |
導入方法
Anaconda環境ならcondaで導入する。
conda install -c conda-forge hoge
conda-forgeが指定されており、condaのdefaultsチャンネルでは入手できないはず。
FITSデータキューブをダウンロードする
astropy.utils.download_fileを使って、手軽にFITSファイルを取得する。download_file()はurl先のファイルをキャッシュに追加してDLされる。ミラーから取得したり (sources)、ローカルファイルからimportしたりもできる。 (import_file_to_cache, import_download_cache)。
公式のデータが用意されており、これをdownload_file()で読み込む。
FITSデータキューブをダウンロードする
from astropy.utils.data import download_file
hi_datafile = download_file(
'http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits',
cache=True, show_progress=True)
FITSのDLが完了した。データはHI4Pi Surveyの中性水素21cm線データらしい。
FITSデータキューブの表示
spectral_cubeを使えば、FITSデータキューブの表示や操作ができる。
FITSデータキューブを変数に格納
astropy.io.fitsとspectral_cubeを使って、簡単にFITSデータキューブを変数に格納。
FITSデータキューブを取得
import astropy.io.fits as iofits
from astropy.utils.data import download_file
from spectral_cube import SpectralCube
# download_fileでurl先ファイルをキャッシュに追加してFITSデータキューブをDL
hi_datafile = download_file(
'http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits',
cache=True, show_progress=True)
# astropy.io.fitsオブジェクトが持つopenメソッドでFITSデータキューブのPrimary HDUをopen
hi_data = iofits.open(hi_datafile)
# SpectralCubeのreadメソッドでFITSデータキューブ情報を取得
cube = SpectralCube.read(hi_data)
print(cube)
出力結果
SpectralCube with shape=(450, 150, 150) and unit=K:
n_x: 150 type_x: GLON-TAN unit_x: deg range: 286.727203 deg: 320.797623 deg
n_y: 150 type_y: GLAT-TAN unit_y: deg range: -50.336450 deg: -28.401234 deg
n_s: 450 type_s: VRAD unit_s: m / s range: -598824.534 m / s: 600409.133 m / s
出力結果を見ると、このデータは(450,150,150)で単位が温度(K)の3次元データキューブだとわかる。また、n_x, n_yはGLON, GLATすなわちGalactic longitude, Galactic latitudeで、銀河座標のこと。また、n_sはVRADすなわちradial velocity(視線速度)のことを表す。
FITSデータキューブのクイック表示
SpectralCubeのquicklookメソッドでFITSデータキューブを簡易表示する。
FITSデータキューブのクイック表示
SpecralCube.read(hoge)でFITSデータキューブを取得した後、次のようにquicklook()でデータを表示できる。
cube[300, :, :].quicklook()
見ての通り、cubeは3次元ndarrayオブジェクトであり、300枚目を指定したければ、[300, :, :]とすればよい。
出力結果
出力結果では以下のwarningとINFOが出てくる。
INFO: Auto-setting vmin to -2.037e+00 [aplpy.core]
INFO: Auto-setting vmax to 2.123e+01 [aplpy.core]
これを見ると、どうやらSpectralCubeはAPLpyを使って表示していることがわかる。APLpyについては以下参照。論文ライクな図を作るためには必須のライブラリで、matplotlibベースで天文向けの図を作成するモジュールが整備されている。
また、公式ドキュメントの出力図はcmap='viridis'になっているようだが、筆者の環境ではグレースケールになった。matplotlibの出力はデフォルトでviridisになるはずだが、原因は不明。
また、視線速度方向のスペクトルを取得したければn_s, n_x, n_yのうちn_sを全選択して、n_x, n_yで空間方向のプロット範囲を決める。
FITSキューブデータの視線速度方向のスペクトル取得
cube[:, 75, 75].quicklook()
上の例では、n_sを全選択し、75ピクセル四方の範囲で視線速度方向のスペクトルを取得している。出力結果は次のようになる。
縦軸が温度(K), 横軸が視線速度(m/s)を表している。
FITSデータキューブをx,y軸でスライス
FITSデータキューブを異なるx,y軸でスライス
得られる図はいわゆる位置-速度図(position-velocity diagram)になる。恐らくSpecralCubeで簡単に表示できるのだろうが、その方法はまだ調べていないので、それ以外の方法でスライスできないか考えてみる。for文を使うのはナンセンス。
思いついたのは、numpy.rot90()を使って、3次元FITSキューブデータの配列ごと回転してしまう方法。x軸にスライスしたFITSキューブデータが欲しければ、予めastropy.io.fitsでFITSデータをndarrayオブジェクトに変換し、
import astropy.io.fits as iofits
import numpy as np
hdu = iofits.open('hoge.fits')[0]
fits_data = hdu.data
numpy.rot90()で配列を回転。hoge.fitsはなんでもいい。ここでは、公式ドキュメントのreduced_TAN_C14.fitsをそのまま使う。
rot_data = np.rot90(fits_data, k=-1, axes=(0,1))
np.rot90()の扱い方に混乱しかけたので、メモ。3次元配列の(v,y,x)に対して、(0,1,2)が対応する。kはデフォルトで1になっており、配列を90°回転させた回数を指定。axes=(a,b)の部分は、第1, 2引数で指定した平面で回転する。
これを用いて、l-v図すなわち縦軸: 銀経, 横軸: 視線速度の方向に回転させてみる。
通常はx-y平面の天体画像がv方向に並んでいるのだから、こんな感じでプロットできる。
import astropy.io.fits as iofits
from astropy.utils.data import download_file
import matplotlib.pyplot as plt
import numpy as np
hi_datafile = download_file(
'http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits',
cache=True, show_progress=True)
fits_data = iofits.open(hi_datafile)[0].data
# v-x平面で回転
vy_diagram = np.rot90(fits_data, axes=(0,1))
# x-y平面で回転
vx_diagram = np.rot90(vy_diagram, axes=(1,2))
# matplotlibで図を出力。見栄えを整える。
plt.imshow(vx_diagram[100])
plt.xlabel(r'$V_\mathrm{rad}$')
plt.ylabel('Galactic longitude')
plt.colorbar(shrink=0.4)
出力結果
うまくいったようだ。ただ、座標は物理座標のままだし、カラーバーもshrinkで無理やり縮めているので、あまりスマートな方法じゃない。
異なる視線速度方向のスペクトルを表示する
cube[:, a, b]のa,bの部分を変えるだけで、指定した空間範囲の視線速度方向のスペクトルを表示することができる。
# 例
cube[:, 50, 50]
出力結果
空間座標(50, 50)のスペクトルが得られた。
サブキューブを作る
特にデータサイズが大きいFITSならば、欲しい空間範囲に制限して、データサイズを節約したほうがよいだろう。まずは、Galactic longitudeとlatitudeの範囲を指定し、その後cubeオブジェクトが持つsubcubeメソッドを使って、サブキューブを作る。
欲しい空間範囲だけ切り取り、サブキューブを作る
空間範囲を指定してサブキューブを作る
空間座標の範囲指定をするときに、astropy.unitsで単位を取得して、座標に単位(degree)を与える。
import astropy.units as u
# lon. lat.の範囲を指定する
lat_range = [-46, -40] * u.deg
lon_range = [306, 295] * u.deg
# 指定範囲の座標でサブキューブを作る(list型で指定)
x0, x1 = lon_range[0], lon_range[1]
y0, y1 = lat_range[0], lat_range[1]
#xとyの最大最小値を指定。(lowとhigh)
sub_cube = cube.subcube(xlo=x0, xhi=x1, ylo=y0, yhi=y1)
print(sub_cube)
出力結果
SpectralCube with shape=(450, 38, 57) and unit=K:
n_x: 57 type_x: GLON-TAN unit_x: deg range: 294.113498 deg: 306.009028 deg
n_y: 38 type_y: GLAT-TAN unit_y: deg range: -46.014280 deg: -40.027398 deg
n_s: 450 type_s: VRAD unit_s: m / s range: -598824.534 m / s: 600409.133 m / s
速度範囲を指定してサブキューブを作る
sub_cubeを作った状態で、sub_cubeオブジェクトが持つ.spectral_slab()メソッドを使って、速度範囲を指定して切り取る。slabとは、幅の広い厚い板のことらしく、FITSデータキューブの速度範囲を指定してサブキューブを作ることは、厚切りのパンをスライスするイメージに近いだろう。
astropy.unitsでkmとsを取得してkm/sの単位を与え、速度範囲の始めと終わりを指定する。
sub_cube_slab = sub_cube.spectral_slab(-300. *u.km / u.s, 300. *u.km / u.s)
print(sub_cube_slab)
出力結果
SpectralCube with shape=(226, 38, 57) and unit=K:
n_x: 57 type_x: GLON-TAN unit_x: deg range: 294.113498 deg: 306.009028 deg
n_y: 38 type_y: GLAT-TAN unit_y: deg range: -46.014280 deg: -40.027398 deg
n_s: 226 type_s: VRAD unit_s: m / s range: -299683.842 m / s: 301268.441 m / s
astropy.unitsで単位指定しない方法
astropy.unitsでいちいち単位指定するのは面倒だから、cubeオブジェクトの持つ.world()メソッドを使って、座標情報を配列として抜き出す。例えば、FITSデータキューブのl,bを取得したければ、以下のようにする。
_, b, _ = cube.world[0, 10:50, 0]
_, _, l = cube.world[0, 0, :]
いきなりアンダースコアがでてきて面食らったが、アンダースコアの部分は座標を抽出しませんという意味だろう。現に変数にアンダースコアを指定するときは、必要の値や一時的な変数として定義される。bをprint()で出力してみると、
[-45.99408714 -45.84754016 -45.70088965 -45.55413761 -45.40728605
-45.26033698 -45.11329243 -44.96615443 -44.81892502 -44.67160624
-44.52420014 -44.37670876 -44.22913418 -44.08147844 -43.93374364
-43.78593183 -43.63804509 -43.49008552 -43.3420552 -43.19395622
-43.04579068 -42.89756068 -42.74926834 -42.60091574 -42.45250502
-42.30403829 -42.15551766 -42.00694526 -41.85832321 -41.70965365
-41.5609387 -41.4121805 -41.26338119 -41.1145429 -40.96566777
-40.81675795 -40.66781558 -40.51884281 -40.36984178 -40.22081464] deg
最後にdegがあるので、単位付きの銀緯情報が得られる。だから、銀経のほうも適切にスライスして与えれば、似たようなことができるはず。ただ、どちらにせよ最大最小を指定しないといけないので、組み込み関数maxやminを使って、xloやxhiの値を指定する必要がある。かえって煩雑になる可能性があるため、こういう用途で.world()メソッドは使わないほうが良いのかもしれない。
積分強度図とは？
.moment()メソッドを使えば、簡単に積分強度図を作ることができる。momentとは、大規模なデータセットを平均値や標準偏差などの重要な情報に圧縮することらしい。通常のデジタルカメラはRGBの3色合成が得られるが、電波天文などでは、フルでスペクトルを取得することができ、FITSデータキューブの情報を圧縮する。天体データがガウス分布に従うなら、ヒストグラムとして圧縮できる。また、FITSデータキューブを速度方向に足し合わせて(積分して)いけば、ガスの質量すなわちmoment0マップを得ることができる。すなわち、moment0とは、
のこと(
さらに、ガスがどこに移動しているか知るため(我々太陽系に遠ざかっているか/近づいているか)、moment1マップを得ることができる。moment1マップとは、スペクトルの強度で重み付けした速度のことだろう。
つまり、FITSキューブをそのまま処理しようとするのではなく、データをmoment0, 1マップに落とし込むことで、何が起きているかを理解することができる。[1]
2次以上のmonentマップ
moment2マップでは、速度分散(velocity dispersion)やスペクトル線の線幅を得ることができるらしい。また、高次のmomentNマップは、次の式で定義できる。分光軸単位で線幅マップを得るには、Linewidth mapsを参照とのこと。
積分強度図と柱密度
.moment()を使って、積分強度図を描いてみる。momentメソッドは、cubeオブジェクトだけでなく、subcubeオブジェクトの属性でもある。
astropy.unitsでmomentマップを取得する
import astropy.io.fits as iofits
import astropy.units as u
from astropy.utils.data import download_file
from spectral_cube import SpectralCube
hi_datafile = download_file(
'http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits',
cache=True, show_progress=True)
hi_data = iofits.open(hi_datafile)
cube = SpectralCube.read(hi_data)
hi_data.close()
# moment0マップを取得する
moment_0 = cube.with_spectral_unit(u.km/u.s).moment(order=0)
# moment1マップを取得する
moment_1 = cube.with_spectral_unit(u.km/u.s).moment(order=1)
print(f'moment_0マップの単位: {moment_0.unit}')
print(f'moment_1マップの単位: {moment_1.unit}')
出力結果
moment_0マップの単位: K km / s
moment_1マップの単位: km / s
moment0を柱密度に変換する
FITSキューブの1枚1枚のチャンネルマップは、奥行き方向の情報を持たないので、視線速度方向にどのぐらい物質が重なっているかを評価する指標を用いる。それが柱密度 (column density)だ。柱密度の導出は、spatial densityやabundance, 運動温度など基本的な物理量を導き出すための第一歩である。具体的には、エネルギー準位
光学的に薄い場合、HIの柱密度は、
hi_column_density = moment_0 * 1.82 * 10**18 / (u.cm * u.cm) * u.s / u.K / u.km
print(hi_column_density)
出力結果
[[1.81440277e+20 1.74025162e+20 1.87063660e+20 ... 2.19830017e+20
2.17239365e+20 2.10436630e+20]
[1.81522339e+20 1.76137711e+20 1.83742249e+20 ... 2.25064026e+20
2.19466385e+20 2.12471756e+20]
[1.90567719e+20 1.90308060e+20 1.86109756e+20 ... 2.22132950e+20
2.26488892e+20 2.17063980e+20]
...
[6.41453735e+20 6.53187988e+20 6.81179016e+20 ... 8.62090088e+20
8.57158630e+20 8.34916870e+20]
[7.09448242e+20 7.19683411e+20 7.29619873e+20 ... 8.46022034e+20
8.62227051e+20 8.40449890e+20]
[7.21431946e+20 7.29688049e+20 7.50491089e+20 ... 8.48804749e+20
8.52107849e+20 8.54180786e+20]] 1 / cm2
となる。単位を見ると、柱密度は面密度と同じ物理次元を持つことがわかる。
momentマップの保存と表示
作成したmoment_0, 1オブジェクトは、.write()メソッドでFITSファイルとして保存することができる。
積分強度図の保存
moment_0 = cube.with_spectral_unit(u.km/u.s).moment(order=0)
moment_1 = cube.with_spectral_unit(u.km/u.s).moment(order=1)
# momentマップをFITSファイルとして保存
moment_0.write('moment_0.fits')
moment_1.write('moment_1.fits')
また、作成したmoment_0, 1オブジェクトをmatplotlibで綺麗に表示してみる。単にplt.imshow()で出力すると、座標が物理座標だから良くない。そこで、matplotlib.pyplotが持つadd_subplotの引数にprojectionで座標を指定する。
projectionを指定してmoment0, 1マップを表示
まずはprojecttionに座標情報を与えるため、moment_1オブジェクトが持つ.wcs()メソッドで座標を取得する。
m1_wcs = moment_1.wcs
print(m1_wcs)
出力結果。FITSのキーワードが表示される。
WCS Keywords
Number of WCS axes: 2
CTYPE : 'GLON-TAN' 'GLAT-TAN'
CRVAL : 303.75 -40.0
CRPIX : 75.55762081784387 77.3723930634757
PC1_1 PC1_2 : 1.0 0.0
PC2_1 PC2_2 : 0.0 1.0
CDELT : -0.1494444443846667 0.1538888888273333
NAXIS : 0 0
公式に沿って、図を出力するプログラムを書く。具体的には、moment_1マップに対して、fk5で赤道座標をオーバーレイして、さらにHI柱密度のコントアを重ねている。
import astropy.io.fits as iofits
import astropy.units as u
from astropy.utils.data import download_file
import matplotlib.pyplot as plt
from spectral_cube import SpectralCube
hi_datafile = download_file(
'http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits',
cache=True, show_progress=True)
hi_data = iofits.open(hi_datafile)
cube = SpectralCube.read(hi_data)
hi_data.close()
# astropy.unitsでmoment0マップの単位を取得する
moment_0 = cube.with_spectral_unit(u.km/u.s).moment(order=0)
# astropy.unitsでmoment1マップの単位を取得する
moment_1 = cube.with_spectral_unit(u.km/u.s).moment(order=1)
print(f'moment_0マップの単位: {moment_0.unit}')
print(f'moment_1マップの単位: {moment_1.unit}')
# moment_0をcolumn densityに変換
hi_column_density = moment_0 * 1.82 * 10**18 / (u.cm * u.cm) * u.s / u.K / u.km
print(hi_column_density)
# figsizeを決める
fig = plt.figure(figsize=(18,12))
# subplotにmoment_1の座標指定
ax = fig.add_subplot(111, projection=moment1.wcs)
# moment_1マップを表示(.hdu.dataメソッドでFITSデータ取得)
im = ax.imshow(moment_1.hdu.data, cmap='RdBu_r', vmin=0, vmax=200)
ax.invert_yaxis()
ax.set_xlabel("Galactic Longitude (degrees)", fontsize=16)
ax.set_ylabel("Galactic Latitude (degrees)", fontsize=16)
cbar = plt.colorbar(im, pad=.07)
cbar.set_label('Velocity (km/s)', size=16)
# 赤道座標をオーバーレイする
overlay = ax.get_coords_overlay('fk5')
overlay.grid(color='white', ls='dotted', lw=2)
overlay[0].set_axislabel('Right Ascension (J2000)', fontsize=16)
overlay[1].set_axislabel('Declination (J2000)', fontsize=16)
# axオブジェクトの.contour()メソッドで柱密度のコントアを重ねる
levels = (1e20, 5e20, 1e21, 3e21, 5e21, 7e21, 1e22)
ax.contour(hi_column_density.hdu.data, cmap='Greys_r', alpha=0.5, levels=levels)
出力結果
公式はサブキューブのslabデータで図を出力しているが、筆者はcubeごと積分して図を出力した。
また、moment_1の部分をmoment_0に変えれば、積分強度図も表示できる。なお、下図のcmapは'jet'で、カラーバーのラベルをIntegrated Intensityに変更している。
位置-速度図を表示
位置-速度図は文字通り、経度と速度の情報のみを持つダイアグラムのこと。さっきは無理やりnp.rot90()で回転させて表示したのだが、ようやくここで解答を学べる。
FITSデータキューブを銀緯軸でスライスする
銀緯のインデックスをint型で指定して、それをfig.add_subplotのslices引数に指定することで、簡単に位置-速度図を描けるようだ。冒頭のimport部分は省略するが、コードはこんな感じ。
# lat.のindexを指定
lat_slice = 18
# slicesにindexを指定し、projectionに座標情報を与える
fig = plt.figure(figsize=(18,12))
ax = fig.add_subplot(111, projection=cube.wcs, slices=('y', lat_slice, 'x'))
# index = 18の位置-速度図を.imshow()で表示
im = ax.imshow(cube[:, lat_slice, :].transpose().data, cmap='jet')
ax.invert_yaxis()
ax.set_xlabel("LSR Velocity (m/s)", fontsize=16)
ax.set_ylabel("Galactic Longitude (degrees)", fontsize=16)
# カラーバーのorientation引数にhorizontalを与えて、バーを水平にする
cbar = plt.colorbar(im, pad=.07, orientation='horizontal')
cbar.set_label('Temperature (K)', size=16)
出力結果
両端は非数データだから、空白になっている。公式の図で余白がない理由は、サブキューブのslabデータだからだろう。
slabデータからSMCだけ分離して位置-速度図を描く
公式からの宿題。新しいスペクトルslabデータを作って、SMC(小マゼラン雲(Small Magellanic Cloud))だけ分離し、異なるインデックス(銀経軸)で位置-速度図を作成する問題。
b v
今度はSMCを銀経軸でスライスして、位置-速度図を作成する。まずは、復習のためsub_cube_slabを作成していく。
import astropy.io.fits as iofits
import astropy.units as u
from astropy.utils.data import download_file
import matplotlib.pyplot as plt
from spectral_cube import SpectralCube
hi_datafile = download_file(
'http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits',
cache=True, show_progress=True)
hi_data = iofits.open(hi_datafile)
cube = SpectralCube.read(hi_data)
hi_data.close()
# sub_cubeの作成
lat_range = [-46, -40] * u.deg
lon_range = [306, 295] * u.deg
x0, x1 = lon_range[0], lon_range[1]
y0, y1 = lat_range[0], lat_range[1]
sub_cube = cube.subcube(xlo=x0, xhi=x1, ylo=y0, yhi=y1)
# sub_cube_slabの作成
sub_cube_slab = sub_cube.spectral_slab(-300. *u.km / u.s, 300. *u.km / u.s)
# 銀経軸でスライスして位置-速度図を描く
# lon.のindexを指定
lon_slice = 18
fig = plt.figure(figsize=(18,12))
ax = fig.add_subplot(111, projection=sub_cube_slab.wcs, slices=('lon_slice', 'y', 'x'))
im = ax.imshow(sub_cube_slab[:, :, lon_slice].transpose().data)
ax.invert_yaxis()
ax.set_xlabel("LSR Velocity (m/s)", fontsize=16)
ax.set_ylabel("Galactic Longitude (degrees)", fontsize=16)
# カラーバーのorientation引数にhorizontalを与えて、バーを水平にする
cbar = plt.colorbar(im, pad=.07, orientation='horizontal')
cbar.set_label('Temperature (K)', size=16)
出力結果
銀経軸でスライスして、銀緯-速度図を得ることができた。
astroqueryでHerschelデータを取得
Herschel 350 micronとHI emissionのデータを比較して、ダストをトレースしたい。
astroquery.ESASkyクラスで、Herschelデータを簡単に取得できる。
astroqueryでESAのHerschelデータを読み込む
import astropy.units as u
from astroquery.esasky import ESASky
result = ESASky.query_region_maps('SMC', radius=1*u.deg, missions='Herschel')
print(result)
出力結果
TableList with 1 tables:
'0:HERSCHEL' with 15 column(s) and 28 row(s)
TableListにどのような情報があるかを知るには、HERSCHELキーを指定すればよい。
print(result['HERSCHEL'].keys())
出力結果
['dec_deg', 'duration', 'end_time', 'filter', 'fov', 'instrument', 'observation_id', 'observation_oid', 'observing_mode_name', 'postcard_url', 'product_url', 'ra_deg', 'start_time', 'stc_s', 'target_name']
result['HERSCHEL']から350 micron画像だけをダウンロードするために、観測時のフィルター情報を取得する必要があるようだ。そこで、result['HERSCHEL']に['filter']を設定し、このフィルター情報を見てみると、
result['HERSCHEL']['filter']
出力結果
filter
microns
-------------
70, 160
250, 350, 500
100, 160
70, 160
70, 160
70, 160
70, 160
250, 350, 500
100, 160
70, 160
...
100, 160
70, 160
250, 350, 500
70, 160
250, 350, 500
100, 160
70, 160
70, 160
70, 160
250, 350, 500
Length = 28 rows
一部省略されているが、250, 350, 500 micronの3つで観測したデータがある。
astroqueryでHerschel 350 micronデータを取得
350 micron画像以外はboolean配列を用いてマスクする。
# HERSCHELクエリを文字列に変換
filters = result['HERSCHEL']['filter'].astype(str)
# filterにforループでbooleanマスクをかける
mask = np.array(['250, 350, 500' == s for s in filters], dtype='bool')
# 350 micronデータ以外をマスクしてTableListを再作成
target_obs = TableList({"HERSCHEL":result['HERSCHEL'][mask]})
350 mircronデータのTableListを作成することができたら、後はESASkyオブジェクトのget_mapsメソッドを使って、TableListからデータをダウンロードするだけ。
IR_images = ESASky.get_maps(target_obs)
しかし、ここで躓いた。get_mapsを実行すると、EOFError: Ran out of inputのエラーが発生してしまうのである。そこで、get_maps()の引数にcache=Falseを与え、リモート先のurlデータをキャッシュに追加しない設定に書き換えた。
出力結果。データのダウンロードは15分程かかった。
INFO: Starting download of HERSCHEL data. (5 files) [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342198590 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8634359&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342205055 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8614152&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342198565 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8613787&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342198566 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8634358&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342205092 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8614195&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading of HERSCHEL data complete. [astroquery.esasky.core]
ダウンロードした350 micronデータのIR_images変数オブジェクト情報を見るために、.info()を指定する。
IR_images['HERSCHEL'][0]['350'].info()
出力結果
Filename: Maps/HERSCHEL/anonymous1637497344/hspirepmw401_25pxmp_0110_m7303_1342198565_1342198566_1462476888800.fits.gz
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 184 ()
1 image 1 ImageHDU 47 (2407, 2141) float64
2 error 1 ImageHDU 47 (2407, 2141) float64
3 coverage 1 ImageHDU 47 (2407, 2141) float64
4 History 1 ImageHDU 23 ()
5 HistoryScript 1 BinTableHDU 39 84R x 1C [326A]
6 HistoryTasks 1 BinTableHDU 46 65R x 4C [1K, 27A, 1K, 9A]
7 HistoryParameters 1 BinTableHDU 74 450R x 10C [1K, 20A, 13A, 196A, 1L, 1K, 1L, 74A, 11A, 41A]
図に出力するために、IR_imagesの座標情報とPrimaryHDUを取得する。
herschel_header = IR_images['HERSCHEL'][0]['350']['image'].header
herschel_wcs = WCS(IR_images['HERSCHEL'][0]['350']['image']) # Extract WCS information
herschel_imagehdu = IR_images['HERSCHEL'][0]['350']['image'] # Extract Image data
print(herschel_wcs)
出力結果
WCS Keywords
Number of WCS axes: 2
CTYPE : 'RA---TAN' 'DEC--TAN'
CRVAL : 18.99666058287459 -71.82876374906319
CRPIX : 1063.0 1517.0
NAXIS : 2407 2141
Herschel 350 micronデータを表示
取得したIR_imagesのWCS情報、PrimaryHDUをもとにして、matplolib用いてHerschel 350 micronデータを表示する。
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
# NaNは0とみなしてマスクをかける
# himage_nan_locs = np.isnan(herschel_imagehdu.data)
herschel_data_nonans = herschel_imagehdu.data
# herschel_data_nonans[himage_nan_locs] = 0
# WCS情報をprojectionに与える
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=herschel_wcs)
# 350 micronデータを.imshow()で表示. (matplotlib.colorsのLogNormでログスケールにする)
im = ax.imshow(herschel_data_nonans, cmap='viridis',
norm=LogNorm(), vmin=2, vmax=50)
ax.set_xlabel("Right Ascension", fontsize = 16)
ax.set_ylabel("Declination", fontsize = 16)
ax.grid(color = 'white', ls = 'dotted', lw = 2)
cbar = plt.colorbar(im, pad=.07)
# herschel_headerから単位を指定し、カラーバーのラベルに表示
cbar.set_label(''.join(['Herschel 350'r'$\mu$m ','(', herschel_header['BUNIT'], ')']), size = 16)
# get_coords_overlay()で銀河座標グリッドをオーバーレイ
overlay = ax.get_coords_overlay('galactic')
overlay.grid(color='r', ls='dotted', lw=1)
overlay[0].set_axislabel('Galactic Longitude', fontsize=14)
overlay[1].set_axislabel('Galactic Latitude', fontsize=14)
出力結果
周りにランダムノイズみたいなのがある。これは何だろうか？
どうやら、ピクセル値が0の場所をLogNorm()で出力すると、当該表示になるようだ。確かにlog0は定義できないので、LogNorm()で表すとこうなるのだろう。
HI 21cmのコントアを重ねる
WCSAxes.get_transform()で座標変換することができる。Herschel 350 micronデータにHI 21cm のコントアを重ねてみる。
Herschel 350 micronデータにHI 21cmコントアを重ねる
# projectionにHerschelのWCS情報を与える
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=herschel_wcs)
# momentマップを.imshow()で表示
im = ax.imshow(herschel_data_nonans, cmap='viridis',
norm=LogNorm(), vmin=5, vmax=50, alpha=.8)
ax.set_xlabel("Right Ascension", fontsize=16)
ax.set_ylabel("Declination", fontsize=16)
ax.grid(color = 'white', ls='dotted', lw=2)
x_lim = ax.get_xlim()
y_lim = ax.get_ylim()
cbar = plt.colorbar(im, fraction=0.046, pad=-0.1)
cbar.set_label(''.join(['Herschel 350'r'$\mu$m ','(', herschel_header['BUNIT'], ')']), size=16)
# 銀河座標のグリッドをオーバーレイ
overlay = ax.get_coords_overlay('galactic')
overlay.grid(color='r', ls='dotted', lw=1)
overlay[0].set_axislabel('Galactic Longitude', fontsize=14)
overlay[1].set_axislabel('Galactic Latitude', fontsize=14)
hi_transform = ax.get_transform(hi_column_density.wcs)
# 柱密度コントアを重ねる(透明度alpha=0.8)
# コントアレベルを決める
levels = (2e21, 3e21, 5e21, 7e21, 8e21, 1e22)
ax.contour(hi_column_density.hdu.data, cmap='Greys_r', alpha=0.8, levels=levels,
transform=hi_transform)
# moment1マップを重ねる(透明度alpha=0.5)
im_hi = ax.imshow(moment_1.hdu.data, cmap='RdBu_r', vmin=0, vmax=200, alpha=0.5, transform=hi_transform)
cbar_hi = plt.colorbar(im_hi, orientation='horizontal', fraction=0.046, pad=0.07)
cbar_hi.set_label('HI 'r'$21$cm Mean Velocity (km/s)', size=16)
ax.set_xlim(x_lim)
ax.set_ylim(y_lim)
出力結果
データの分解能に合わせて再投影
reprojectを使えば、画像の分解能を落として再投影できる。しかも、データのフラックス値を保存して情報を失わずに再投影できるようだ。HerschelデータはHI 21cmよりも高い空間分解能持つので、正確に強度やフラックスを比較するには、reprojectでHerschelデータの分解能を落とし、画像を再投影する必要がある。
まずは両者の分解能を調べてみる。
HerschelとHI 21cmデータの分解能を調べる
画像の分解能とはすなわち参照点(各座標の画素)の増分のことだから、FITSヘッダからCDELTキーワードを調べてみればよい。CDELT1がx方向の増分、CDELT2がy方向の増分を表す。
print('Herschel Resolution (dx,dy) = ', herschel_header['cdelt1'], herschel_header['cdelt2'])
print('HI Resolution (dx,dy) = ', hi_column_density.hdu.header['cdelt1'], hi_column_density.hdu.header['cdelt1'])
出力結果
Herschel Resolution (dx,dy) = -0.002777777777778 0.002777777777777778
HI Resolution (dx,dy) = -0.14944444438467 -0.14944444438467
出力結果を見てみると、Herschelデータの方が10倍以上の空間分解能を持っていることがわかる。
reproject_interpで分解能を落とす
reproject_interpを用いれば、分解能を落とし再投影することができる。具体的には、reproject_interpに再投影するHDUオブジェクトと、WCS情報を含むヘッダを与えてやる。
rescaled_herschel_data, _ = reproject_interp(herschel_imagehdu,
hi_column_density.hdu.header)
rescalingしたデータを新規のHDUに格納する。
rescaled_herschel_imagehdu = fits.PrimaryHDU(data = rescaled_herschel_data,
header = hi_column_density.hdu.header)
ここで取得したrescaled_herschel_imagehduは、HI 21cmデータに合わせて分解能が落ちている。
分解能を落としたデータを再投影
rescaled_herschel_imagehduが取得できれば、後はこれをプロットするだけ。
# NaNがあったらboolean maskをかけて0にする
image_nan_locs = np.isnan(rescaled_herschel_imagehdu.data)
rescaled_herschel_data_nonans = rescaled_herschel_imagehdu.data
rescaled_herschel_data_nonans[image_nan_locs] = 0
# rescaleしたHerschelデータヘッダから、WCS情報をprojectionに指定
fig = plt.figure(figsize = (18,12))
ax = fig.add_subplot(111,projection = WCS(rescaled_herschel_imagehdu))
# rescalingしたデータを表示 (表示はLogNorm())
im = ax.imshow(rescaled_herschel_data_nonans, cmap = 'viridis',
norm = LogNorm(), vmin = 5, vmax = 50, alpha = .8)
#im = ax.imshow(rescaled_herschel_imagehdu.data, cmap = 'viridis',
# norm = LogNorm(), vmin = 5, vmax = 50, alpha = .8)
ax.set_xlabel("Galactic Longitude", fontsize = 16)
ax.set_ylabel("Galactic Latitude", fontsize = 16)
ax.grid(color = 'white', ls = 'dotted', lw = 2)
x_lim = ax.get_xlim()
y_lim = ax.get_ylim()
cbar = plt.colorbar(im, fraction=0.046, pad=-0.1)
cbar.set_label(''.join(['Herschel 350'r'$\mu$m ','(', herschel_header['BUNIT'], ')']), size = 16)
# 赤道座標のグリッドをオーバーレイ
overlay = ax.get_coords_overlay('fk5')
overlay.grid(color='black', ls='dotted', lw = 1)
overlay[0].set_axislabel('Right Ascension', fontsize = 14)
overlay[1].set_axislabel('Declination', fontsize = 14)
hi_transform = ax.get_transform(hi_column_density.wcs)
# 柱密度マップを重ねる
levels = (2e21, 3e21, 5e21, 7e21, 8e21, 1e22)
ax.contour(hi_column_density.hdu.data, cmap = 'Greys_r', alpha = 0.8, levels = levels,
transform = hi_transform)
# moment1マップを重ねる
im_hi = ax.imshow(moment_1.hdu.data, cmap = 'RdBu_r', vmin = 0, vmax = 200, alpha = 0.5, transform = hi_transform)
# HIの平均速度カラーバーを追加
cbar_hi = plt.colorbar(im_hi, orientation = 'horizontal', fraction=0.046, pad=0.07)
cbar_hi.set_label('HI 'r'$21$cm Mean Velocity (km/s)', size = 16)
ax.set_xlim(x_lim)
ax.set_ylim(y_lim)
出力結果
復習: チャレンジ問題を解く
ページ最下部のChallengeに取り組む。ここでは、reprojectとWCSを使って異なる座標投影でマッピングしなおし、どのように画像をresamplingするかが目的だろう。銀河座標表示のHI 21cmデータをHerschelデータの赤道座標情報を使ってresamplingしてみる。復習のため細かく説明を与えた。
まずは必要となるものをimport。前述したように、HIデータをDLするためにastroqueryが必要で、FITSデータキューブを扱うためにspectral_cubeが必要。さらに、データをreprojectするには、reprojectも必要。また、積分強度図の単位を指定するため、astropy.unitをimportしておく。全部importすると以下のようになる。
from matplotlib.colors import LogNorm
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from astropy.utils.data import download_file
from astropy.wcs import WCS
import astropy.units as u
from astropy.utils.data import download_file
from astroquery.esasky import ESASky
from astroquery.utils import TableList
from reproject import reproject_interp
from spectral_cube import SpectralCube
まずは、HIデータをdownload_fileでDLしHDUを取得しておく。今まで通りSMCのデータでいいだろう。
hi_datafile = download_file(
'http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits',
cache=True, show_progress=True)
hi_data = fits.open(hi_datafile)
cube = SpectralCube.read(hi_data)
hi_data.close()
sub_cube_slabデータを生成する。まずは適切な銀河座標範囲を指定しsub_cubeを作成する。その後、SMCを含む速度範囲のみを抽出してslabデータを作成する。
# 座標範囲の指定
lat_range = [-46, -40] * u.deg
lon_range = [306, 295] * u.deg
# サブキューブのx,yの最大最小を指定
sub_cube = cube.subcube(xlo=lon_range[0], xhi=lon_range[1], ylo=lat_range[0], yhi=lat_range[1])
# サブキューブをslabデータにする
sub_cube_slab = sub_cube.spectral_slab(-300. *u.km / u.s, 300. *u.km / u.s)
次に、sub_cube_slabからmoment0マップを作成しておく。また、astropy.unitsで積分強度図の単位も取得しておく。今回はカラーバーを付けるつもりはないが、カラーバーを付けるときはその都度手入力するのではなく、astropy.unitsでの単位を指定するほうがよいだろう。
moment_0 = sub_cube_slab.with_spectral_unit(u.km/u.s).moment(order=0)
astropy.unitsで積分強度の単位を取得してみる。
print(f'Moment_0の単位: {moment_0.unit}')
出力結果
Moment_0の単位: K km / s
さらにESAのHerschelデータを取得する。取得にはastroqueryを用いればよい。なお、350 micronデータを取得する場合、HERSCHELキーワードで指定されたクエリにboolean maskをかけ、350 micronデータのみ取得すればよい。取得したmask配列をTableListに渡す。
result = ESASky.query_region_maps('SMC', radius=1*u.deg, missions='Herschel')
filters = result['HERSCHEL']['filter'].astype(str)
mask = np.array(['250, 350, 500' == s for s in filters], dtype='bool')
target_obs = TableList({"HERSCHEL":result['HERSCHEL'][mask]})
astroqueryのESASky.get_mapsを使ってデータをダウンロードする。ダウンロードには時間がかかる。(約15分)
download_dirには各自のディレクトリを指定すること。未指定の場合、カレントディレクトリにMapsディレクトリが生成されHerschel FITSデータがDLされる。
IR_images = ESASky.get_maps(target_obs, download_dir='/hoge/piyo/', cache=False)
出力結果
INFO: Starting download of HERSCHEL data. (5 files) [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342198590 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8634359&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342205055 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8614152&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342198565 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8613787&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342198566 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8634358&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading Observation ID: 1342205092 from http://archives.esac.esa.int/hsa/whsa-tap-server/data?RETRIEVAL_TYPE=STANDALONE&observation_oid=8614195&DATA_RETRIEVAL_ORIGIN=UI [astroquery.esasky.core]
INFO: [Done] [astroquery.esasky.core]
INFO: Downloading of HERSCHEL data complete. [astroquery.esasky.core]
INFO: Maps available at /hoge/piyo/. [astroquery.esasky.core]
データがDLできたら、HerschelのImage HDUを取得する。
herschel_imagehdu = IR_images['HERSCHEL'][0]['350']['image']
ここでデータの非数を0でマスクしておく。
himage_nan_locs = np.isnan(herschel_imagehdu.data)
herschel_data_nonans = herschel_imagehdu.data
herschel_data_nonans[himage_nan_locs] = 0
作成したHIデータの積分強度図とHerschel 350 micronデータをmatplotlibで表示する。
fig = plt.figure(figsize=(15,7))
ax1 = plt.subplot(121, projection=WCS(herschel_imagehdu.header))
ax1.imshow(herschel_data_nonans,norm=LogNorm(),vmin=2,vmax=50)
# .coordsで座標labelを付ける
ax1.coords['ra'].set_axislabel('Right Ascension')
ax1.coords['dec'].set_axislabel('Declination')
ax1.set_title(r'$Herschel\;\;\;350{\mu}m$')
ax2 = plt.subplot(122, projection=moment_0.wcs)
ax2.imshow(moment_0.hdu.data,norm=LogNorm())
ax2.coords['glon'].set_axislabel('Galactic Longitude')
ax2.coords['glat'].set_axislabel('Galactic Latitude')
# 見栄えが良くなるようHIデータのlabelとtickはreverseしておく
ax2.coords['glat'].set_axislabel_position('r')
ax2.coords['glat'].set_ticklabel_position('r')
ax2.set_title('HI 21cm emission')
出力結果
ではreprojectでHIデータを赤道座標にreproject方法は至ってシンプル。reproject_interpの第1引数にreproject対象のHDUを指定し、第2引数にreproject用のWCS情報を含むヘッダを与える。
rescaled_hi_data, footprint = reproject_interp(moment_0.hdu,
# reproject the Herschal image to match the HI data
herschel_imagehdu.header)
あとはmatplotlibを使って図を表示するのみ。
fig = plt.figure(figsize=(20,7))
ax1 = plt.subplot(131, projection=WCS(herschel_imagehdu.header))
ax1.imshow(herschel_imagehdu.data, origin='lower', norm=LogNorm(), vmin=2, vmax=50)
ax1.coords['ra'].set_axislabel('Right Ascension')
ax1.coords['dec'].set_axislabel('Declination')
ax1.set_title(r'$Herschel\;\;\;350{\mu}m$')
ax2 = plt.subplot(132, projection=WCS(herschel_imagehdu.header))
ax2.imshow(rescaled_hi_data, origin='lower', norm=LogNorm())
ax2.coords['ra'].set_axislabel('Right Ascension')
ax2.coords['dec'].set_axislabel('Declination')
ax2.set_title('Reprojected HI data')
ax3 = plt.subplot(133, projection=WCS(herschel_imagehdu.header))
ax3.imshow(footprint, origin='lower',)
ax3.coords['ra'].set_axislabel('Right Ascension')
ax3.coords['dec'].set_axislabel('Declination')
ax3.coords['dec'].set_axislabel_position('r')
ax3.coords['dec'].set_ticklabel_position('r')
ax3.set_title('footprint')
出力結果
左からESAのHerschel 350 micronデータ、reprojectしたHIデータ、reprojectのfootprint。reprojectできたようだ。.subplot()で図を作成しているが、DRYなコードを目指すならあまり良い方法ではない。関数化したほうが解析には便利。このChallenge問題には解答が見当たらないが、恐らくこんな感じだろう。
Discussion