Writing a FrameSet to a PyFITS header

[rmjarvis]

I'm trying to figure out how to take a WCS that I have in a FrameSet (in python) and write it into a FITS header. Following the documentation in the help for PyFITSAdapter, I came up with the following:

import starlink.Ast
import starlink.Atl
import pyfits

f = pyfits.open('1904-66_TAN.fits')
# File from http://www.atnf.csiro.au/people/mcalabre/WCS/example_data.html

fc = starlink.Ast.FitsChan( starlink.Atl.PyFITSAdapter(f[0]) )
wcs = fc.read()
f.close()

hdu = pyfits.PrimaryHDU()
fc2 = starlink.Ast.FitsChan( None, starlink.Atl.PyFITSAdapter(hdu) )
fc2.write(wcs)
fc2.writefits()

print hdu.header

But at the end of this, the header is empty.

Before the writefits() command, the header is still its original value with keys like SIMPLE, BITPIX, NAXIS. But after writefits(), it is empty. So that function is clearly doing something. Just not what I expected.

Could you please explain what I am doing wrong here? Thanks so much for your help!

This was caused by a bug in the PyFITSAdapter class. Also note that the "fc2" FitsChan will use "NATIVE" encoding by default (i.e. it will use an AST-specific keywords for the WCS information). If you want standard FITS-WCS keywords, you need to append "Encoding=FITS-WCS" to the end of the FitsChan constructor call.

[dsberry]

This should be fixed by commit 4b4cbacf94. It was a bug within the Atl.PyFITSAdapter class that caused blank keywords to be handled incorrectly. Version 2.6 of pyast includes this fix and is now available on pypi.

[rmjarvis]

Thanks for addressing this so quickly. Unfortunately however, this doesn't quite fix the problem for me.

It does copy a bunch of the header information from the original. But it doesn't get everything. And in particular, it doesn't write the WCS keys, which are the only ones I actually care about. Here is a modified script for you to try:

import starlink.Ast
import starlink.Atl
import pyfits

f = pyfits.open('1904-66_TAN.fits')
h1 = f[0].header
fc = starlink.Ast.FitsChan( starlink.Atl.PyFITSAdapter(f[0]) )
wcs = fc.read()
f.close()

hdu = pyfits.PrimaryHDU()
fc2 = starlink.Ast.FitsChan( None, starlink.Atl.PyFITSAdapter(hdu) )
fc2.write(wcs)
fc2.writefits()
h2 = hdu.header

for key in h1:
    if key not in h2:
        print 'h2 is missing key ',key

Among others, I find that it is missing CTYPE, CRPIX, CRVAL*, etc.

[rmjarvis]

Ah, never mind. I was misinterpreting the output. I see now that it is writing its own keys, not copying the original keys.

Is there any way to get it to convert its version back to standard fits conventions that can be read by ds9, etc.? Or is that just functionality that starlink doesn't have?

[dsberry]

The easiest way to do this is probably just to modify the "fc2 = ..." line to be

fc2 = starlink.Ast.FitsChan( None, starlink.Atl.PyFITSAdapter(hdu), "Encoding=FITS-WCS" )

Funnily enough, Bill Joye is planning that the next version of DS9 will be able to read the AST native encoding.

If you want the "fc2.write(wcs)" operation to create standard FITS-WCS keywords such as CTYPE, CRVAL, etc, then you need to indicate that "fc2" should use FITS-WCS encoding. An empty FitsChan will use AST native encoding by default (a keyword scheme that overcomes many of the severe restrictions imposed by FITS-WCS). See the following links for more information about encodings:

http://starlink.jach.hawaii.edu/devdocs/sun211.htx/node170.html http://starlink.jach.hawaii.edu/devdocs/sun211.htx/node155.html http://starlink.jach.hawaii.edu/devdocs/sun211.htx/node495.html

[rmjarvis]

Ah, perfect! That's exactly what I was looking for. Thanks so much!