import matplotlib.pyplot as plt
import numpy as np
from astropy.stats import sigma_clipped_stats
from astropy.visualization import simple_norm
from photutils.aperture import CircularAperture, RectangularAperture
from photutils.datasets import (load_simulated_hst_star_image,
                                make_noise_image)
from photutils.detection import DAOStarFinder

hdu = load_simulated_hst_star_image()
data = hdu.data + make_noise_image(hdu.data.shape,
                                   distribution='gaussian',
                                   mean=10.0, stddev=5.0, seed=0)
mean, median, std = sigma_clipped_stats(data, sigma=3.0)
threshold = 5.0 * std
daofind = DAOStarFinder(threshold, fwhm=2.5, sharpness_range=(0.2, 1.5))
mask = np.zeros(data.shape, dtype=bool)
mask[650:851, 600:851] = True
mask[250:451, 150:551] = True
sources = daofind(data - median, mask=mask)
positions = np.transpose((sources['x_centroid'], sources['y_centroid']))
apertures = CircularAperture(positions, r=10.0)
fig, ax = plt.subplots()
norm = simple_norm(data, 'sqrt', percent=99)
axim = ax.imshow(data, norm=norm, origin='lower')
ax.set_title('Star finder with a mask to exclude regions')
p1 = apertures.plot(ax=ax, color='red')
rect1 = RectangularAperture((725, 750), 250, 200, theta=0)
rect2 = RectangularAperture((350, 350), 400, 200, theta=0)
p2 = rect1.plot(ax=ax, color='white', ls='dashed')
p3 = rect2.plot(ax=ax, color='white', ls='dashed')