如何在astropy Python中修剪.fits图像并保留世界坐标以进行绘图？

Question

如何在astropy Python中修剪.fits图像并保留世界坐标以进行绘图？

6

这个问题困扰了我一段时间。我正在尝试处理大量数据，这些数据以.fits文件的形式存在（大约为11000x9000像素）。我需要做的是在天空上创建一个“缩放”的RA / Dec坐标图（最好使用astropy.wcs），其中包含来自一个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
import numpy as np

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]
wmap = WCS(hdu.header)
data = hdu.data

fig = plt.figure()
ax1 = fig.add_subplot(121, projection=wmap)
ax2 = fig.add_subplot(122, projection=wmap)
# Scale input image
bottom, top = 0., 12000.
data = (((top - bottom) * (data - data.min())) / (data.max() - data.min())) + bottom


'''First plot'''
ax1.imshow(data, origin='lower', cmap='gist_heat_r')

# Now plot contours
xcont = np.arange(np.size(data, axis=1))
ycont = np.arange(np.size(data, axis=0))
colors = ['forestgreen','green', 'limegreen']
levels = [2000., 7000., 11800.]

ax1.contour(xcont, ycont, data, colors=colors, levels=levels, linewidths=0.5, smooth=16)
ax1.set_xlabel('RA')
ax1.set_ylabel('Dec')
ax1.set_title('Full image')

''' Second plot ''' 
datacut = data[250:650, 250:650]
ax2.imshow(datacut, origin='lower', cmap=cmap)
ax2.contour(xcont, ycont, data, colors=colors, levels=levels, linewidths=0.5, smooth=16)

ax2.set_xlabel('RA')
ax2.set_ylabel('')
ax2.set_title('Sliced image')
plt.show()

我尝试使用切片块的WCS坐标来解决这个问题，但我不确定能否将其传递到任何地方！

pixcoords = wcs.wcs_pix2world(zip(*[range(250,650),range(250,650)]),1)

- FriskyGrub

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- MSeifert · Accepted Answer

好消息是：您可以简单地切割您的astropy.WCS，这使得您的任务相对容易。

...

wmapcut = wmap[250:650, 250:650] # sliced here
datacut = data[250:650, 250:650]
ax2 = fig.add_subplot(122, projection=wmapcut) # use sliced wcs as projection
ax2.imshow(datacut, origin='lower', cmap='gist_heat_r')
# contour has to be sliced as well
ax2.contour(np.arange(datacut.shape[0]), np.arange(datacut.shape[1]), datacut, 
            colors=colors, levels=levels, linewidths=0.5, smooth=16)
...

如果您的文件具有不同的WCS，则可能需要进行一些重投影（例如，请参见reproject）。