logo标签列表

在裁剪fits文件后，坐标并没有保持不变。

pythonastronomyastropyphotutils
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)

原始图像（左）和裁剪后的图像（右）如下所示：

enter image description here

帧中心星的（赤经，赤纬）赤道坐标为：

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)
1
Gabriel - 啊，好的，谢谢您让我知道。现在当我点击photutils标签时，会出现“ photutils标签没有使用指南，您能帮助我们创建吗？”并且当我点击“详细信息”时，信息很混乱，我不明白您提议了什么。 我不是SO的超级用户，如果您是并且知道该在何处报告，请报告。这似乎是他们可以改进其流程或新建议标签的事情。 - Christoph
我所提供的标签描述正在等待同行审查。我以为你能看到这个。你能在Meta上发布一个关于这个问题的问题吗？http://meta.stackoverflow.com/ - Gabriel
@Gabriel，不要将解决方案编辑到您的问题中，请发布一个包含您的解决方案（例如Cutout2D）的答案，并再次从问题中删除解决方案部分。您甚至可以接受它作为答案！ - MSeifert
我为了完整起见添加它们。Christoph可以将代码添加到他自己的答案中，这是最好解决问题的答案。 - Gabriel
3个回答
3
您IP地址为143.198.54.68，由于运营成本限制，当前对于免费用户的使用频率限制为每个IP每72小时10次对话，如需解除限制，请点击左下角设置图标按钮（手机用户先点击左上角菜单按钮）。 - Gabriel
请阅读Cutout2D的文档字符串。它说你必须单独传递wcs，即from astropy.wcs import WCS; cutout = Cutout2D(data=hdu.data, wcs=WCS(hdu.header)。然后查看输出的cutout对象。根据文档，它应该有一个新的data和wcs属性供您剪裁。如果需要，您可以再次将wcs转换为标头，并从数据和标头制作HDU。这有帮助吗？我现在离线了，但如果需要，明天可以发布工作代码。实际上，在Cutout2D文档中应该有一个WCS示例，有人应该发送一个拉取请求。 - Christoph
在测试中有两个使用Cutout2D和wcs的例子：https://github.com/astropy/astropy/blob/master/astropy/nddata/tests/test_utils.py#L441 这些例子应该被加入文档中，但现在对你可能会有所帮助。 - Christoph
我已经在我的问题中添加了两种实现方法，使用cutout_footprint和Cutout2D。你的解决方案确实更简单，所以我将接受这个答案。谢谢你们两个！ - Gabriel
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。

这将坐标更接近原始值，但仍然有偏差：12:10:37 +39:25:13。您尝试过这个解决方案吗？它是否按预期工作？同时会触发一个警告：WARNING: FITSFixedWarning: RADECSYS='FK5' / WORLD COORDINATE FRAME the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs] WARNING: FITSFixedWarning: 'datfix' made the change 'Changed '13/03/95' to '1995-03-13''. [astropy.wcs.wcs] - Gabriel
我认为原因是photutils在使用cutout_footprint时采用了(x, y)约定。您可以尝试在切片WCS时使用[280-50: 280+50，267-25: 267+25]吗？ - MSeifert
是的，就是这样。Cutout2D类看起来也很有趣，我会去看看的。谢谢MSeifert！ - Gabriel
1
@Gabriel 我刚刚发现了一些奇怪的事情，你的中心是 (x, y) 但你的盒子是 (y, x)。我不确定这是否正确。 - MSeifert
MSeifert没错。这显然是cutout_footprint的工作原理。如果你注意到我在原问题中右侧的裁剪图像，你会看到裁剪是使用(y, x)框应用的。一开始这看起来很奇怪。也许我会在该项目的Github上询问。 - Gabriel
1
@MSeifert - 我已经发布了一个答案，建议直接使用Cutout2D。手动进行WCS更改更加复杂且容易出错，特别是像这里讨论的那样切换x/y。 Cutout2D有所帮助，尽管您仍然需要阅读文档以正确获取x / y顺序。 - Christoph
1
由于需要使用额外的软件包，这里提供另一种解决方案：ccdproc，特别是基于astropy.nddata.NDData的CCDData类。它可以简化文件读取过程。
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，只有ccdproc.CCDData中的read和write部分是显式实现的。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接