Plotting lines over an image using the same projection

Question

I want to make a plot using .fits files (astronomical images) and I am experiencing two issues which I think they are related:

Using this example from astropy:

from matplotlib import pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
from astropy.utils.data import download_file

fits_file = 'http://data.astropy.org/tutorials/FITS-images/HorseHead.fits'
image_file = download_file(fits_file, cache=True)
hdu = fits.open(image_file)[0]
wcs = WCS(hdu.header)

fig = plt.figure()
fig.add_subplot(111, projection=wcs)
plt.imshow(hdu.data, origin='lower', cmap='cubehelix')
plt.xlabel('RA')
plt.ylabel('Dec')
plt.show()

I can generate this image:

Now I would like to plot some points using the same coordinates as the image:

plt.scatter(85, -2, color='red')

However, when I do this:

I am ploting at the pixel coordinantes. Furthermore, the image no longer matches the frame size (although the coordinates seem fine)

Any advice on how to deal with these issues?

Answer 1

It is very easy to plot given coordinates. All you have to do is apply a transform.

I copied your example and added comments where I changed something and why.

from matplotlib import pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
from astropy.utils.data import download_file

fits_file = 'http://data.astropy.org/tutorials/FITS-images/HorseHead.fits'
image_file = download_file(fits_file, cache=True)

# Note that it's better to open the file with a context manager so no
# file handle is accidentally left open.
with fits.open(image_file) as hdus:
    img = hdus[0].data
    wcs = WCS(hdus[0].header)

fig = plt.figure()

# You need to "catch" the axes here so you have access to the transform-function.
ax = fig.add_subplot(111, projection=wcs)
plt.imshow(img, origin='lower', cmap='cubehelix')
plt.xlabel('RA')
plt.ylabel('Dec')

# Apply a transform-function:
plt.scatter(85, -2, color='red', transform=ax.get_transform('world'))

And the result is:

Note that if you want the Canvas to only show the region of the image just apply the limits again afterwards:

# Add a scatter point which is in the extend of the image:
plt.scatter(85.3, -2.5, color='red', transform=ax.get_transform('world'))

plt.ylim(0, img.shape[0])
plt.xlim(0, img.shape[1])

which gives:

A side note as well here. AstroPy has a great units support so instead of converting arcmins and arcsecs to degrees you can just define the "unit". You still need the transform though:

from astropy import units as u
x0 = 85 * u.degree + 20 * u.arcmin
y0 = -(2 * u.degree + 25 * u.arcmin)
plt.scatter(x0, y0, color='red', transform=ax.get_transform('world'))