2016-08-18 3 views
1

私は.fitsファイル(天文画像)を使用してプロットを作りたいと私は、彼らが関連していると思う二つの問題を経験しています:astropyからこの例を使用してプロットライン

を:

from matplotlib import pyplot as plt 
from astropy.io import fits 
from astropy.wcs import WCS 
from astropy.utils.data import download_file 

fits_file = 'http://data.astropy.org/tutorials/FITS-images/HorseHead.fits' 
image_file = download_file(fits_file, cache=True) 
hdu = fits.open(image_file)[0] 
wcs = WCS(hdu.header) 

fig = plt.figure() 
fig.add_subplot(111, projection=wcs) 
plt.imshow(hdu.data, origin='lower', cmap='cubehelix') 
plt.xlabel('RA') 
plt.ylabel('Dec') 
plt.show() 

私はこの画像を生成することができます。

enter image description here

今、私は同じ座標を使用して、いくつかの点をプロットしたいと思いますしかし

plt.scatter(85, -2, color='red') 

、私はこれを行う:画像として

enter image description here

私は、ピクセルcoordinantesでplotingています。さらに、イメージはもはやフレームのサイズに一致しません(座標はうまく見えますが)

これらの問題に対処する方法についてのアドバイスはありますか?

答えて

3

与えられた座標をプロットするのはとても簡単です。あなたがする必要があるのは、transformを適用することだけです。

あなたの例をコピーし、私が何かを変更した理由とその理由を追加しました。

from matplotlib import pyplot as plt 
from astropy.io import fits 
from astropy.wcs import WCS 
from astropy.utils.data import download_file 

fits_file = 'http://data.astropy.org/tutorials/FITS-images/HorseHead.fits' 
image_file = download_file(fits_file, cache=True) 

# Note that it's better to open the file with a context manager so no 
# file handle is accidentally left open. 
with fits.open(image_file) as hdus: 
    img = hdus[0].data 
    wcs = WCS(hdus[0].header) 

fig = plt.figure() 

# You need to "catch" the axes here so you have access to the transform-function. 
ax = fig.add_subplot(111, projection=wcs) 
plt.imshow(img, origin='lower', cmap='cubehelix') 
plt.xlabel('RA') 
plt.ylabel('Dec') 

# Apply a transform-function: 
plt.scatter(85, -2, color='red', transform=ax.get_transform('world')) 

そして結果は次のとおりです。

enter image description here

注あなたはキャンバスのみ画像の領域を表示したい場合は、単にその後再び制限を適用することを:

# Add a scatter point which is in the extend of the image: 
plt.scatter(85.3, -2.5, color='red', transform=ax.get_transform('world')) 

plt.ylim(0, img.shape[0]) 
plt.xlim(0, img.shape[1]) 

をそれは以下を与える:

enter image description here

サイドノートもここにあります。 AstroPyは素晴らしいユニットサポートを持っているので、arcminsとarcsecsを度に変換するのではなく、単に "ユニット"を定義することができます。あなたはまだかかわらず、変換が必要です。

from astropy import units as u 
x0 = 85 * u.degree + 20 * u.arcmin 
y0 = -(2 * u.degree + 25 * u.arcmin) 
plt.scatter(x0, y0, color='red', transform=ax.get_transform('world')) 

enter image description here

+1

うわー!返信とすべての偉大なアドバイスをありがとう! – Delosari

関連する問題