将 WCS 坐标分配给 FITS 图像

Question

我一直在疯狂搜索文档，但找不到此文档的答案。

我在 python 中生成 FITS 图像，需要将 WCS 坐标分配给图像。我知道有很多方法可以通过将点源与已知目录进行匹配来做到这一点，但在这种情况下，我正在生成尘埃图，因此点源匹配不起作用（据我所知）。

所以图像是一个形状为 (240,240) 的 2D Numpy 数组。它是这样写的（x 和 y 坐标分配有点奇怪，它以某种方式工作）：

H, xedges, yedges = np.histogram2d(glat, glon, bins=[ybins, xbins], weights=Av)
count, x, y = np.histogram2d(glat, glon, bins=[ybins, xbins])
H/=count
hdu = pyfits.PrimaryHDU(H)
hdu.writeto(filename)

>>> print H.shape
(240,240)

这一切都很好。分配银河坐标似乎你需要做的就是：

glon_coords = np.linspace(np.amin(glon), np.amax(glon), 240)
glat_coords = np.linspace(np.amin(glat), np.amax(glat), 240)

但是我不明白 FITS 图像是如何存储这些坐标的，所以我不知道怎么写。我也尝试在 SAO DS9 中分配它们，但没有运气。我只需要一种直接将这些坐标分配给图像的方法。

感谢您提供的任何帮助。

Answer 1

我建议您开始使用astropy。对于您的项目，astropy.wcs 包可以帮助您编写 FITS WCS 标头，并且astropy.io.fits API 与您现在使用的 pyfits 基本相同。此外，帮助页面非常好，我要做的就是翻译他们的 WCS 构建页面以匹配您的示例。

对于您的问题：FITS 不会用坐标“标记”每个像素。我想可以创建一个像素查找表或类似的东西，但actual WCS 是一个algorithmic X，Y 像素到天体坐标的转换（在你的情况下是“银河”）。 A nice page is here.

我要给你举的例子在这里：

http://docs.astropy.org/en/latest/wcs/index.html#building-a-wcs-structure-programmatically

这是我为您的项目提供的未经测试的伪代码：

# untested code

from __future__ import division # confidence high

# astropy
from astropy.io import fits as pyfits
from astropy import wcs

# your code
H, xedges, yedges = np.histogram2d(glat, glon, bins=[ybins, xbins], weights=Av)
count, x, y = np.histogram2d(glat, glon, bins=[ybins, xbins])
H/=count

# characterize your data in terms of a linear translation from XY pixels to 
# Galactic longitude, latitude. 

# lambda function given min, max, n_pixels, return spacing, middle value.
linwcs = lambda x, y, n: ((x-y)/n, (x+y)/2)

cdeltaX, crvalX = linwcs(np.amin(glon), np.amax(glon), len(glon))
cdeltaY, crvalY = linwcs(np.amin(glat), np.amax(glat), len(glat))

# wcs code ripped from 
# http://docs.astropy.org/en/latest/wcs/index.html

w = wcs.WCS(naxis=2)

# what is the center pixel of the XY grid.
w.wcs.crpix = [len(glon)/2, len(glat)/2]

# what is the galactic coordinate of that pixel.
w.wcs.crval = [crvalX, crvalY]

# what is the pixel scale in lon, lat.
w.wcs.cdelt = numpy.array([cdeltX, cdeltY])

# you would have to determine if this is in fact a tangential projection. 
w.wcs.ctype = ["GLON-TAN", "GLAT-TAN"]

# write the HDU object WITH THE HEADER
header = w.to_header()
hdu = pyfits.PrimaryHDU(H, header=header)
hdu.writeto(filename)