Writing axes to FITS file

Question

Is there a way to write existing xy axes to a FITS file along with the data itself in Python?

For example here is some simple code saving a matrix to a FITS file named TestFITS:

import numpy as np
from astropy.io import fits

test_matrix = np.random.uniform(0,1,[5,3])
x = np.arange(5,5+len(test_matrix[:,0]))
y = np.arange(5,5+len(test_matrix[0,:]))
hdu = fits.PrimaryHDU(test_matrix)
hdu.writeto('TestFITS')

But if I wished to save x and y to the file as well could that be done?

Answer 1

You could save them as one-dimensional ImageHDUs in two extensions, next to the PrimaryHDU:

import numpy as np
from astropy.io import fits

test_matrix = np.random.uniform(0,1,[5,3])
x = np.arange(5,5+len(test_matrix[:,0]))
y = np.arange(5,5+len(test_matrix[0,:]))
fits.HDUList([
    fits.PrimaryHDU(test_matrix),
    fits.ImageHDU(x, name='X'),
    fits.ImageHDU(y, name='Y'),
    ]).writeto('testxy.fits')

(The name parameter is not necessary, but can be a nice convenience.)