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.

def add_plot(fn: Path, cm, ax, alpha=1):
    """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.
    """

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