Координаты не сохраняются после кадрирования припадки файл

Question

Рассмотрим следующий код (скачать test.fits):

from astropy.io import fits 
from photutils.utils import cutout_footprint 

# Read fits file. 
hdulist = fits.open('test.fits') 
hdu_data = hdulist[0].data 
hdulist.close() 

# Some center and box to crop 
xc, yc, xbox, ybox = 267., 280., 50., 100. 
# Crop image. 
hdu_crop = cutout_footprint(hdu_data, (xc, yc), (ybox, xbox))[0] 
# Add comment to header 
prihdr = hdulist[0].header 
prihdr['COMMENT'] = "= Cropped fits file") 
# Write cropped frame to new fits file. 
fits.writeto('crop.fits', hdu_crop, prihdr) 

Оригинал (слева) и обрезано (справа) изображения выглядят следующим образом:

(ra, dec) equatorial coordinates для звезды в центре кадра:

Original frame: 12:10:32 +39:24:17 
Cropped frame: 12:12:07 +39:06:50 

Почему координаты разные в обрезанной кадре?

---

Это два способа решить эту проблему, используя два разных метода.

from astropy.io import fits 
from photutils.utils import cutout_footprint 
from astropy.wcs import WCS 
from astropy.nddata.utils import Cutout2D 
import datetime 

# Read fits file. 
hdulist = fits.open('test.fits') 
hdu = hdulist[0].data 
# Header 
hdr = hdulist[0].header 
hdulist.close() 

# Some center and box to crop 
xc, yc, xbox, ybox = 267., 280., 50., 100. 

# First method using cutout_footprint 
# Crop image. 
hdu_crop = cutout_footprint(hdu, (xc, yc), (ybox, xbox))[0] 
# Read original WCS 
wcs = WCS(hdr) 
# Cropped WCS 
wcs_cropped = wcs[yc - ybox // 2:yc + ybox // 2, xc - xbox // 2:xc + xbox // 2] 
# Update WCS in header 
hdr.update(wcs_cropped.to_header()) 
# Add comment to header 
hdr['COMMENT'] = "= Cropped fits file ({}).".format(datetime.date.today()) 
# Write cropped frame to new fits file. 
fits.writeto('crop.fits', hdu_crop, hdr) 

# Second method using Cutout2D 
# Crop image 
hdu_crop = Cutout2D(hdu, (xc, yc), (xbox, ybox), wcs=WCS(hdr)) 
# Cropped WCS 
wcs_cropped = hdu_crop.wcs 
# Update WCS in header 
hdr.update(wcs_cropped.to_header()) 
# Add comment to header 
hdr['COMMENT'] = "= Cropped fits file ({}).".format(datetime.date.today()) 
# Write cropped frame to new fits file. 
fits.writeto('crop.fits', hdu_crop.data, hdr) 

Answer 1

photutils.utils.cutout_footprint только вырезает пиксели, он не обновляет WCS, который используется для преобразования координат между пикселями и мирами.

Используйте вместо этого astropy.nddata.utils.Cutout2D.

Answer 2

Координаты изменены, потому что вы нарезали изображение, но не изменили информацию WCS (особенно значения опорных пикселей).

Один из способов будут использовать astropy.WCS:

from astropy.wcs import WCS 
wcs = WCS(hdulist[0].header) 
wcs_cropped = wcs[280-50 : 280+50 , 267-25 : 267+25] 

, а затем скопировать эти обновленные туалеты в ваш заголовок:

prihdr.update(wcs_cropped.to_header()) 

перед сохранением файла.

Я не уверен, что такое cutout_footprint, возможно, вам нужно изменить индексы среза при создании wcs_cropped.

Существует функциональность удобства в astropy.nddata с именем Cutout2D, которая по умолчанию обновляет WCS.

Answer 3

Другой ответ, потому что он требует дополнительного пакета: ccdproc особенно класс CCDData, который основан на astropy.nddata.NDData:

Это упрощает чтение файла:

from ccdproc import CCDData 
ccd = CCDData.read(filename, unit='erg/cm**2/s/AA') 

Агрегат должен быть указан, так как блок заголовка не совместим с модулем astropy.units.

Важно о CCDData является то, что вы (в основном) не нужно заботиться о units, wcsheader, и masks сами, они сохраняются как атрибуты и самые основные операции сохранения (и обновление) их соответственно.Например, нарезка:

xc, yc, xbox, ybox = 267., 280., 50., 100. 
ccd_cropped = ccd[int(yc - ybox // 2) : int(yc + ybox // 2), 
        int(xc - xbox // 2) : int(xc + xbox // 2)] 

Это нарезанный ccd_cropped уже обновил информацию WCS, так что вы можете просто продолжать его обработку (например, сохранить его еще раз):

ccd_cropped.write('cropped.fits') 

И он должен иметь правильно обновленные координаты. Отрезая часть на самом деле можно также с помощью astropy.nddata.NDDataRef, только read и write часть реализуются явно в ccdproc.CCDData