3

考虑以下代码（下载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)
4

3 回答 3

3

photutils.utils.cutout_footprint只剪切像素，它不会更新用于在像素和世界坐标之间转换的 WCS。

请改用astropy.nddata.utils.Cutout2D。

于 2016-11-26T23:13:04.140 回答
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.nddatanamed中有一个方便的功能，默认情况下Cutout2D会更新WCS。

于 2016-11-26T18:37:34.197 回答
1

另一个答案，因为它需要一个额外的包：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，只有readandwrite部分在ccdproc.CCDData

于 2016-11-28T19:46:51.927 回答