Эта проблема беспокоит меня уже некоторое время. Я пытаюсь обработать большой объем данных в виде файла .fits (порядка 11000x9000 пикселей). Что мне нужно сделать, так это создать график координат RA / Dec (в идеале с использованием astropy.wcs) для многих объектов на небе с контурами из одного файла соответствия и шкалой серого (или тепловой картой из другого).
Моя проблема в том, что всякий раз, когда я разрезаю данные изображения (в интересующую меня область), я теряю связь с координатами неба. Это означает, что нарезанное изображение находится не в правильном месте.
Я адаптировал пример из документации по астропии для избавьте вас от боли в моих данных. (Примечание: я хочу, чтобы контуры покрывали большую площадь, чем изображение, какое бы решение для этого ни было, оно должно работать с обоими данными)
Вот код, с которым у меня возникли проблемы:
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)
