Plot a FITS image over a grid of the entire sky

Question

I would like to plot a FITS image over a sketch of the entire sky.

I got this far in plotting the entire sky:

import matplotlib.pyplot as plt
from astropy.utils.data import get_pkg_data_filename
from astropy.io import fits

image_file = get_pkg_data_filename('tutorials/FITS-images/HorseHead.fits')
image_data = fits.getdata(image_file, ext=0)

plt.subplot(111, projection='aitoff')
plt.grid(True)
plt.imshow(image_data, cmap='gray')
plt.show()

But I can't seem to get the FITS image properly aligned with the grid

---

The above code results in the following But I'm trying to get something more like where the blue square is the actual position of the image in image_file

Answer 1

The image_data is just an array of pixel values, and doesn't actually contain any information about the positions and orientations of the pixels on a sky map. You have to extract that information from the FITS file too.

An example would be (based in information from here and here):

import matplotlib.pyplot as plt
from astropy.utils.data import get_pkg_data_filename
from astropy.io import fits
from astropy.wcs import WCS
from astropy.visualization.wcsaxes.frame import EllipticalFrame

image_file = get_pkg_data_filename('tutorials/FITS-images/HorseHead.fits')
image_header = fits.open(image_file)[0].header  # extract header info
image_data = fits.getdata(image_file, ext=0)

# use the WCS class to get coordinate info and projection axes to use
wcs = WCS(image_header)

# make the axes using the EllipticalFrame class
ax = plt.subplot(projection=wcs, frame_class=EllipticalFrame)

# add the image
im = ax.imshow(image_data, origin='lower')

# add a grid
overlay = ax.get_coords_overlay('fk5')
overlay.grid(color='white', ls='dotted')

As the particular image example doesn't extend across the whole sky you'll only see a small patch (but you can extend the plot using ax.set_xlim() and ax.set_ylim() as required, although I'm not sure what the units are, and it's also worth noting that this image actually just covers a very small patch of sky), but if your image was over the whole sky it would look like an Aitoff projection as shown here.

I think this only works with the latest version of astropy (v3.1), Matplotlib (v3.0.2), and therefore requires Python >= 3.5.

You may also want to look at healpy, although that can't read in the astropy example FITS file.