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

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) экваториальные координаты звезды в центре кадра:

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]

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

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, wcs, header а также 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