考虑以下代码(下载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)
原始图像(左)和裁剪后的图像(右)如下所示:
帧中心星的(赤经,赤纬)赤道坐标为:
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)