code: def add_plot(fn: Path, cm, ax, alpha=1):
      
  
      with fits.open(fn, mode='readonly', memmap=False) as f:
          img = f[0].data
  
          yPix, xPix = f[0].shape[-2:]
          x, y = np.meshgrid(range(xPix), range(yPix))  # pixel indices to find RA/dec of
          xy = np.column_stack((x.ravel(order='C'), y.ravel(order='C')))
  
          radec = wcs.WCS(f[0].header).all_pix2world(xy, 0)
  
      ra = radec[:, 0].reshape((yPix, xPix), order='C')
      dec = radec[:, 1].reshape((yPix, xPix), order='C')
  
      ax.set_title(fn.name)
      ax.pcolormesh(ra, dec, img, alpha=alpha, cmap=cm, norm=LogNorm())
      ax.set_ylabel('Right Ascension [deg.]')
      ax.set_xlabel('Declination [deg.]')
docs: Astrometry.net makes file ".new" with the image and the WCS SIP 2-D polynomial fit coefficients in the FITS header
  
      We use DECL as "x" and RA as "y".
      pcolormesh() is used as it handles arbitrary pixel shapes.
      Note that pcolormesh() cannot tolerate NaN in X or Y (NaN in C is OK).
      This is handled in https://github.com/scivision/pcolormesh_nan.py.