import numpy as np
import matplotlib.pyplot as plt

from astropy.stats import gaussian_fwhm_to_sigma
from astropy.convolution import Gaussian2DKernel
from astropy.wcs import WCS
import astropy.units as u

from photutils import detect_threshold, source_properties, detect_sources
from photutils import SkyEllipticalAperture
from photutils import datasets

# get a mock image
hdu = datasets.load_star_image()
threshold = detect_threshold(hdu.data, snr=3.)
image_wcs = WCS(hdu.header)

# detect sources
sigma = 2.0 * gaussian_fwhm_to_sigma
kernel = Gaussian2DKernel(sigma, x_size=3, y_size=3)
kernel.normalize()
segm = detect_sources(hdu.data, threshold, npixels=5, filter_kernel=kernel)

# get properties of detected objects
props = source_properties(hdu.data, segm, wcs=image_wcs)


# ---
# create a sky aperture using detected properties

# the unit of "semimajor_axis_sigma" should be "pix"
print(props[0].semimajor_axis_sigma.unit)

# this works
aperture_1 = SkyEllipticalAperture(positions=props[0].sky_centroid,
        a=props[0].semimajor_axis_sigma.value * u.pix,
        b=props[0].semiminor_axis_sigma.value * u.pix,
        theta=props[0].orientation)

# but this doesn't
aperture_2 = SkyEllipticalAperture(positions=props[0].sky_centroid,
        a=props[0].semimajor_axis_sigma,
        b=props[0].semiminor_axis_sigma,
        theta=props[0].orientation)