Assign WCS coordinates to a FITS image

Question

I've been searching the documentation like crazy and cannot find an answer for this one.

I'm generating FITS images in python and need to assign WCS coordinates to the image. I know there are lots of ways to do this by matching up point sources with a known catalog, but in this case I'm generating a dust map, so point source matching won't work (as far as I can tell).

So the image is a 2D Numpy array of shape (240,240). It's written like so (the x and y coordinate assignments are a little weird, it works somehow):

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)

That all works fine on its own. For assigning galactic coordinates seems like all you'd need to do is something like:

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

But I don't understand how the FITS image stores these coordinates, so I don't know how to write them. I've tried assigning them in SAO DS9, too, but no luck. I just need a straightforward way of assigning these coordinates to the image.

Thanks for any help you can provide.

Answer 1

I would advise you to start using astropy. For the purposes of your project the astropy.wcs package can help you write a FITS WCS header, and the astropy.io.fits API is basically identical to the pyfits one you are using now. Moreover the help pages are excellent, and all I'm about to do is translate their WCS building page to match your example.

To your question: FITS does not "tag" each pixel with a coordinate. I suppose it is possible to create a pixel lookup table or something like that, but the actual WCS is an algorithmic translation of X,Y pixels to astrometric coordinates (in your case "Galactic"). A nice page is here.

The example I would point you to is here:

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

And here is my untested pseudocode for your project:

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