import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
from regions.core import PixCoord
from regions.shapes.circle import CirclePixelRegion

filename = get_pkg_data_filename('photometry/M6707HH.fits')
hdulist = fits.open(filename)
hdu = hdulist[0]
center = PixCoord(158.5, 1053.5)
aperture = CirclePixelRegion(center, 4.)
mask = aperture.to_mask(mode='exact')
data = mask.cutout(hdu.data)
weighted_data = mask.multiply(hdu.data)

fig, axes = plt.subplots(nrows=1, ncols=3)
axes[0].set_title('Mask', size=9)
axes[0].imshow(mask.data, origin='lower',
               extent=mask.bbox.extent)
axes[1].set_title('Data cutout', size=9)
axes[1].imshow(data, origin='lower',
               extent=mask.bbox.extent)
axes[2].set_title('Data cutout multiplied by mask', size=9)
axes[2].imshow(weighted_data, origin='lower',
               extent=mask.bbox.extent)
hdulist.close()