Overlay the MOC of Spitzer on a 2MASS cutout

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from astropy.visualization import simple_norm
from astropy.wcs import WCS
from mocpy import MOC

# load 2MASS cutout covering the galactic plane
hdu = fits.open(
    "http://alasky.u-strasbg.fr/hips-image-services/hips2fits?hips=CDS%2FP%2F2MASS%2FK&width=1200&height=700&fov=30&projection=TAN&coordsys=galactic&rotation_angle=0.0&object=gal%20center&format=fits",
)

# load Spitzer MOC
moc = MOC.from_fits("http://skies.esac.esa.int/Spitzer/IRAC1_bright_ISM/Moc.fits")

# create WCS from 2MASS image header
twomass_wcs = WCS(header=hdu[0].header)

fig = plt.figure(111, figsize=(15, 10))
ax = fig.add_subplot(1, 1, 1, projection=twomass_wcs)
moc.fill(ax=ax, wcs=twomass_wcs, alpha=0.3, color="red")
im = ax.imshow(
    hdu[0].data,
    origin="lower",
    norm=simple_norm(hdu[0].data, "sqrt", min_cut=0, max_cut=150),
)

# compute skycoords for every pixel of the image
width = hdu[0].header["NAXIS1"]
height = hdu[0].header["NAXIS2"]

xv, yv = np.meshgrid(np.arange(0, width), np.arange(0, height))
skycoords = twomass_wcs.pixel_to_world(xv, yv)

ra, dec = skycoords.icrs.ra.deg, skycoords.icrs.dec.deg

mask_in_moc = moc.contains(ra * u.deg, dec * u.deg)

/tmp/ipykernel_326933/2840231297.py:1: DeprecationWarning: This method is deprecated and has been replaced by contains_lonlat
  mask_in_moc = moc.contains(ra * u.deg, dec * u.deg)

img = hdu[0].data
img_test = img.copy()

img_test[~mask_in_moc] = 0

fig = plt.figure(111, figsize=(15, 10))
ax = fig.add_subplot(1, 1, 1, projection=twomass_wcs)

im = ax.imshow(
    img_test,
    origin="lower",
    norm=simple_norm(hdu[0].data, "sqrt", min_cut=0, max_cut=150),
)

# number of pixels in the mask
np.sum(img[mask_in_moc])

8690767.0

# number of pixels in the entire image
np.sum(img)

14976694.0