如何修剪.fits圖像並保存世界座標以繪製python的空間圖？

Question

這個問題一直困擾着我一段時間。我正在嘗試處理大量的.fits文件形式的數據（大小爲11000x9000像素）。我需要做的是創建一個放大的RA/Dec座標圖（理想情況下使用astropy.wcs）在天空中的許多對象，以及一個擬合文件的輪廓和灰度（或來自另一個的熱圖）。

我的問題是，無論何時我從圖像中分割數據（到我感興趣的區域），我都會失去與天空座標的關聯。這意味着切片圖像不在正確的位置。

我調整了the astropy docs的一個示例，以節省您對數據的痛苦。 （注：我想輪廓覆蓋比圖像更多的面積，無論解決方案是本應該在這兩個數據工作）

這裏是我遇到的麻煩的代碼：

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) 

Answer 1

好消息是：你可以簡單地切你astropy.WCS，以及使你的任務relativly瑣碎：

... 

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你可能需要做的一些reprojection（例如reproject）