# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Base class for profiles.
"""
import abc
import warnings
from functools import cached_property
import astropy.units as u
import numpy as np
from astropy.utils.exceptions import AstropyUserWarning
from photutils.aperture.core import _update_method_subpixels_docstring
from photutils.aperture.photometry import AperturePhotometry
from photutils.utils._deprecation import deprecated_positional_kwargs
from photutils.utils._quantity_helpers import process_quantities
from photutils.utils._stats import nanmax, nansum
__all__ = ['ProfileBase']
[docs]
@_update_method_subpixels_docstring
class ProfileBase(metaclass=abc.ABCMeta):
# numpydoc ignore: PR01,PR02,PR04,PR07
"""
Abstract base class for profile classes.
Parameters
----------
data : 2D `~numpy.ndarray`
The 2D data array. The data should be background-subtracted.
xycen : tuple of 2 floats
The ``(x, y)`` pixel coordinate of the source center.
radii : 1D float `~numpy.ndarray`
An array of radii defining the profile apertures. ``radii`` must
be strictly increasing with a minimum value greater than or
equal to zero, and contain at least 2 values. The radial spacing
does not need to be constant. See the subclass documentation for
details on how ``radii`` is interpreted.
error : 2D `~numpy.ndarray`, optional
The 1-sigma errors of the input ``data``. ``error`` is assumed
to include all sources of error, including the Poisson error of
the sources (see `~photutils.utils.calc_total_error`). ``error``
must have the same shape as the input ``data``.
mask : 2D bool `~numpy.ndarray`, optional
A boolean mask with the same shape as ``data`` where a `True`
value indicates the corresponding element of ``data`` is masked.
Masked data are excluded from all calculations.
<method_subpixels_descriptions>
"""
# Define axis labels used by `~photutils.profiles.ProfileBase.plot`.
# Subclasses may override these.
_xlabel = 'Radius (pixels)'
_ylabel = 'Profile'
# The user-facing name of the ``radii`` parameter, used in
# validation error messages. Subclasses that rename the parameter
# (e.g., ``half_sizes``) may override this.
_radii_name = 'radii'
def __init__(self, data, xycen, radii, *, error=None, mask=None,
method='exact', subpixels=5):
(data, error), unit = process_quantities((data, error),
('data', 'error'))
if error is not None and error.shape != data.shape:
msg = 'error must have the same shape as data'
raise ValueError(msg)
self.data = data
self.unit = unit
self.xycen = xycen
self.radii = self._validate_radii(radii)
self.error = error
self.mask = self._compute_mask(data, error, mask)
self.method = method
self.subpixels = subpixels
# The total normalization applied to the profile. This is the
# only mutable state of the class. All normalization-dependent
# attributes (e.g., ``profile``) are derived from immutable
# cached values divided by this value, and `normalize` and
# `unnormalize` update it with a single atomic attribute
# store, making the class safe for concurrent reads during
# normalization changes.
self.normalization_value = 1.0
def _validate_radii(self, radii):
"""
Validate and return the radii array.
"""
name = self._radii_name
if isinstance(radii, u.Quantity):
msg = (f'{name} must be a plain array of pixel values, '
'not a Quantity')
raise TypeError(msg)
radii = np.array(radii)
if radii.ndim != 1 or radii.size < 2:
msg = f'{name} must be a 1D array and have at least two values'
raise ValueError(msg)
if radii.min() < 0:
msg = f'minimum {name} must be >= 0'
raise ValueError(msg)
if not np.all(radii[1:] > radii[:-1]):
msg = f'{name} must be strictly increasing'
raise ValueError(msg)
return radii
def _compute_mask(self, data, error, mask):
"""
Compute the mask array, automatically masking non-finite data or
error values.
"""
badmask = ~np.isfinite(data)
if error is not None:
badmask |= ~np.isfinite(error)
if mask is not None:
if mask.shape != data.shape:
msg = 'mask must have the same shape as data'
raise ValueError(msg)
# Keep only non-finite values not already masked by the user
badmask &= ~mask
combined_mask = mask | badmask # all masked pixels
else:
combined_mask = badmask
if np.any(badmask):
msg = ('Input data contains non-finite values (e.g., NaN '
'or inf) that were automatically masked.')
warnings.warn(msg, AstropyUserWarning)
return combined_mask
@property
@abc.abstractmethod
def radius(self):
"""
The profile radius in pixels as a 1D `~numpy.ndarray`.
"""
@property
@abc.abstractmethod
def _raw_profile(self):
"""
The raw (unnormalized) profile as a 1D `~numpy.ndarray`.
"""
@property
@abc.abstractmethod
def _raw_profile_error(self):
"""
The raw (unnormalized) profile errors as a 1D `~numpy.ndarray`.
If no ``error`` array was provided, an empty array with shape
``(0,)`` is returned.
"""
@property
def profile(self):
"""
The profile as a 1D `~numpy.ndarray`.
The returned values reflect the current profile normalization
(see `normalize`).
"""
return self._raw_profile / self.normalization_value
@property
def profile_error(self):
"""
The profile errors as a 1D `~numpy.ndarray`.
The returned values reflect the current profile normalization
(see `normalize`). If no ``error`` array was provided, an empty
array with shape ``(0,)`` is returned.
"""
return self._raw_profile_error / self.normalization_value
@cached_property
def _circular_apertures(self):
"""
A list of `~photutils.aperture.CircularAperture` objects.
The first element may be `None`.
"""
from photutils.aperture import CircularAperture
apertures = []
for radius in self.radii:
if radius <= 0.0:
apertures.append(None)
else:
apertures.append(CircularAperture(self.xycen, radius))
return apertures
def _compute_photometry(self, apertures):
"""
Compute aperture fluxes, flux errors, and areas for the given
apertures using a single batched photometry call.
Parameters
----------
apertures : list
A list of aperture objects sharing the same position.
Leading elements may be `None` (e.g., for a zero radius), in
which case the corresponding flux, error, and area are set
to zero.
Returns
-------
flux : `~numpy.ndarray`
The aperture fluxes.
flux_err : `~numpy.ndarray`
The aperture flux errors.
areas : `~numpy.ndarray`
The aperture areas.
"""
n_none = sum(aperture is None for aperture in apertures)
result = AperturePhotometry(
self.data, apertures[n_none:], error=self.error,
mask=self.mask, method=self.method, subpixels=self.subpixels)
fluxes = result.flux
areas = result.area.to_value(u.pix ** 2)
flux_errs = (result.flux_err if self.error is not None
else np.array([]))
if n_none > 0:
zeros = np.zeros(n_none)
fluxes = np.concatenate((zeros, fluxes))
areas = np.concatenate((zeros, areas))
if self.error is not None:
flux_errs = np.concatenate((zeros, flux_errs))
if self.unit is not None:
fluxes = fluxes << self.unit
flux_errs = flux_errs << self.unit
return fluxes, flux_errs, areas
@cached_property
def _photometry(self):
"""
The aperture fluxes, flux errors, and areas as a function of
radius.
"""
return self._compute_photometry(self._circular_apertures)
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
def normalize(self, method='max'):
"""
Normalize the profile.
The normalization is computed from the raw (unnormalized)
profile values, so repeated calls do not accumulate. The
most recent call determines the normalization. Because both
normalization methods scale linearly with the profile values,
this is equivalent to normalizing an already-normalized profile.
Parameters
----------
method : {'max', 'sum'}, optional
The method used to normalize the profile:
* ``'max'`` (default):
The profile is normalized such that its maximum value is
1.
* ``'sum'``:
The profile is normalized such that the sum of its values
is 1.
"""
if method == 'max':
func = nanmax
elif method == 'sum':
func = nansum
else:
msg = "invalid method, must be 'max' or 'sum'"
raise ValueError(msg)
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
normalization = func(self._raw_profile)
if normalization == 0 or not np.isfinite(normalization):
msg = ('The profile cannot be normalized because the max or '
'sum is zero or non-finite.')
warnings.warn(msg, AstropyUserWarning)
else:
# A single atomic attribute store. Concurrent readers see
# either the old or the new normalization, never a mixed
# state.
self.normalization_value = normalization
[docs]
def unnormalize(self):
"""
Unnormalize the profile back to the original state before any
calls to `normalize`.
"""
self.normalization_value = 1.0
@staticmethod
def _trim_to_monotonic(xarr, profile, name):
"""
Trim arrays to the first monotonically increasing region.
This is used by interpolation methods that require a
monotonically increasing profile.
Parameters
----------
xarr : 1D `~numpy.ndarray`
The x-axis values (e.g., radius or half-size).
profile : 1D `~numpy.ndarray`
The profile values.
name : str
A descriptive name for the profile used in the error
message.
Returns
-------
xarr, profile : tuple of `~numpy.ndarray`
The trimmed arrays.
"""
finite_mask = np.isfinite(profile)
if not np.all(finite_mask):
# Keep only the leading finite segment
first_nonfinite = np.argmin(finite_mask)
xarr = xarr[:first_nonfinite]
profile = profile[:first_nonfinite]
# np.diff produces an array of length n-1: diff[i] represents
# the step from profile[i] to profile[i+1]. A value <= 0 means
# the profile stopped increasing at that step.
diff = np.diff(profile) <= 0
if np.any(diff):
# idx is an index into the *diff* array, not the profile
# array. diff[idx] <= 0 means the drop occurs between
# profile[idx] and profile[idx+1], so profile[idx] is
# the last good value. We therefore need profile[:idx+1]
# (inclusive) to retain it.
idx = np.argmax(diff) # first non-monotonic step in diff-space
xarr = xarr[:idx + 1]
profile = profile[:idx + 1]
if len(xarr) < 2:
msg = (f'The {name} profile is not monotonically '
'increasing even at the smallest values -- cannot '
'interpolate. Try using different input values '
'(especially the starting values) and/or using the '
'"exact" aperture overlap method.')
raise ValueError(msg)
return xarr, profile
def __repr__(self):
cls_name = self.__class__.__name__
n_radii = len(self.radii)
normalized = self.normalization_value != 1.0
return (f'{cls_name}(xycen={self.xycen}, n_radii={n_radii}, '
f'normalized={normalized})')
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
def plot(self, ax=None, **kwargs):
"""
Plot the profile.
Parameters
----------
ax : `matplotlib.axes.Axes` or `None`, optional
The matplotlib axes on which to plot. If `None`, then the
current `~matplotlib.axes.Axes` instance is used.
**kwargs : dict, optional
Any keyword arguments accepted by `matplotlib.pyplot.plot`.
Returns
-------
lines : list of `~matplotlib.lines.Line2D`
A list of lines representing the plotted data.
"""
import matplotlib.pyplot as plt
if ax is None:
ax = plt.gca()
profile = self.profile
unit = None
if isinstance(profile, u.Quantity):
unit = profile.unit
profile = profile.value
lines = ax.plot(self.radius, profile, **kwargs)
ax.set_xlabel(self._xlabel)
ylabel = self._ylabel
# A normalized profile is dimensionless, so the unit is
# included only when the profile has a physical unit
if unit is not None and unit != u.dimensionless_unscaled:
ylabel = f'{ylabel} ({unit})'
ax.set_ylabel(ylabel)
return lines
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
def plot_error(self, ax=None, **kwargs):
"""
Plot the profile errors.
Parameters
----------
ax : `matplotlib.axes.Axes` or `None`, optional
The matplotlib axes on which to plot. If `None`, then the
current `~matplotlib.axes.Axes` instance is used.
**kwargs : dict, optional
Any keyword arguments accepted by
`matplotlib.pyplot.fill_between`.
Returns
-------
poly : `matplotlib.collections.PolyCollection` or `None`
A `~matplotlib.collections.PolyCollection` containing the
plotted polygons, or `None` if no errors were input.
"""
if len(self.profile_error) == 0:
msg = 'Errors were not input'
warnings.warn(msg, AstropyUserWarning)
return None
import matplotlib.pyplot as plt
if ax is None:
ax = plt.gca()
# Set default fill_between facecolor.
# facecolor must be first key, otherwise it will override color
# kwarg (i.e., cannot use setdefault here)
if 'facecolor' not in kwargs:
kws = {'facecolor': (0.5, 0.5, 0.5, 0.3)}
kws.update(kwargs)
else:
kws = kwargs
profile = self.profile
profile_error = self.profile_error
if isinstance(profile, u.Quantity):
profile = profile.value
if isinstance(profile_error, u.Quantity):
profile_error = profile_error.value
ymin = profile - profile_error
ymax = profile + profile_error
return ax.fill_between(self.radius, ymin, ymax, **kws)
