裁剪后坐标不保存fi

2024-09-27 09:32:13 发布

您现在位置:Python中文网/ 问答频道 /正文

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

原始(左)和裁剪(右)图像如下所示:

enter image description here

位于画面中央的恒星的(ra,dec)equatorial coordinates是:

^{pr2}$

为什么裁剪框中的坐标不同?在


这是两种解决问题的方法,使用两种不同的方法。在

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)

Tags: tocropimporthdrfileheaderwcsfits
3条回答

{{a1}特别是基于另一个cda2}的包:

它简化了文件的读取:

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

必须指定单元,因为头单元与astropy.units模块不兼容。在

关于CCDData的重要一点是,您(大多数情况下)不需要亲自处理unitswcsheader和{},它们被保存为属性,最基本的操作会相应地保存(和更新)它们。例如切片:

^{pr2}$

此切片ccd_cropped已更新WCS信息,因此您可以继续处理它(如再次保存):

ccd_cropped.write('cropped.fits')

它应该有正确更新的坐标。切片部分实际上也可以使用^{},只有read和{}部分在ccdproc.CCDData中显式实现

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

请改用astropy.nddata.utils.Cutout2D。在

坐标发生了变化,因为您切片了图像,但没有更改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复制到您的页眉:

^{pr2}$

保存文件之前。在

我不确定cutout_footprint的作用,所以在创建wcs_cropped时,您可能需要更改切片索引。在

astropy.nddata中有一个名为^{}的便利功能,它在默认情况下更新WCS。在

相关问题 更多 >

    热门问题