python を用いた fits画像をWCS空間で上手に重ねる方法

環境構築

pythonにmatplotlibとastropyがインストールされていればOK。anconda3系が入っていればデフォルトでOK。

銀経銀緯で2つの画像を重なる例

可視光とX線のNGC4051の画像を例として、それぞれ全体像のプロットと、中心付近を拡大したプロットの例を示します。

ポイント

画像のfitsファイルからWCSクラスを生成し、それをmatplotlibのprojectionというオプションに入れるだけ。ただし、これだけだと初期の画像サイズにしか表示されない。
ある場所を拡大したい場合は、画像ファイルと同じように、WCSオブジェクトにカットをかける。サンプルコードの、hriwcs[hricx-hridx:hricx+hridx,hricy-hridx:hricy+hridx] の箇所。これを matplotlib の wcs オプションに与えればよい。
２つの画像を重ねて表示したい場合は、難しい問題になる。その場合は、 WCS座標でイメージをマージしpythonでプロットする方法に簡単にまとめてます。

コードとファイル置き場

qiita_plot_fits_2image.py
NGC4051の可視とX線イメージ
- http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/NGC4051_DSS.fits
- http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/NGC4051_HRI.fits

画像の取り方は、pythonのastroqueryを用いて検索し、skyviewでfits画像を取得する方法を参考にされたい。

サンプルコード

qiita_plot_fits_2image.py


import matplotlib 
matplotlib.use('TkAgg')
import glob 
import matplotlib.pyplot as plt
import sys
from matplotlib.colors import LogNorm
from astropy.wcs import WCS
from astropy.io import fits

name="NGC4051"

hri = glob.glob(name + '*HRI*.fits')
dss = glob.glob(name + '*DSS*.fits')	

if len(hri) == 1 and len(dss) == 1:
	pass
else:
	print("need to store HRI and DSS fits")

F = plt.figure(figsize=(10,8))

hrifilename = hri[0]
hriname = hrifilename.replace(".fits",".png")
hrihdu = fits.open(hrifilename)[0]
hriwcs = WCS(hrihdu.header)
hridata = hrihdu.data

hrixlen, hriylen = hridata.shape
hricx = int(0.5 * hrixlen)
hricy = int(0.5 * hriylen)
hridx = int(hrixlen*0.1)
hriwcscut = hriwcs[hricx-hridx:hricx+hridx,hricy-hridx:hricy+hridx]

dssfilename = dss[0]
dssname = dssfilename.replace(".fits",".png")
dsshdu = fits.open(dssfilename)[0]
dsswcs = WCS(dsshdu.header)
dssdata = dsshdu.data

dssxlen, dssylen = dssdata.shape
dsscx = int(0.5 * dssxlen)
dsscy = int(0.5 * dssylen)
dssdx = int(dssxlen*0.1)
dsswcscut = dsswcs[dsscx-dssdx:dsscx+dssdx,dsscy-dssdx:dsscy+dssdx]

plt.figtext(0.45,0.93, name, size="large")
plt.figtext(0.15,0.9, "X-ray, Rosat HRI")
plt.figtext(0.57,0.9, "Optical, UK Shimidt")				

plt.subplot(221, projection=dsswcs)
#		plt.subplot(221, projection=hriwcs)		
try:
	plt.imshow(hridata, origin='lower', norm=LogNorm())
	plt.colorbar()					
except:
	print("ERROR, couldn't plot log-z scale")			
	plt.close()

plt.grid(color='white', ls='solid')
plt.xlabel('Galactic Longitude')
plt.ylabel('Galactic Latitude')

plt.subplot(223, projection=dsswcscut)
#		plt.subplot(223, projection=hriwcscut)		
plt.imshow(hridata[hricx-hridx:hricx+hridx,hricy-hridx:hricy+hridx], origin='lower')
plt.colorbar()		
plt.grid(color='white', ls='solid')
plt.xlabel('Galactic Longitude')
plt.ylabel('Galactic Latitude')


plt.subplot(222, projection=dsswcs)
try:
	plt.imshow(dssdata, origin='lower', norm=LogNorm())
	plt.colorbar()								
except: 
	print("ERROR, couldn't plot log-z scale")
	plt.imshow(dssdata, origin='lower')
	plt.colorbar()		
	
plt.grid(color='white', ls='solid')
plt.xlabel('Galactic Longitude')

plt.subplot(224, projection=dsswcscut)
plt.imshow(dssdata[dsscx-dssdx:dsscx+dssdx,dsscy-dssdx:dsscy+dssdx], origin='lower')
plt.colorbar()		
plt.grid(color='white', ls='solid')
plt.xlabel('Galactic Longitude')

plt.savefig(name  + ".png")
plt.close()

出力結果

上２つの全体像は z 軸を log で表示し、下２つはリニアで表示している。z 軸を log 表示するときにデータが空やゼロしかないとこけるので、try except で回避するようにしている。

参考例

pythonのastroqueryを用いて検索し、skyviewでfits画像を取得する方法を用いて、X線(ROSAT)と可視光(DSS)の画像を取得して作成したgifアニメーションの例。gif アニメを作るには、convert -delay 100 -loop 1 *.png xrayopt.gif で生成した。

X線と可視光 (~200天体) : http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/xrayopt.gif
可視光 (~400天体): http://www-x.phys.se.tmu.ac.jp/~syamada/qiita/qiita_wcs/opt.gif