python - 如何叠加两个具有不同 WCS/分辨率的 .fits 图像？

Question

我正在尝试从一个物体的图像中绘制轮廓，假设在波段 A 中，在同一物体的较低分辨率图像的顶部，假设在波段 Z 中。

两张图片都在大文件中，所以我必须为每张图片创建一个 2D 剪切图。然而，带 Z 中的图像比带 A 的像素少得多。因此，当我尝试将它们绘制在同一张图中时，带 Z 中的图像根据像素尺寸绘制，因此在左下角非常小图的(见下面链接中的图)。我需要两个图像都根据天文角度尺寸而不是像素来绘制。

import numpy as np
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
from astropy.nddata import Cutout2D
from reproject import reproject_interp

# Read in .fits images and WCS information
image_a = fits.open("image_a.fits")
image_z = fits.open("image_z.fits")
wcs_a = WCS(image_a.header).celestial
wcs_z = WCS(image_z.header).celestial

# Define object RA and DEC
obj_coords = SkyCoord(ra=135.19081*u.degree, dec=0.68991393*u.degree, frame='fk5')

# Create 2D cutouts of the object in each band in a 6 by 6 arcsec box
size = u.Quantity((6, 6), u.arcsec)
stamp_a = Cutout2D(image_a.data, obj_coords, size, wcs=wcs_a)
stamp_z = Cutout2D(image_z.data, obj_coords, size, wcs=wcs_z)

# Plot the the stamp in band Z with contours from band A
ax = plt.subplot(projection = wcs_z)
plt.imshow(stamp_z.data, cmap='gray', origin='lower', alpha=1)
ax.contour(stamp_a.data, levels = [2*sigma,3*sigma,4*sigma,5*sigma], 
colors='green')

Resulting plot (click me)

我还看到了有关使用重新投影的内容。我试过了，但我不认为我完全理解它，因为它只会产生一个很小的像素。这是我尝试过的:

array, footprint = reproject_interp((stamp_a.data, wcs_a), image_z.header)
plt.imshow(array, origin='lower')

如果我缺少任何信息来回答问题，请告诉我。在使用了几年的 IDL 之后，我上周才切换到 Python，所以如果它有点粗糙，我深表歉意。

Answer 1

您使用方式的问题reproject是你传递了 (stamp_a.data, wcs_a)，但是 wcs_a 是来自原始图像的 WCS，而不是来自 stamp。

您可以从 Cutout2D 中获取与您的图章匹配的 WCS 对象图像。我认为更改为 (stamp_a.data, stamp_a.wcs) 会给你一个正确的结果。

你看过astropy.visualisation.wcsaxes了吗？然而？你可以看一个例子 here如何绘制图像，以及如何从具有不同 WCS 的第二幅图像中绘制轮廓。如果这对您有用，它将比手动使用 Cutout2D 和 reproject 更简单。

我不确定性能；如果您的图像很大，您仍然可能想像现在一样进行剪裁。如果内存/执行时间不是问题，您可以只绘制整个图像，然后使用 wcsaxes 设置您想要的范围(请参阅他们的教程文档)。

是自己重新投影还是让 matplotlib 和 astropy.visualisation.wcsaxes 在幕后进行的问题也是一个品味问题。直接使用 reproject 的代码会更多一些也更复杂，但可以让您更好地控制所使用的精确重投影方法(例如插值或较慢的通量守恒 reproject_exact)和可能更容易完全理解正在发生的事情。