title: How do I overlay two .fits images each with different WCS / resolutions?
body: <p>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. </p>
  
  <p>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. </p>
  
  <pre><code>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')
  </code></pre>
  
  <p><a href="https://i.stack.imgur.com/WNiMl.png" rel="noreferrer">Resulting plot (click me)</a></p>
  
  <p>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:</p>
  
  <pre><code>array, footprint = reproject_interp((stamp_a.data, wcs_a), image_z.header)
  plt.imshow(array, origin='lower')
  </code></pre>
  
  <p>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.</p>