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
from photutils.datasets import (load_simulated_hst_star_image,
                                make_noise_image)
from photutils.detection import find_peaks

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 = median + (5.0 * std)
tbl = find_peaks(data, threshold, box_size=11)
positions = np.transpose((tbl['x_peak'], tbl['y_peak']))
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')
patches = apertures.plot(color='red')