我建议您开始使用astropy。就您的项目而言,astropy.wcs包可以帮助您编写 FITS WCS 标头,并且astropy.io.fits API 与您现在使用的pyfits API 基本相同。此外,帮助页面非常好,我要做的就是翻译他们的 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)