为FITS图像分配WCS坐标

5
我一直在疯狂搜索文档,但没有找到答案。 我正在使用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中分配它们,但没有成功。我只需要一个简单直接的方法来将这些坐标分配给图像。

感谢您提供的任何帮助。

1个回答

5
我建议你开始使用astropy。对于你的项目,astropy.wcs包可以帮助你编写FITS WCS头文件,而astropy.io.fits API基本上与你现在使用的pyfits相同。此外,帮助页面非常好,我将翻译他们的WCS构建页面以匹配你的示例。
关于你的问题:FITS文件并不会为每个像素打上坐标标签。我想创建一个像素查找表之类的东西是可能的,但实际的WCS是将X,Y像素算法转换为天体测量学坐标(在你的情况下是“银河系”)的算法。这里有一个很好的页面
我要给你指引的例子在这里:

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)

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