4

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 actual result But I'm trying to get something more like expected result where the blue square is the actual position of the image in image_file

usernumber
  • 1,958
  • 1
  • 21
  • 58

1 Answers1

4

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

plot of horsehead nebula

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.

Matt Pitkin
  • 3,989
  • 1
  • 18
  • 32
  • But how do I extend the grid beyond the frame of the image? – usernumber Feb 15 '19 at 12:42
  • 1
    You can use `ax.set_xlim()` and `ax.set_ylim()` although I'm not sure what units the limits are in (potentially something like arcseconds). So, if you do `ax.set_xlim((-180*60, 180*60))` and `ax.set_ylim((-90*60, 90*60))` you'll start to see more (but not the whole) of the sky. Note that this image actually covers a very small bit of the sky, so if you have the whole sky it will just look like a dot. – Matt Pitkin Feb 15 '19 at 13:06
  • But that leaves the image centered on the horsehead nebula rather than on the coordinate (0,0) point – usernumber Feb 15 '19 at 13:25
  • I'm afraid it's hit the limit of my knowledge. You could also have a look at trying [aplpy](https://aplpy.readthedocs.io/en/stable/), although the documentation is a bit limited and I couldn't get it to work. Or it might be worth asking the question on the astropy forum http://www.astropy.org/contribute.html – Matt Pitkin Feb 15 '19 at 13:53