# Get the data
from astropy.io import fits
hdulist = fits.open('https://lambda.gsfc.nasa.gov/data/map/dr3/skymaps/5yr//wmap_band_imap_r9_5yr_K_v3.fits')

# Set up the HEALPix projection
from astropy_healpix import HEALPix
from astropy.coordinates import Galactic
nside = hdulist[1].header['NSIDE']
order = hdulist[1].header['ORDERING']
hp = HEALPix(nside=nside, order=order, frame=Galactic())

# Sample a 300x200 grid in RA/Dec
from astropy import units as u
ra = np.linspace(-15., 15., 300) * u.deg
dec = np.linspace(-10., 10., 200) * u.deg
ra_grid, dec_grid = np.meshgrid(ra, dec)

# Set up Astropy coordinate objects
from astropy.coordinates import SkyCoord
coords = SkyCoord(ra_grid.ravel(), dec_grid.ravel(), frame='icrs')

# Interpolate values
temperature = hdulist[1].data['temperature']
tmap = hp.interpolate_bilinear_skycoord(coords, temperature)
tmap = tmap.reshape((200, 300))

# Make a plot of the interpolated temperatures
plt.figure(figsize=(9, 5))
im = plt.imshow(tmap, extent=[-1, 1, -10, 10], cmap=plt.cm.RdYlBu, origin='lower', aspect='auto')
plt.colorbar(im)
plt.xlabel('Right ascension (ICRS)')
plt.ylabel('Declination (ICRS)')
plt.show()