from astropy.visualization import (ImageNormalize,
                                   AsymmetricPercentileInterval)
from astropy.wcs import WCS
from matplotlib import pyplot as plt
from ligo.skymap.plot import mellinger
from reproject import reproject_interp

ax = plt.axes(projection='astro hours aitoff')
backdrop = mellinger()
backdrop_wcs = WCS(backdrop.header).dropaxis(-1)
interval = AsymmetricPercentileInterval(45, 98)
norm = ImageNormalize(backdrop.data, interval)
backdrop_reprojected = np.asarray([
    reproject_interp((layer, backdrop_wcs), ax.header)[0]
    for layer in norm(backdrop.data)])
backdrop_reprojected = np.rollaxis(backdrop_reprojected, 0, 3)
ax.imshow(backdrop_reprojected)