import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
from photutils.aperture import EllipticalAperture
from photutils.datasets import make_4gaussians_image
from photutils.morphology import data_properties

slc = np.s_[40:80, 75:105]
data = make_4gaussians_image()[slc]  # extract single object
mask = data < 50
cat = data_properties(data, mask=mask)
columns = ['label', 'x_centroid', 'y_centroid', 'semimajor_axis',
           'semiminor_axis', 'orientation']
tbl = cat.to_table(columns=columns)
r = 2.5  # approximate isophotal extent
xypos = (cat.x_centroid, cat.y_centroid)
a = cat.semimajor_axis.value * r
b = cat.semiminor_axis.value * r
theta = cat.orientation.to_value(u.radian)
apertures = EllipticalAperture(xypos, a, b, theta=theta)

fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(data, origin='lower')
apertures.plot(ax=ax, color='C3', lw=2)

dx_major = a * np.cos(theta)
dy_major = a * np.sin(theta)
color = 'C1'
width = 0.2
ax.arrow(cat.x_centroid, cat.y_centroid, dx_major, dy_major, color=color,
         length_includes_head=True, width=width)
theta2 = theta + np.pi / 2
dx_minor = b * np.cos(theta2)
dy_minor = b * np.sin(theta2)
ax.arrow(cat.x_centroid, cat.y_centroid, dx_minor, dy_minor, color=color,
         length_includes_head=True, width=width)