How do I overlay two .fits images each with different WCS / resolutions?

Question

I'm trying to plot contours from one image of an object, let's say in band A, on top of a lower resolution image of that same object, let's say in band Z.

Both images are in large fits files, so I have to create a 2D cutout of each. However, the image in band Z has far fewer pixels than that of band A. So, when I try to plot them in the same figure, the image from band Z plots according to pixel dimensions and is therefore very tiny in the bottom left corner of the figure (see plot in link below). I need both images to plot on the basis of astronomical angular dimensions, not pixels.

import numpy as np
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
from astropy.nddata import Cutout2D
from reproject import reproject_interp

# Read in .fits images and WCS information
image_a = fits.open("image_a.fits")
image_z = fits.open("image_z.fits")
wcs_a = WCS(image_a.header).celestial
wcs_z = WCS(image_z.header).celestial

# Define object RA and DEC
obj_coords = SkyCoord(ra=135.19081*u.degree, dec=0.68991393*u.degree, frame='fk5')

# Create 2D cutouts of the object in each band in a 6 by 6 arcsec box
size = u.Quantity((6, 6), u.arcsec)
stamp_a = Cutout2D(image_a.data, obj_coords, size, wcs=wcs_a)
stamp_z = Cutout2D(image_z.data, obj_coords, size, wcs=wcs_z)

# Plot the the stamp in band Z with contours from band A
ax = plt.subplot(projection = wcs_z)
plt.imshow(stamp_z.data, cmap='gray', origin='lower', alpha=1)
ax.contour(stamp_a.data, levels = [2*sigma,3*sigma,4*sigma,5*sigma], 
colors='green')

Resulting plot (click me)

I've also seen something about using reproject. I tried it but don't think I totally understand it as it just results in a tiny little pixel. This is what I tried:

array, footprint = reproject_interp((stamp_a.data, wcs_a), image_z.header)
plt.imshow(array, origin='lower')

Let me know if there is any info I'm missing to answer the question. I've just made the switch to python last week after using IDL for a few years, so I apologize if it's a bit rough.

Answer 1

The problem with how you use reproject is that you pass (stamp_a.data, wcs_a), but wcs_a is the WCS from the original image, not from the stamp.

You can get a WCS object that matches your stamp from the Cutout2D image. I think changing to (stamp_a.data, stamp_a.wcs) will give you a correct result.

Have you had a look at astropy.visualisation.wcsaxes yet? You can see an example here how to plot an image, and overplot contours from a second image with a different WCS. If this works for you, it will be simpler than using Cutout2D and reproject manually.

I'm not sure concerning performance; if your images are huge you still might want to do the cutouts like you do now. If memory / execution time isn't a problem, you can just plot the whole image and then set the range you want with wcsaxes (see their tutorial docs).

The question whether to do the reproject yourself or let matplotlib and astropy.visualisation.wcsaxes behind the scenes is also a matter of taste. Using reproject directly is a bit more code and more complex, but will give you more control over the exact reprojection method used (e.g. interpolation or the slower flux-conserving reproject_exact) and might be easier to fully understand what is going on.