# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Tools for calculating properties of sources defined by an Aperture.
"""
import inspect
import warnings
from copy import deepcopy
import astropy.units as u
import numpy as np
from astropy.coordinates import SkyCoord
from astropy.nddata import NDData, StdDevUncertainty
from astropy.stats import (SigmaClip, biweight_location, biweight_midvariance,
mad_std)
from astropy.utils import lazyproperty
from astropy.utils.exceptions import AstropyUserWarning
from photutils.aperture import Aperture, SkyAperture, region_to_aperture
from photutils.aperture._batch_photometry import (FLAG_COL_BBOX_CLIPPED,
FLAG_COL_MASKED,
FLAG_COL_NONFINITE_DATA,
FLAG_COL_NONFINITE_ERROR,
FLAG_COL_NPIX, FLAG_COL_SEG,
FLAG_COL_UNCORRECTED,
FLAG_COL_VALID,
batch_aperture_sums)
from photutils.aperture._batch_stats import (batch_aperture_gather,
batch_biweight, batch_gini,
batch_mad, batch_mean_var,
batch_moments, batch_order_stats,
batch_sigma_clip_center,
batch_sigma_clip_sum,
batch_sort_values)
from photutils.aperture._segmentation import (SEG_METHOD_CODES,
make_segmentation_exclusion,
process_segmentation_inputs)
from photutils.aperture.core import (_aperture_metadata,
_update_method_subpixels_docstring)
from photutils.aperture.flags import (APERTURE_FLAGS, _counts_to_flag_bits,
decode_aperture_flags)
from photutils.morphology import gini as gini_func
from photutils.utils._deprecation import (create_empty_deprecated_qtable,
deprecated, deprecated_getattr,
deprecated_positional_kwargs)
from photutils.utils._misc import _get_meta
from photutils.utils._moments import _image_moments
from photutils.utils._parameters import validate_table_columns
from photutils.utils._quantity_helpers import process_quantities
__all__ = ['ApertureStats']
# Scale factor that converts the median absolute deviation to a robust
# estimate of the standard deviation (1 / scipy.stats.norm.ppf(0.75)).
# This must match ``astropy.stats.mad_std`` and the value in
# ``photutils.aperture._batch_stats``.
_MAD_STD_SCALE = 1.482602218505602
# Remove in 4.0
_DEPRECATED_ATTRIBUTES: dict = {
'covar_sigx2': 'covariance_xx',
'covar_sigxy': 'covariance_xy',
'covar_sigy2': 'covariance_yy',
'cxx': 'ellipse_cxx',
'cxy': 'ellipse_cxy',
'cyy': 'ellipse_cyy',
'data_sumcutout': 'data_sum_cutout',
'error_sumcutout': 'error_sum_cutout',
'get_id': 'select_id',
'get_ids': 'select_ids',
'semimajor_sigma': 'semimajor_axis',
'semiminor_sigma': 'semiminor_axis',
'xcentroid': 'x_centroid',
'ycentroid': 'y_centroid',
}
# Public attributes that are never collapsed to a scalar for a scalar
# instance because they describe the whole object rather than a single
# per-source value (see ``ApertureStats.__getattribute__``).
_SCALAR_EXCLUDE = frozenset({'default_columns', 'isscalar', 'n_apertures',
'properties'})
class _UncachedLazyProperty(lazyproperty):
"""
A property that is discovered as a source property (like
`~astropy.utils.lazyproperty`) but is recomputed on every access
instead of being cached.
This is used for `ApertureStats.flags`, whose value can change
over the object's lifetime. It reflects whether covariance-derived
properties have been computed (see `ApertureStats.flags`). Caching
would freeze the first-computed value, so the getter runs on every
access. Subclassing `~astropy.utils.lazyproperty` keeps it in the
`ApertureStats.properties` list and the ``to_table`` machinery.
"""
def __get__(self, obj, owner=None):
if obj is None:
return self
return self.fget(obj)
[docs]
@_update_method_subpixels_docstring
class ApertureStats:
# numpydoc ignore: PR01,PR02,PR04,PR07
"""
Class to create a catalog of statistics for pixels within an
aperture.
Note that this class returns the statistics of the input
``data`` values within the aperture. It does not convert data
in surface brightness units to flux or counts. Conversion from
surface-brightness units should be performed before using this
function.
Parameters
----------
data : 2D `~numpy.ndarray`, `~astropy.units.Quantity`, \
`~astropy.nddata.NDData`
The 2D array from which to calculate the source properties.
For accurate source properties, ``data`` should be
background-subtracted. Non-finite ``data`` values (NaN and inf)
are automatically masked.
aperture : `~photutils.aperture.Aperture` or supported `~regions.Region`
The aperture or region to apply to the data. The aperture
or region object may contain more than one position. If the
input ``aperture`` is a `~photutils.aperture.SkyAperture` or
`~regions.SkyRegion` object, then a WCS must be input using
the ``wcs`` keyword. Region objects are converted to aperture
objects.
error : 2D `~numpy.ndarray` or `~astropy.units.Quantity`, optional
The total error array corresponding to the input ``data``
array. ``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``. If ``data`` is a
`~astropy.units.Quantity` array then ``error`` must be a
`~astropy.units.Quantity` array (and vice versa) with identical
units. Non-finite ``error`` values (NaN and +/- inf) are not
automatically masked, unless they are at the same position of
non-finite values in the input ``data`` array. Such pixels can
be masked using the ``mask`` keyword.
mask : 2D `~numpy.ndarray` (bool), 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. Non-finite
values (NaN and inf) in the input ``data`` are automatically
masked.
wcs : WCS object or `None`, optional
A world coordinate system (WCS) transformation that
supports the `astropy shared interface for WCS
<https://docs.astropy.org/en/stable/wcs/wcsapi.html>`_ (e.g.,
`astropy.wcs.WCS`, `gwcs.wcs.WCS`). ``wcs`` is required if
the input ``aperture`` is a `~photutils.aperture.SkyAperture`
or `~regions.SkyRegion` object. If `None`, then all sky-based
properties will be set to `None`.
sigma_clip : `None` or `astropy.stats.SigmaClip` instance, optional
A `~astropy.stats.SigmaClip` object that defines the sigma
clipping parameters. If `None` then no sigma clipping will
be performed.
sum_method : {'exact', 'center', 'subpixel'}, optional
The method used to determine the pixel weights (the fraction
of the pixel area covered by the aperture). This method
is used only for calculating the ``sum``, ``sum_error``,
``sum_aper_area``, ``data_sum_cutout``, and ``error_sum_cutout``
properties. All other properties use the "center" aperture mask
method. The following methods are available:
<method_bullets>
<subpixels_description>
ddof : int, optional
The delta degrees of freedom used when computing the ``var``
and ``std`` properties. The divisor used in the calculation
is ``N - ddof``, where ``N`` is the number of unmasked pixels
within the aperture. The default is ``ddof=0``, which gives the
population variance and standard deviation. Use ``ddof=1`` to
obtain the sample (unbiased) variance and standard deviation.
This keyword affects only the ``var`` and ``std`` properties.
All other properties are unaffected, including the ``mean_err``
and ``median_err`` standard errors, which are always computed
using the sample standard deviation. Apertures with ``N <=
ddof`` unmasked pixels have an undefined ``var`` and ``std`` and
are set to NaN.
local_bkg : float, `~numpy.ndarray`, `~astropy.units.Quantity`, or `None`
The per-pixel local background values to subtract from the data
before performing measurements. If input as an array, the order
of ``local_bkg`` values corresponds to the order of the input
``aperture`` positions. ``local_bkg`` must have the same length
as the input ``aperture`` or must be a scalar value, which
will be broadcast to all apertures. If `None`, then no local
background subtraction is performed. If the input ``data`` has
units, then ``local_bkg`` must be a `~astropy.units.Quantity`
with the same units.
<segmentation_descriptions>
Notes
-----
``data`` should be background-subtracted for accurate source
properties. In addition to global background subtraction, local
background subtraction can be performed using the ``local_bkg``
keyword values.
`~regions.Region` objects are converted to `Aperture` objects using
the :func:`region_to_aperture` function.
The returned statistics are measured for the pixels within the
input aperture at its input position. This class does not change
the position of the input aperture. This class returns the centroid
value of the pixels within the input aperture, but the input
aperture is not recentered at the measured centroid position
when making the measurements. If desired, you can create a new
`Aperture` object using the measured centroid and then re-run
`~photutils.aperture.ApertureStats`.
All properties other than the sum-related ones described below are
calculated using the "center" aperture-mask method, which assigns
aperture weights of either 0 or 1, so the ``data`` pixel values
are used directly and without weighting. This choice reflects a
fundamental limitation because, unlike the mean or variance, order
statistics (``min``, ``max``, ``median``) and robust estimators
(``mad_std``, ``biweight_location``, ``biweight_midvariance``) have
no standard, unambiguous definition for pixels with fractional
(partial) aperture weights. Accordingly, these quantities cannot be
rigorously computed from a weighted aperture footprint.
The input ``sum_method`` and ``subpixels`` keywords are
used to determine the aperture-mask method only for the
sum-related properties: ``sum``, ``sum_err``, ``sum_aper_area``,
``data_sum_cutout``, and ``error_sum_cutout`` (also listed in the
``SUM_FOOTPRINT_PROPERTIES`` class attribute). All other properties,
including ``mean``, ``median``, ``std``, and the morphological
properties, always use the "center" aperture-mask method regardless
of ``sum_method``. The default is ``sum_method='exact'``, which
produces exact aperture-weighted photometry.
The sum-related properties have their own separate `sum_flags`
quality flags, distinct from the `flags` property used by all
other properties. To check for any quality issue across both
footprints, combine the two flag columns with a bitwise OR, e.g.,
``aperstats.flags | aperstats.sum_flags``.
The calculated statistics are always float64, regardless of the
input ``data`` dtype (`~astropy.units.Quantity` values with float64
dtype if the input ``data`` has units).
Examples
--------
>>> from photutils.datasets import make_4gaussians_image
>>> from photutils.aperture import CircularAperture, ApertureStats
>>> data = make_4gaussians_image()
>>> aper = CircularAperture((150, 25), 8)
>>> aperstats = ApertureStats(data, aper)
>>> print(aperstats.x_centroid)
149.99080259251238
>>> print(aperstats.y_centroid)
24.97484633000507
>>> print(aperstats.centroid)
[149.99080259 24.97484633]
>>> print(aperstats.mean, aperstats.median)
47.76300955780609 31.913789514433084
>>> print(aperstats.std)
39.193655383492974
>>> print(aperstats.sum)
9286.709206410273
>>> print(aperstats.sum_aper_area)
201.0619298297468 pix2
>>> # More than one aperture position
>>> aper2 = CircularAperture(((150, 25), (90, 60)), 10)
>>> aperstats2 = ApertureStats(data, aper2)
>>> print(aperstats2.x_centroid)
[149.98470724 89.97893946]
>>> print(aperstats2.sum)
[10177.62548482 36653.97704059]
"""
# Properties computed using the ``sum_method`` aperture-mask
# footprint (set by the ``sum_method``/``subpixels`` keywords).
# All other properties use the "center" aperture-mask method.
SUM_FOOTPRINT_PROPERTIES = ('sum', 'sum_err', 'sum_aper_area',
'data_sum_cutout', 'error_sum_cutout',
'sum_flags')
def __init__(self, data, aperture, *, error=None, mask=None, wcs=None,
sigma_clip=None, sum_method='exact', subpixels=5, ddof=0,
local_bkg=None, segmentation_image=None, labels=None,
mask_method='none'):
if isinstance(data, NDData):
data, error, mask, wcs = self._unpack_nddata(data, error, mask,
wcs)
inputs = (data, error, local_bkg)
names = ('data', 'error', 'local_bkg')
inputs, unit = process_quantities(inputs, names)
(data, error, local_bkg) = inputs
self._data = self._validate_array(data, 'data', shape=False)
self._data_unit = unit
self._input_aperture = self._validate_aperture(aperture)
aperture_meta = _aperture_metadata(aperture) # use input aperture
if isinstance(aperture, SkyAperture) and wcs is None:
msg = 'A wcs is required when using a SkyAperture'
raise ValueError(msg)
# Convert region to aperture if necessary
if not isinstance(aperture, Aperture):
aperture = region_to_aperture(aperture)
self.aperture = aperture
self._error = self._validate_array(error, 'error')
self._mask = self._validate_array(mask, 'mask')
self._wcs = wcs
if sigma_clip is not None and not isinstance(sigma_clip, SigmaClip):
msg = 'sigma_clip must be a SigmaClip instance'
raise TypeError(msg)
self.sigma_clip = sigma_clip
self.sum_method = sum_method
self.subpixels = subpixels
if (isinstance(ddof, bool)
or not isinstance(ddof, (int, np.integer)) or ddof < 0):
msg = 'ddof must be a non-negative integer'
raise ValueError(msg)
self.ddof = ddof
self._local_bkg = np.zeros(self.n_apertures) # no local bkg
if local_bkg is not None:
local_bkg = np.atleast_1d(local_bkg)
if local_bkg.ndim != 1:
msg = 'local_bkg must be a 1D array'
raise ValueError(msg)
n_local_bkg = len(local_bkg)
if n_local_bkg not in (1, self.n_apertures):
msg = ('local_bkg must be scalar or have the same length '
'as the input aperture')
raise ValueError(msg)
local_bkg = np.broadcast_to(local_bkg, self.n_apertures)
if np.any(~np.isfinite(local_bkg)):
msg = ('local_bkg must not contain any non-finite '
'(e.g., inf or NaN) values')
raise ValueError(msg)
self._local_bkg = local_bkg # always an iterable
self._ids = np.arange(self.n_apertures) + 1
self.default_columns = ['id', 'x_centroid', 'y_centroid',
'sky_centroid', 'sum', 'sum_err',
'sum_aper_area', 'sum_flags',
'center_aper_area', 'min', 'max', 'mean',
'median', 'mode', 'std', 'mad_std', 'var',
'biweight_location', 'biweight_midvariance',
'fwhm', 'semimajor_axis', 'semiminor_axis',
'orientation', 'eccentricity', 'flags']
self.meta = _get_meta()
self.meta.update(aperture_meta)
# Validate the segmentation-masking inputs and resolve the
# per-aperture source labels
self._mask_method = mask_method
seg_positions = np.atleast_2d(self._pixel_aperture.positions)
(self._segmentation,
self._seg_labels) = process_segmentation_inputs(
segmentation_image, labels, mask_method,
seg_positions, self._data.shape)
@staticmethod
def _unpack_nddata(data, error, mask, wcs):
nddata_attr = {'error': error, 'mask': mask, 'wcs': wcs}
for key, value in nddata_attr.items():
if value is not None:
msg = (f'The {key!r} keyword will be ignored. Its value '
'is obtained from the input NDData object.')
warnings.warn(msg, AstropyUserWarning)
mask = data.mask
wcs = data.wcs
if isinstance(data.uncertainty, StdDevUncertainty):
if data.uncertainty.unit is None:
error = data.uncertainty.array
else:
error = data.uncertainty.array * data.uncertainty.unit
if data.unit is not None:
data = u.Quantity(data.data, unit=data.unit)
else:
data = data.data
return data, error, mask, wcs
@staticmethod
def _validate_aperture(aperture):
try:
from regions import Region
aper_types = (Aperture, Region)
except ImportError:
aper_types = Aperture
if not isinstance(aperture, aper_types):
msg = 'aperture must be an Aperture or Region object'
raise TypeError(msg)
return aperture
def _validate_array(self, array, name, *, ndim=2, shape=True):
if name == 'mask' and array is np.ma.nomask:
array = None
if array is not None:
array = np.asanyarray(array)
if array.ndim != ndim:
msg = f'{name} must be a {ndim}D array'
raise ValueError(msg)
if shape and array.shape != self._data.shape:
msg = f'data and {name} must have the same shape'
raise ValueError(msg)
return array
@property
def _lazyproperties(self):
"""
A list of all class lazyproperties (even in superclasses).
The result is cached on the class to avoid repeated
introspection via `inspect.getmembers`.
"""
cls = self.__class__
attr = '_cached_lazyproperties'
# Subclasses get their own lazyproperty list
if attr not in cls.__dict__:
def islazyproperty(obj):
return isinstance(obj, lazyproperty)
setattr(cls, attr,
[i[0] for i in inspect.getmembers(
cls, predicate=islazyproperty)])
return getattr(cls, attr)
@property
def properties(self):
"""
A sorted list of the built-in source properties.
"""
lazyproperties = [name for name in self._lazyproperties if not
name.startswith('_')]
# isscalar and n_apertures are scalar values for the whole
# object, not per-source values, so they are not valid table
# columns
lazyproperties.remove('isscalar')
lazyproperties.remove('n_apertures')
lazyproperties.sort()
return lazyproperties
def __getitem__(self, index):
if self.isscalar:
msg = (f'A scalar {self.__class__.__name__!r} object cannot '
'be indexed')
raise TypeError(msg)
newcls = object.__new__(self.__class__)
# Attributes defined in __init__ that are copied directly to the
# new class
init_attr = ('_data', '_data_unit', '_error', '_mask', '_wcs',
'sigma_clip', 'sum_method', 'subpixels', 'ddof',
'default_columns', 'meta', '_segmentation',
'_mask_method')
for attr in init_attr:
setattr(newcls, attr, getattr(self, attr))
# aperture determines isscalar (needed below)
newcls.aperture = self.aperture[index]
# Keep _ids as a 1D array so a scalar instance's id is always
# backed by a length-1 iterable (see the id property).
newcls._ids = np.atleast_1d(self._ids[index])
# Slice the per-aperture segmentation labels
if self._seg_labels is None:
newcls._seg_labels = None
else:
newcls._seg_labels = np.atleast_1d(self._seg_labels[index])
# Slice evaluated lazyproperty objects
keys = set(self.__dict__.keys()) & set(self._lazyproperties)
keys.add('_local_bkg') # iterable defined in __init__
# The packed gather buffers and their reductions are not
# per-source sliceable; the sliced object recomputes them
# lazily from its sliced inputs.
keys.discard('_batch_inputs')
keys.discard('_fast_gather')
keys.discard('_fast_sum')
keys.discard('_sorted_values')
keys.discard('_order_stats')
keys.discard('_mean_var')
keys.discard('_mad')
keys.discard('_biweight')
keys.discard('_gini')
for key in keys:
value = self.__dict__[key]
# Do not insert attributes that are always scalar (e.g.,
# isscalar, n_apertures), i.e., not an array/list for each
# source
if np.isscalar(value):
continue
try:
# Keep most _<attrs> as length-1 iterables
if (newcls.isscalar and key.startswith('_')
and key != '_pixel_aperture'):
if isinstance(value, np.ndarray):
val = value[:, np.newaxis][index]
else:
val = [value[index]]
else:
val = value[index]
except TypeError:
# Apply fancy indices (e.g., array/list or bool mask) to
# lists.
# See https://numpy.org/doc/stable/release/1.20.0-notes.html
# #arraylike-objects-which-do-not-define-len-and-getitem
arr = np.empty(len(value), dtype=object)
arr[:] = list(value)
val = arr[index].tolist()
newcls.__dict__[key] = val
return newcls
def __str__(self):
cls_name = f'<{self.__class__.__module__}.{self.__class__.__name__}>'
with np.printoptions(threshold=25, edgeitems=5):
fmt = [f'Length: {self.n_apertures}']
return f'{cls_name}\n' + '\n'.join(fmt)
def __repr__(self):
return self.__str__()
def __len__(self):
if self.isscalar:
msg = f'Scalar {self.__class__.__name__!r} object has no len()'
raise TypeError(msg)
return self.n_apertures
def __iter__(self):
for item in range(len(self)):
yield self.__getitem__(item)
# Remove in 4.0
def __getattr__(self, name):
return deprecated_getattr(self, name, _DEPRECATED_ATTRIBUTES,
since='3.0', until='4.0')
def __getattribute__(self, name):
# Collapse the leading position axis of the public per-source
# output attributes to a scalar when a single scalar aperture
# position is input (e.g., ``CircularAperture((10, 20), r=5)``).
# The scalar conversion is applied centrally here instead of
# being repeated on each individual property (see ``_array`` for
# the array-preserving accessor used by the table-building and
# flag-decoding machinery).
value = super().__getattribute__(name)
if (not name.startswith('_')
and name not in _SCALAR_EXCLUDE
and isinstance(value, (np.ndarray, SkyCoord, list, tuple))
and self.isscalar):
try:
if len(value) == 1:
return value[0]
except TypeError: # value has no len
pass
return value
def _array(self, name):
"""
Return the full (never scalar-collapsed) value for a public
output attribute, regardless of whether the instance is scalar.
This bypasses the scalar conversion performed in
``__getattribute__`` so that the table-building and
flag-decoding machinery always operates on arrays.
"""
return object.__getattribute__(self, name)
@lazyproperty
def isscalar(self):
"""
Whether the instance is scalar (e.g., a single aperture
position).
"""
return self._pixel_aperture.isscalar
[docs]
def copy(self):
"""
Return a deep copy of this object.
Returns
-------
result : `ApertureStats`
A deep copy of this object.
"""
return deepcopy(self)
@lazyproperty
def _null_object(self):
"""
Return `None` values.
"""
return np.array([None] * self.n_apertures)
@lazyproperty
def _null_value(self):
"""
Return np.nan values.
"""
values = np.empty(self.n_apertures)
values.fill(np.nan)
return values
@property
def id(self):
"""
The aperture identification number(s).
"""
return self._ids
@property
@deprecated('3.1', alternative="the 'id' attribute", until='4.0')
def ids(self):
"""
The aperture identification number(s).
.. deprecated:: 3.1
Use the `id` attribute instead.
"""
return self.id
[docs]
def select_id(self, id_num):
"""
Return a new `ApertureStats` object for the input ID number
only.
Parameters
----------
id_num : int
The aperture ID number.
Returns
-------
result : `ApertureStats`
A new `ApertureStats` object containing only the source with
the input ID number.
"""
return self.select_ids(id_num)
[docs]
def select_ids(self, id_nums):
"""
Return a new `ApertureStats` object for the input ID numbers
only.
Parameters
----------
id_nums : list, tuple, or `~numpy.ndarray` of int
The aperture ID number(s).
Returns
-------
result : `ApertureStats`
A new `ApertureStats` object containing only the sources with
the input ID numbers.
"""
for id_num in np.atleast_1d(id_nums):
if id_num not in self._array('id'):
msg = f'{id_num} is not a valid source ID number'
raise ValueError(msg)
sorter = np.argsort(self.id)
indices = sorter[np.searchsorted(self.id, id_nums, sorter=sorter)]
return self[indices]
[docs]
@deprecated_positional_kwargs(since='3.0', until='4.0')
def to_table(self, *, columns=None):
"""
Create a `~astropy.table.QTable` of source properties.
Parameters
----------
columns : str, list of str, `None`, optional
Names of columns, in order, to include in the output
`~astropy.table.QTable`. The allowed column names are any of
the `ApertureStats` properties. If ``columns`` is `None`,
then a default list of scalar-valued properties (as defined
by the ``default_columns`` attribute) will be used.
Returns
-------
table : `~astropy.table.QTable`
A table of sources properties with one row per source.
Raises
------
ValueError
If any name in ``columns`` is not a valid (or deprecated)
column name.
"""
if columns is None:
table_columns = self.default_columns
else:
# id is not included in self.properties because it is not
# a lazyproperty
allowed_columns = set(self.properties) | set(self.default_columns)
# Remove 2D cutout images from the allowed columns
allowed_columns = {col for col in allowed_columns
if '_cutout' not in col}
deprecated_names = _DEPRECATED_ATTRIBUTES.copy()
# These are not valid column names, but are deprecated
# attributes that are still accessible in 3.x. They will be
# removed in 4.0.
invalid = ('data_sumcutout', 'error_sumcutout', 'get_id',
'get_ids')
for name in invalid:
deprecated_names.pop(name)
table_columns = validate_table_columns(
columns, allowed_columns, deprecated_names=deprecated_names)
# Replace with QTable in 4.0
tbl = create_empty_deprecated_qtable(
_DEPRECATED_ATTRIBUTES, since='3.0', until='4.0')
tbl.meta.update(self.meta) # keep tbl.meta type
# Evaluate the flag columns last (while preserving the
# requested column order in the output table). The ``flags``
# column reports the ``'singular_covariance'`` bit only if a
# covariance-derived shape property has already been computed
# (see `flags`). Evaluating the flag columns after all other
# columns ensures that any shape columns present in the same
# table are computed first, so the table's ``flags`` column
# consistently reflects the other columns it is shown alongside.
flag_columns = {'flags', 'sum_flags'}
eval_order = ([col for col in table_columns
if col not in flag_columns]
+ [col for col in table_columns if col in flag_columns])
values_map = {}
for column in eval_order:
values = getattr(self, column)
# Column assignment requires an object with a length
if self.isscalar:
values = (values,)
values_map[column] = values
for column in table_columns:
# Use the canonical (non-deprecated) column name so that
# assigning into ``tbl`` does not trigger a second,
# redundant deprecation warning (the deprecated attribute
# access above already warned once).
canonical_column = _DEPRECATED_ATTRIBUTES.get(column, column)
tbl[canonical_column] = values_map[column]
return tbl
@lazyproperty
def n_apertures(self):
"""
The number of positions for the input aperture.
"""
if self.isscalar:
return 1
return len(self._pixel_aperture)
@lazyproperty
def _pixel_aperture(self):
"""
The input aperture as a PixelAperture.
"""
if isinstance(self.aperture, SkyAperture):
return self.aperture.to_pixel(self._wcs)
return self.aperture
@lazyproperty
def _batch_inputs(self):
"""
The validated inputs for the fast Cython batch driver, or `None`.
Returns a tuple of the contiguous arrays and scalar parameters
shared by the center-value and ``sum_method`` gathers (see
`_fast_gather` and `_fast_sum`), together with the fast
sigma-clip specification (`None` when no clipping is performed).
`None` is returned (and the mask-based code path is used) when
the aperture or inputs are not supported by the batch driver,
for example when an unsupported ``sigma_clip`` is used, the
aperture does not opt in to the batch driver, or the data/mask
dtypes are not supported. The set of supported inputs matches
`~photutils.aperture.PixelAperture._batch_photometry`.
"""
# A non-None sigma_clip is only supported when it maps to the
# fast clipping kernel; otherwise fall back to the mask path.
clip_spec = self._fast_clip_spec()
if self.sigma_clip is not None and clip_spec is None:
return None
aper = self._pixel_aperture
# Use the batch driver only if the aperture's own class defines
# the _batch_shape_params hook (see _batch_photometry).
if '_batch_shape_params' not in type(aper).__dict__:
return None
spec = aper._batch_shape_params()
if spec is None:
return None
data = self._data
error = self._error
def _supported(arr):
return (type(arr) is np.ndarray and arr.dtype.kind in 'fiub'
and arr.dtype.itemsize <= 8)
if not _supported(data) or (error is not None
and not _supported(error)):
return None
mask = self._mask
if mask is not None and (not isinstance(mask, np.ndarray)
or mask.dtype != bool
or mask.shape != data.shape):
return None
# Fold non-finite ``data`` values into the mask plane so the
# batch kernels can skip the per-pixel finiteness test (and
# defer the pixel-value load to contributing pixels only, like
# the ``photometry`` method). This matches the mask-based
# path, which masks non-finite data before any segmentation
# correction. The plane distinguishes the two causes for the
# flag counts: bit 1 (value 1) marks input-masked pixels and bit
# 2 (value 2) marks non-finite data pixels. Any nonzero value
# excludes the pixel.
plane = None
if mask is not None:
plane = mask.astype(np.uint8)
if data.dtype.kind == 'f':
nonfinite = ~np.isfinite(data)
if nonfinite.any():
if plane is None:
plane = np.zeros(data.shape, dtype=np.uint8)
plane[nonfinite] = 2
else:
plane[nonfinite & (plane == 0)] = 2
mask = plane
if mask is not None:
mask = np.ascontiguousarray(mask)
seg_arr = None
labels_arr = None
seg_code = 0
if self._segmentation is not None and self._mask_method != 'none':
seg_arr = np.ascontiguousarray(self._segmentation, dtype=np.intp)
labels_arr = np.ascontiguousarray(self._seg_labels, dtype=np.intp)
seg_code = SEG_METHOD_CODES[self._mask_method]
shape_code, params = spec
sum_use_exact, sum_subpixels = aper._translate_mask_method(
self.sum_method, self.subpixels)
ext_x, ext_y = aper._xy_extents
if error is not None:
error = np.ascontiguousarray(error, dtype=np.float64)
return (np.ascontiguousarray(data, dtype=np.float64), error, mask,
np.ascontiguousarray(aper._positions, dtype=np.float64),
shape_code, np.array(params, dtype=np.float64),
float(ext_x), float(ext_y), sum_use_exact, sum_subpixels,
np.ascontiguousarray(self._local_bkg, dtype=np.float64),
seg_arr, labels_arr, seg_code, clip_spec)
@lazyproperty
def _fast_gather(self):
"""
The fast Cython "center"-method value gather, or `None`.
When supported, this walks each aperture bounding box once and
returns a packed buffer of the unmasked "center"-method pixel
values (and their cutout coordinates). The higher-tier
statistics (order statistics and image moments) are reduced
lazily from this buffer. The ``sum_method`` aperture sum, error
variance, and area are computed separately and lazily by
`_fast_sum`, so they are not paid for unless requested.
When a supported ``sigma_clip`` is set, the packed center buffer
is sigma-clipped per source before being returned, so the
downstream reductions operate on the clipped data transparently.
`None` is returned (and the mask-based code path is used) when
the fast batch driver is unavailable (see `_batch_inputs`).
The result is a tuple ``(values, local_x, local_y, starts,
counts, sum_aper, var_aper, sum_area, overlap, ...)`` (a
14-tuple shared with `_fast_sum`, whose last entry is the
per-source flag-count array). Only the center-value entries
(indices 0-4), ``overlap`` (index 8), and the flag counts (index
13) are populated here. The ``sum_method`` entries are `None`.
"""
inputs = self._batch_inputs
if inputs is None:
return None
(data, _error, mask, positions, shape_code, params, ext_x, ext_y,
_sum_use_exact, _sum_subpixels, local_bkg, seg_arr, labels_arr,
seg_code, clip_spec) = inputs
(values, lx, ly, starts, counts, overlap,
flag_counts) = batch_aperture_gather(
data, mask, positions, shape_code, params, ext_x, ext_y,
local_bkg, seg_arr, labels_arr, seg_code)
gather = (values, lx, ly, starts, counts, None, None, None, overlap,
None, None, None, None, flag_counts)
if clip_spec is not None:
gather = self._apply_center_clip(gather, clip_spec)
return gather
@lazyproperty
def _fast_sum(self):
"""
The fast Cython ``sum_method`` aperture gather, or `None`.
When supported, this walks each aperture bounding box once and
returns the ``sum_method`` aperture sum, error variance, and
area. This is computed only when the ``sum``, ``sum_err``, or
``sum_aper_area`` properties are requested, so the more
expensive exact/subpixel overlap fractions are not evaluated for
the value-statistics-only use cases.
When a supported ``sigma_clip`` is set, the packed ``sum_method``
members are sigma-clipped per source and the aperture sum,
variance, and area are recomputed over the survivors.
`None` is returned (and the mask-based code path is used) when
the fast batch driver is unavailable (see `_batch_inputs`).
The result is a tuple with the same layout as `_fast_gather`;
only the ``sum_aper`` (index 5), ``var_aper`` (index 6),
``sum_area`` (index 7), and ``overlap`` (index 8) entries are
populated.
"""
inputs = self._batch_inputs
if inputs is None:
return None
(data, error, mask, positions, shape_code, params, ext_x, ext_y,
sum_use_exact, sum_subpixels, local_bkg, seg_arr, labels_arr,
seg_code, clip_spec) = inputs
emit_sum = 1 if clip_spec is not None else 0
(sums, sum_var, area, overlap, starts, sum_values, sum_fracs,
sum_errsq, scounts, flag_counts) = batch_aperture_sums(
data, error, mask, positions, shape_code, params, ext_x, ext_y,
sum_use_exact, sum_subpixels, seg_arr, labels_arr, seg_code,
local_bkg, emit_sum)
gather = (None, None, None, starts, None, sums, sum_var, area,
overlap, sum_values, sum_fracs, sum_errsq, scounts,
flag_counts)
if clip_spec is not None:
gather = self._apply_sum_clip(gather, clip_spec)
return gather
def _fast_clip_spec(self):
"""
The fast sigma-clip parameters, or `None`.
Returns ``(sigma_lower, sigma_upper, maxiters, cenfunc_code,
stdfunc_code)`` when ``sigma_clip`` is a `SigmaClip` instance
supported by the fast clipping kernel, and `None` otherwise (in
which case the mask-based path is used). The fast path supports
the string ``cenfunc`` values 'median'/'mean', the string
``stdfunc`` values 'std'/'mad_std', and no spatial growing.
"""
sc = self.sigma_clip
if sc is None or not isinstance(sc, SigmaClip):
return None
if not (isinstance(sc.cenfunc, str) and sc.cenfunc in ('median',
'mean')):
return None
if not (isinstance(sc.stdfunc, str) and sc.stdfunc in ('std',
'mad_std')):
return None
if sc.grow: # False or 0 is supported; spatial growing is not
return None
maxiters = sc.maxiters
if maxiters is None or maxiters == np.inf:
maxiters = -1
else:
maxiters = int(maxiters)
cenfunc_code = 0 if sc.cenfunc == 'median' else 1
stdfunc_code = 0 if sc.stdfunc == 'std' else 1
return (float(sc.sigma_lower), float(sc.sigma_upper), maxiters,
cenfunc_code, stdfunc_code)
def _apply_center_clip(self, gather, clip_spec):
"""
Sigma-clip the packed center buffer and return a gather tuple
with the same layout whose center buffer reflects the per-source
clipped data.
"""
sigma_lower, sigma_upper, maxiters, cenfunc, stdfunc = clip_spec
(values, local_x, local_y, starts, counts, _sum, _var, _area,
overlap, _sv, _sf, _se, _sc, flag_counts) = gather
(cvalues, clx, cly, cstarts, ccounts) = batch_sigma_clip_center(
values, local_x, local_y, starts, counts, sigma_lower,
sigma_upper, maxiters, cenfunc, stdfunc)
return (cvalues, clx, cly, cstarts, ccounts, _sum, _var, _area,
overlap, None, None, None, None, flag_counts)
def _apply_sum_clip(self, gather, clip_spec):
"""
Sigma-clip the packed ``sum_method`` member buffers and return a
gather tuple with the same layout whose aperture sum, variance,
and area reflect the per-source clipped data.
"""
sigma_lower, sigma_upper, maxiters, cenfunc, stdfunc = clip_spec
(values, local_x, local_y, starts, counts, _sum, _var, _area,
overlap, sum_values, sum_fracs, sum_errsq, sum_counts,
flag_counts) = gather
has_error = 1 if self._error is not None else 0
(csum, cvar, carea) = batch_sigma_clip_sum(
sum_values, sum_fracs, sum_errsq, starts, sum_counts,
sigma_lower, sigma_upper, maxiters, cenfunc, stdfunc, has_error)
return (values, local_x, local_y, starts, counts, csum, cvar, carea,
overlap, None, None, None, None, flag_counts)
@lazyproperty
def _sorted_values(self):
"""
The packed per-source ascending-sorted center pixel values, or
`None`.
When the fast gather path is active (see `_fast_gather`), each
source's packed pixel values are sorted once. The sorted buffer
is cached and shared by the order statistics (``min``, ``max``,
``median``), ``mad_std``, and the biweight estimators, so the
per-source sort is performed only once. `None` is returned when
the fast path is unavailable.
"""
gather = self._fast_gather
if gather is None:
return None
values, starts, counts = gather[0], gather[3], gather[4]
return (batch_sort_values(values, starts, counts), starts, counts)
@lazyproperty
def _order_stats(self):
"""
The per-source ``(min, max, median)`` arrays, or `None`.
Reduced from the cached sorted buffer (see `_sorted_values`).
`None` when the fast path is unavailable.
"""
sorted_values = self._sorted_values
if sorted_values is None:
return None
return batch_order_stats(*sorted_values)
@lazyproperty
def _mean_var(self):
"""
The per-source ``(mean, var)`` arrays, or `None`.
The population (``ddof=0``) variance is returned. Computed
directly from the packed center buffer (no sort required).
`None` when the fast path is unavailable.
"""
gather = self._fast_gather
if gather is None:
return None
values, starts, counts = gather[0], gather[3], gather[4]
return batch_mean_var(values, starts, counts)
@lazyproperty
def _mad(self):
"""
The per-source unscaled median absolute deviation, or `None`.
Reduced from the cached sorted buffer (see `_sorted_values`).
`None` when the fast path is unavailable.
"""
sorted_values = self._sorted_values
if sorted_values is None:
return None
return batch_mad(*sorted_values)
@lazyproperty
def _biweight(self):
"""
The per-source ``(biweight_location, biweight_midvariance)``
arrays, or `None`.
Reused from the cached sorted buffer, median (`_order_stats`),
and unscaled MAD (`_mad`). `None` when the fast path is
unavailable.
"""
sorted_values = self._sorted_values
if sorted_values is None:
return None
_, _, median = self._order_stats
return batch_biweight(*sorted_values, median, self._mad)
@lazyproperty
def _gini(self):
"""
The per-source Gini coefficient, or `None`.
Computed from the packed center buffer (the absolute values are
sorted internally). `None` when the fast path is unavailable.
"""
gather = self._fast_gather
if gather is None:
return None
values, starts, counts = gather[0], gather[3], gather[4]
return batch_gini(values, starts, counts)
def _finalize_value_stat(self, fast_result, stat_func, *,
square_unit=False, apply_unit=True):
"""
Return a per-source value statistic, using the fast Cython
reduction when available and otherwise the mask-based per-source
code path.
Parameters
----------
fast_result : `~numpy.ndarray` or `None`
The per-source statistic from the fast reduction, or `None`
to use the mask-based fallback.
stat_func : callable
The fallback callable applied to each source's 1D array of
unmasked pixel values when the fast path is unavailable.
square_unit : bool, optional
Whether the statistic has squared data units (e.g. variance).
apply_unit : bool, optional
Whether to apply the data unit at all (`False` for
dimensionless statistics such as the Gini coefficient).
"""
if fast_result is not None:
result = fast_result.copy()
else:
result = np.array([stat_func(arr)
for arr in self._data_values_center])
if apply_unit:
unit = self._data_unit
if unit is not None:
if square_unit:
unit = unit**2
result <<= unit
return result
@lazyproperty
def _aperture_masks_center(self):
"""
The aperture masks (`ApertureMask`) generated with the 'center'
method, always as an iterable.
"""
aperture_masks = self._pixel_aperture.to_mask(method='center')
if self.isscalar:
aperture_masks = (aperture_masks,)
return aperture_masks
@lazyproperty
def _aperture_masks(self):
"""
The aperture masks (`ApertureMask`) generated with the
``sum_method`` method, always as an iterable.
"""
aperture_masks = self._pixel_aperture.to_mask(method=self.sum_method,
subpixels=self.subpixels)
if self.isscalar:
aperture_masks = (aperture_masks,)
return aperture_masks
@lazyproperty
def _overlap_slices(self):
"""
The aperture mask overlap slices with the data, always as an
iterable.
The overlap slices are the same for all aperture mask methods.
"""
overlap_slices = []
for apermask in self._aperture_masks_center:
(slc_large, slc_small) = apermask.get_overlap_slices(
self._data.shape)
overlap_slices.append((slc_large, slc_small))
return overlap_slices
@lazyproperty
def _data_cutouts(self):
"""
The local-background-subtracted unmasked data cutouts using the
aperture bounding box, always as an iterable.
"""
cutouts = []
for (slices, local_bkg) in zip(self._overlap_slices,
self._local_bkg, strict=True):
if slices[0] is None:
cutout = None # no aperture overlap with the data
else:
# Copy is needed to preserve input data because masks are
# applied to these cutouts later
cutout = (self._data[slices[0]].astype(float, copy=True)
- local_bkg)
cutouts.append(cutout)
return cutouts
def _make_aperture_cutouts(self, aperture_masks, *, count_clipped=True):
"""
Make aperture-weighted cutouts for the data and variance, and
cutouts for the total mask and aperture mask weights.
Parameters
----------
aperture_masks : list of `ApertureMask`
A list of `ApertureMask` objects.
count_clipped : bool, optional
Whether to count the number of sigma-clipped pixels per
source. Only `_footprint_flag_inputs` for the "center"
footprint uses this count, so it is skipped (left at 0)
for the ``sum_method`` footprint to avoid the wasted
computation.
Returns
-------
result : list of `~numpy.ndarray`
A list of cutout arrays for the data, variance,
mask and weight arrays for each source (aperture
position), followed by the overlap indicator, per-source
flag counts (see the ``FLAG_COL_*`` constants in
`photutils.aperture._batch_photometry`, with semantics
identical to the batch drivers), and the number of
sigma-clipped pixels (always 0 when ``count_clipped`` is
`False`).
"""
data_cutouts = []
variance_cutouts = []
mask_cutouts = []
weight_cutouts = []
overlaps = []
flag_counts = []
n_clipped = []
positions = np.atleast_2d(self._pixel_aperture.positions)
for idx, (data_cutout, apermask, slices) in enumerate(
zip(self._data_cutouts, aperture_masks,
self._overlap_slices, strict=True)):
fc_row = np.zeros(8, dtype=np.intp)
n_clipped_ = 0
slc_large, slc_small = slices
if slc_large is None: # aperture does not overlap the data
overlap = False
data_cutout = np.array([np.nan])
variance_cutout = np.array([np.nan])
mask_cutout = np.array([False])
weight_cutout = np.array([np.nan])
else:
# Create a mask of non-finite ``data`` values combined
# with the input ``mask`` array
nonfinite_mask = ~np.isfinite(data_cutout)
if self._mask is not None:
user_mask = self._mask[slc_large]
data_mask = nonfinite_mask | user_mask
else:
user_mask = None
data_mask = nonfinite_mask
error_cutout = (None if self._error is None
else self._error[slc_large])
# Apply segmentation-based masking and/or symmetric
# neighbor correction
exclude = None
affected = None
if (self._segmentation is not None
and self._mask_method != 'none'):
segm_cutout = self._segmentation[slc_large]
cutout_xycen = (positions[idx, 0] - slc_large[1].start,
positions[idx, 1] - slc_large[0].start)
(data_cutout, error_cutout, exclude,
affected) = make_segmentation_exclusion(
self._mask_method, segm_cutout,
self._seg_labels[idx], data=data_cutout,
error=error_cutout, base_mask=data_mask,
cutout_xycen=cutout_xycen)
pre_seg_mask = data_mask
data_mask = data_mask | exclude
overlap = True
aperweight_cutout = apermask.data[slc_small]
weight_cutout = aperweight_cutout * ~data_mask
# Per-source pixel counts for the quality flags,
# matching the batch-driver semantics (computed before
# any sigma clipping)
weighted = aperweight_cutout > 0
fc_row[FLAG_COL_NPIX] = np.count_nonzero(weighted)
fc_row[FLAG_COL_BBOX_CLIPPED] = (
aperweight_cutout.shape != apermask.data.shape)
if user_mask is not None:
fc_row[FLAG_COL_MASKED] = np.count_nonzero(
weighted & user_mask)
fc_row[FLAG_COL_NONFINITE_DATA] = np.count_nonzero(
weighted & nonfinite_mask & ~user_mask)
else:
fc_row[FLAG_COL_NONFINITE_DATA] = np.count_nonzero(
weighted & nonfinite_mask)
if affected is not None:
pre_valid = weighted & ~pre_seg_mask
fc_row[FLAG_COL_SEG] = np.count_nonzero(
affected & pre_valid)
if self._mask_method == 'correct':
# In 'correct' mode, the excluded pixels are
# exactly the uncorrectable neighbor pixels
fc_row[FLAG_COL_UNCORRECTED] = np.count_nonzero(
exclude & pre_valid)
fc_row[FLAG_COL_VALID] = np.count_nonzero(
weighted & ~data_mask)
if error_cutout is not None:
fc_row[FLAG_COL_NONFINITE_ERROR] = np.any(
~np.isfinite(error_cutout) & weighted & ~data_mask)
# Apply the aperture mask; for "exact" and "subpixel"
# this is an expanded boolean mask using the aperture
# mask zero values
mask_cutout = (aperweight_cutout == 0) | data_mask
data_cutout = data_cutout.copy()
if self.sigma_clip is None:
# data_cutout will have zeros where mask_cutout is True
data_cutout *= ~mask_cutout
else:
# To input a mask, SigmaClip needs a MaskedArray
data_cutout_ma = np.ma.masked_array(data_cutout,
mask=mask_cutout)
data_sigclip = self.sigma_clip(data_cutout_ma)
# Define a mask of only the sigma-clipped pixels
sigclip_mask = data_sigclip.mask & ~mask_cutout
if count_clipped:
n_clipped_ = np.count_nonzero(sigclip_mask)
weight_cutout *= ~sigclip_mask
mask_cutout = data_sigclip.mask
data_cutout = data_sigclip.filled(0.0)
# Need to apply the aperture weights
data_cutout *= aperweight_cutout
if self._error is None:
variance_cutout = None
else:
# Apply the exact weights and total mask;
# error_cutout will have zeros where mask_cutout is True
variance = error_cutout**2
variance_cutout = (variance * aperweight_cutout
* ~mask_cutout)
data_cutouts.append(data_cutout)
variance_cutouts.append(variance_cutout)
mask_cutouts.append(mask_cutout)
weight_cutouts.append(weight_cutout)
overlaps.append(overlap)
flag_counts.append(fc_row)
n_clipped.append(n_clipped_)
# Use zip (instead of np.transpose) because these may contain
# arrays that have different shapes
return list(zip(data_cutouts, variance_cutouts, mask_cutouts,
weight_cutouts, overlaps, flag_counts, n_clipped,
strict=True))
@lazyproperty
def _aperture_cutouts_center(self):
"""
Aperture-weighted cutouts for the data, variance, total mask,
and aperture weights using the "center" aperture mask method.
"""
return self._make_aperture_cutouts(self._aperture_masks_center)
@lazyproperty
def _aperture_cutouts(self):
"""
Aperture-weighted cutouts for the data, variance, total mask,
and aperture weights using the input ``sum_method`` aperture
mask method.
"""
# The sigma-clipped pixel count is not tracked for this
# footprint: only the "center" footprint's flags use it (see
# `_footprint_flag_inputs`), so counting it here would be wasted
# computation.
return self._make_aperture_cutouts(self._aperture_masks,
count_clipped=False)
@lazyproperty
def _mask_cutout_center(self):
"""
Boolean mask cutouts representing the total mask.
The total mask is combination of the input ``mask``, non-finite
``data`` values, the cutout aperture mask using the "center"
method, and the sigma-clip mask.
"""
return list(zip(*self._aperture_cutouts_center, strict=True))[2]
@lazyproperty
def _mask_cutout(self):
"""
Boolean mask cutouts representing the total mask.
The total mask is combination of the input ``mask``,
non-finite ``data`` values, the cutout aperture mask using the
``sum_method`` method, and the sigma-clip mask.
"""
return list(zip(*self._aperture_cutouts, strict=True))[2]
def _make_masked_array_center(self, array):
"""
Return a list of cutout masked arrays using the ``_mask_cutout``
mask.
Units are not applied.
"""
return [np.ma.masked_array(arr, mask=mask)
for arr, mask in zip(array, self._mask_cutout_center,
strict=True)]
def _make_masked_array(self, array):
"""
Return a list of cutout masked arrays using the
``_mask_sumcutout`` mask.
Units are not applied.
"""
return [np.ma.masked_array(arr, mask=mask)
for arr, mask in zip(array, self._mask_cutout, strict=True)]
@lazyproperty
def data_cutout(self):
"""
A 2D aperture-weighted cutout from the data using the aperture
mask with the "center" method as a `~numpy.ma.MaskedArray`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
return self._make_masked_array_center(
list(zip(*self._aperture_cutouts_center, strict=True))[0])
@lazyproperty
def data_sum_cutout(self):
"""
A 2D aperture-weighted cutout from the data using the aperture
mask with the input ``sum_method`` method as a
`~numpy.ma.MaskedArray`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
return self._make_masked_array(list(zip(*self._aperture_cutouts,
strict=True))[0])
@lazyproperty
def _variance_cutout_center(self):
"""
A 2D aperture-weighted variance cutout using the aperture mask
with the input "center" method as a `~numpy.ma.MaskedArray`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
if self._error is None:
return self._null_object
return self._make_masked_array_center(
list(zip(*self._aperture_cutouts_center, strict=True))[1])
@lazyproperty
def _variance_cutout(self):
"""
A 2D aperture-weighted variance cutout using the aperture mask
with the input ``sum_method`` method as a
`~numpy.ma.MaskedArray`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
if self._error is None:
return self._null_object
return self._make_masked_array(list(zip(*self._aperture_cutouts,
strict=True))[1])
@lazyproperty
def error_sum_cutout(self):
"""
A 2D aperture-weighted error cutout using the aperture mask with
the input ``sum_method`` method as a `~numpy.ma.MaskedArray`.
The cutout does not have units due to current limitations of
masked quantity arrays.
The mask is `True` for pixels from the input ``mask``,
non-finite ``data`` values (NaN and inf), sigma-clipped pixels
within the aperture, and pixels where the aperture mask has zero
weight.
"""
if self._error is None:
return self._null_object
return [np.sqrt(var) for var in self._variance_cutout]
@lazyproperty
def _weight_cutout_center(self):
"""
A 2D `~numpy.ma.MaskedArray` cutout from the aperture mask
weights array using the aperture bounding box.
The aperture mask weights are for the "center" method.
The mask is `True` for pixels outside the aperture mask, pixels
from the input ``mask``, non-finite ``data`` values (NaN and
inf), and sigma-clipped pixels.
"""
return self._make_masked_array_center(
list(zip(*self._aperture_cutouts_center, strict=True))[3])
@lazyproperty
def _weight_cutout(self):
"""
A 2D `~numpy.ma.MaskedArray` cutout from the aperture mask
weights array using the aperture bounding box.
The aperture mask weights are for the ``sum_method`` method.
The mask is `True` for pixels outside the aperture mask, pixels
from the input ``mask``, non-finite ``data`` values (NaN and
inf), and sigma-clipped pixels.
"""
return self._make_masked_array(list(zip(*self._aperture_cutouts,
strict=True))[3])
@lazyproperty
def _moment_data_cutout(self):
"""
A list of 2D `~numpy.ndarray` cutouts from the data.
Masked pixels are set to zero in these arrays (zeros do not
contribute to the image moments). The aperture mask weights are
for the "center" method.
These arrays are used to derive moment-based properties.
"""
data = deepcopy(self._array('data_cutout'))
cutouts = []
for arr in data:
if arr.size == 1 and np.isnan(arr[0]): # no aperture overlap
arr_ = np.empty((2, 2))
arr_.fill(np.nan)
else:
arr_ = arr.data
arr_[arr.mask] = 0.0
cutouts.append(arr_)
return cutouts
@lazyproperty
def _all_masked(self):
"""
True if all pixels within the aperture are masked.
"""
return np.array([np.all(mask) for mask in self._mask_cutout_center])
@lazyproperty
def _overlap(self):
"""
True if there is no overlap of the aperture with the data.
"""
return list(zip(*self._aperture_cutouts_center, strict=True))[4]
def _footprint_flag_inputs(self, footprint):
"""
The per-source flag inputs for the given aperture footprint.
Parameters
----------
footprint : {'center', 'sum'}
The aperture footprint: the "center" mask method used by the
value statistics, or the ``sum_method`` mask method used by
the sum properties.
Returns
-------
flag_counts : `~numpy.ndarray`
The per-source flag counts (see the ``FLAG_COL_*``
constants in `photutils.aperture._batch_photometry`).
overlap : `~numpy.ndarray` (bool)
Whether the aperture bounding box overlaps the data.
n_kept : `~numpy.ndarray` or `None`
For the center footprint, the per-source number of valid
pixels remaining after sigma clipping; `None` for the sum
footprint.
"""
if footprint == 'center':
gather = self._fast_gather
if gather is not None:
return (gather[13], np.asarray(gather[8]),
np.asarray(gather[4]))
cutouts = self._aperture_cutouts_center
else:
gather = self._fast_sum
if gather is not None:
return gather[13], np.asarray(gather[8]), None
cutouts = self._aperture_cutouts
flag_counts = np.array([cut[5] for cut in cutouts])
overlap = np.array([cut[4] for cut in cutouts])
n_kept = None
if footprint == 'center':
n_clipped = np.array([cut[6] for cut in cutouts])
n_kept = flag_counts[:, FLAG_COL_VALID] - n_clipped
return flag_counts, overlap, n_kept
def _footprint_flags(self, footprint):
"""
The per-source flag bits for the given aperture footprint.
The bounding-box-clipped candidate sources are resolved to the
precise outside-weight test using per-source aperture masks.
"""
flag_counts, overlap, _ = self._footprint_flag_inputs(footprint)
if footprint == 'center':
method, subpixels = 'center', 1
else:
method, subpixels = self.sum_method, self.subpixels
candidates = flag_counts[:, FLAG_COL_BBOX_CLIPPED].astype(bool)
w_out = self._pixel_aperture._resolve_outside_weights(
self._data.shape, method=method, subpixels=subpixels,
candidates=candidates)
return _counts_to_flag_bits(flag_counts, overlap, w_out)
@lazyproperty
def _base_flags(self):
"""
The count-based value-statistics flag bits (1D int array).
These are the flag bits that do not depend on any lazily
computed covariance property: the "center"-method
footprint bits plus the sigma-clip and ``ddof`` bits. The
``singular_covariance`` bit is folded in by the `flags` property
only if a covariance-derived property has been computed.
"""
# The gather kernel and the center-method cutouts do not
# evaluate error values, so the non-finite-error bit is defined
# by the sum footprint only (see `sum_flags`).
flags = (self._footprint_flags('center')
& ~APERTURE_FLAGS.NON_FINITE_ERROR)
flag_counts, _, n_kept = self._footprint_flag_inputs('center')
n_valid = flag_counts[:, FLAG_COL_VALID]
if self.sigma_clip is not None:
n_clipped = n_valid - n_kept
flags[n_clipped > 0] |= APERTURE_FLAGS.SIGMA_CLIPPED
flags[(n_valid > 0)
& (n_kept == 0)] |= APERTURE_FLAGS.ALL_CLIPPED
if self.ddof > 0:
w_in = flag_counts[:, FLAG_COL_NPIX]
flags[(w_in > 0)
& (n_kept <= self.ddof)] |= APERTURE_FLAGS.TOO_FEW_PIXELS
return flags
@lazyproperty
def _singular_covariance_mask(self):
"""
Boolean mask (1D) marking sources with a singular or nearly
singular covariance matrix.
A source is flagged when the minor-axis variance (the smaller
eigenvalue of its normalized second-moment covariance matrix)
falls below ``1/12``. This flags both unresolved, nearly
point-like sources, where the covariance determinant drops below
``(1/12)**2``, and rank-1 degenerate sources, where one axis is
unresolved while the other is extended (which the determinant
alone would miss). Sources with undefined moments (no overlap or
fully masked) have a non-finite determinant and are not flagged
here. They are already reported by the overlap and masking bits.
"""
covar = self._raw_covariance
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
covar_det = np.linalg.det(covar)
# Smaller eigenvalue (the minor-axis variance) of each 2x2
# symmetric covariance matrix, via the closed form
# lambda = tr/2 -/+ sqrt((tr/2)**2 - det). The discriminant
# ((lambda1 - lambda2)/2)**2 is non-negative for a real
# symmetric matrix; clip tiny negative rounding to zero.
half_trace = 0.5 * (covar[:, 0, 0] + covar[:, 1, 1])
disc = np.maximum(half_trace**2 - covar_det, 0.0)
min_eigval = half_trace - np.sqrt(disc)
# ``1/12`` is the variance of a uniform distribution across a
# single pixel, and hence the smallest second moment a resolved
# source can have given finite pixel size. The determinant test
# (det < (1/12)**2) flags the isotropic case where both axes are
# unresolved; the eigenvalue test (min_eigval < 1/12)
# additionally flags rank-1 degeneracy (one unresolved axis) and
# covariance matrices that are not positive semidefinite
# (det < 0).
delta = 1.0 / 12
finite = np.isfinite(covar_det) & np.isfinite(min_eigval)
return finite & ((covar_det < delta**2) | (min_eigval < delta))
[docs]
@_UncachedLazyProperty
@_update_method_subpixels_docstring
def flags(self):
# numpydoc ignore: RT01
"""
The bitwise quality flags for the value statistics.
The flags are evaluated on the "center"-method footprint used
by the value statistics (e.g., ``mean``, ``median``, ``std``).
The sum properties (``sum``, ``sum_err``, and ``sum_aper_area``)
have their own separate `sum_flags`. The ``'sigma_clipped'``,
``'all_clipped'``, and ``'too_few_pixels'`` flags are evaluated
on this footprint. To check for any quality issue across both
footprints, combine the two flag columns with a bitwise OR
(e.g., ``flags | sum_flags``).
The ``'singular_covariance'`` bit is special. It reports whether
a source's covariance matrix is singular or nearly singular,
a condition that is only knowable once a covariance-derived
shape property (e.g., ``semimajor_axis``, ``orientation``,
``eccentricity``) has been computed. To avoid forcing that
computation, the bit is included only if such a property has
already been evaluated on this object; otherwise it is omitted.
This means the value of ``flags`` reflects the measurements
requested so far, so accessing a shape property and then
re-reading ``flags`` may set additional bits. The default
`to_table` always evaluates the shape properties, so its
``flags`` column always reflects the ``'singular_covariance'``
bit.
See `~photutils.aperture.decode_aperture_flags` for decoding
flag values. The flags are:
<flag_descriptions>
"""
# A fresh array is returned on every access (never mutating the
# cached `_base_flags`), so concurrent readers and previously
# returned arrays are never modified in place.
flags = self._base_flags.copy()
if '_covariance' in self.__dict__:
flags[self._singular_covariance_mask] |= (
APERTURE_FLAGS.SINGULAR_COVARIANCE)
return flags
@lazyproperty
@_update_method_subpixels_docstring
def sum_flags(self):
# numpydoc ignore: RT01
"""
The bitwise quality flags for the sum properties.
The flags are evaluated on the ``sum_method`` footprint
used by the sum properties (``sum``, ``sum_err``, and
``sum_aper_area``). The value statistics have their own separate
`flags`. The ``'non_finite_error'`` flag is evaluated on this
footprint. The ``'sigma_clipped'``, ``'all_clipped'``, and
``'too_few_pixels'`` flags apply only to the value statistics
and are never set here. To check for any quality issue across
both footprints, combine the two flag columns with a bitwise OR
(e.g., ``flags | sum_flags``).
See `~photutils.aperture.decode_aperture_flags` for decoding
flag values. The flags are:
<flag_descriptions>
"""
return self._footprint_flags('sum')
[docs]
def decode_flags(self, *, column='flags', return_bit_values=False):
"""
Decode the source quality flags into individual components.
This is a convenience method that calls
`~photutils.aperture.decode_aperture_flags` with the `flags` or
`sum_flags` property.
Parameters
----------
column : {'flags', 'sum_flags'}, optional
Which quality flags to decode: ``'flags'`` for the value
statistics (default) or ``'sum_flags'`` for the sum
properties.
return_bit_values : bool, optional
If `True`, return the decoded bit flags (integers) instead
of the flag names (strings).
Returns
-------
decoded : list of list of str or list of list of int
A list of the active flag names (or bit values) for each
source.
See Also
--------
photutils.aperture.decode_aperture_flags
Examples
--------
>>> import numpy as np
>>> from photutils.aperture import ApertureStats, CircularAperture
>>> data = np.ones((25, 25))
>>> mask = np.zeros(data.shape, dtype=bool)
>>> mask[12, 12] = True
>>> aper = CircularAperture([(12.0, 12.0), (0.0, 12.0)], r=3.0)
>>> aperstats = ApertureStats(data, aper, mask=mask)
>>> for names in aperstats.decode_flags():
... print(names)
['masked_pixels']
['partial_overlap']
"""
if column not in ('flags', 'sum_flags'):
msg = "column must be 'flags' or 'sum_flags'"
raise ValueError(msg)
return decode_aperture_flags(self._array(column),
return_bit_values=return_bit_values)
def _get_values(self, array):
"""
Get a 1D array of unmasked aperture-weighted values from the
input array.
An array with a single NaN is returned for completely-masked
sources.
"""
if self.isscalar:
array = (array,)
return [arr.compressed() if len(arr.compressed()) > 0
else np.array([np.nan]) for arr in array]
@lazyproperty
def _data_values_center(self):
"""
A 1D array of unmasked aperture-weighted data values using the
"center" method.
An array with a single NaN is returned for completely-masked
sources.
"""
return self._get_values(self.data_cutout)
@lazyproperty
def moments(self):
"""
Spatial moments up to 3rd order of the source.
"""
gather = self._fast_gather
if gather is not None:
values, lx, ly, starts, counts = gather[:5]
overlap = gather[8]
zeros = np.zeros(self.n_apertures)
mom = batch_moments(values, lx, ly, starts, counts, zeros, zeros)
# No-overlap sources have NaN moments (the mask-based path
# uses an all-NaN cutout); all-masked overlapping sources
# have zero moments (an all-zero cutout).
mom[~overlap] = np.nan
return mom
return np.array([_image_moments(arr, order=3)
for arr in self._moment_data_cutout])
@lazyproperty
def moments_central(self):
"""
Central moments (translation invariant) of the source up to 3rd
order.
"""
cutout_centroid = self._array('cutout_centroid')
gather = self._fast_gather
if gather is not None:
values, lx, ly, starts, counts = gather[:5]
cen_x = np.ascontiguousarray(cutout_centroid[:, 0])
cen_y = np.ascontiguousarray(cutout_centroid[:, 1])
mom = batch_moments(values, lx, ly, starts, counts, cen_x, cen_y)
# Empty sources (no overlap or fully masked) have a NaN
# centroid, so their central moments are NaN (matching the
# mask-based path).
mom[counts == 0] = np.nan
return mom
return np.array([_image_moments(arr, center=(xcen_, ycen_), order=3)
for arr, xcen_, ycen_ in
zip(self._moment_data_cutout, cutout_centroid[:, 0],
cutout_centroid[:, 1], strict=True)])
@lazyproperty
def cutout_centroid(self):
"""
The ``(x, y)`` coordinate, relative to the cutout data, of the
centroid within the aperture.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
moments = self._array('moments')
# Ignore divide-by-zero RuntimeWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
y_centroid = moments[:, 1, 0] / moments[:, 0, 0]
x_centroid = moments[:, 0, 1] / moments[:, 0, 0]
return np.transpose((x_centroid, y_centroid))
@lazyproperty
def centroid(self):
"""
The ``(x, y)`` coordinate of the centroid.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
origin = np.transpose((self.bbox_xmin, self.bbox_ymin))
return self.cutout_centroid + origin
@lazyproperty
def x_centroid(self):
"""
The ``x`` coordinate of the centroid.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
return np.transpose(self._array('centroid'))[0]
@lazyproperty
def y_centroid(self):
"""
The ``y`` coordinate of the centroid.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
return np.transpose(self._array('centroid'))[1]
@lazyproperty
def sky_centroid(self):
"""
The sky coordinate of the centroid of the unmasked pixels within
the aperture, returned as a `~astropy.coordinates.SkyCoord`
object.
The output coordinate frame is the same as the input ``wcs``.
`None` if ``wcs`` is not input.
"""
if self._wcs is None:
return self._null_object
return self._wcs.pixel_to_world(self.x_centroid, self.y_centroid)
@lazyproperty
def sky_centroid_icrs(self):
"""
The sky coordinate in the International Celestial Reference
System (ICRS) frame of the centroid of the unmasked pixels
within the aperture, returned as a
`~astropy.coordinates.SkyCoord` object.
`None` if ``wcs`` is not input.
"""
if self._wcs is None:
return self._null_object
return self.sky_centroid.icrs
@lazyproperty
def _bbox(self):
"""
The `~photutils.aperture.BoundingBox` of the aperture, always as
an iterable.
"""
apertures = self._pixel_aperture
if self.isscalar:
apertures = (apertures,)
return [aperture.bbox for aperture in apertures]
@lazyproperty
def bbox(self):
"""
The `~photutils.aperture.BoundingBox` of the aperture.
Note that the aperture bounding box is calculated using the
exact size of the aperture, which may be slightly larger than
the aperture mask calculated using the "center" method.
"""
return self._bbox
@lazyproperty
def _bbox_bounds(self):
"""
The bounding box x and y minimum and maximum bounds.
"""
bbox = self._array('bbox')
return np.array([(bbox_.ixmin, bbox_.ixmax - 1,
bbox_.iymin, bbox_.iymax - 1)
for bbox_ in bbox])
@lazyproperty
def bbox_xmin(self):
"""
The minimum ``x``-pixel index of the bounding box.
"""
return np.transpose(self._bbox_bounds)[0]
@lazyproperty
def bbox_xmax(self):
"""
The maximum ``x``-pixel index of the bounding box.
Note that this value is inclusive, unlike numpy slice indices.
"""
return np.transpose(self._bbox_bounds)[1]
@lazyproperty
def bbox_ymin(self):
"""
The minimum ``y``-pixel index of the bounding box.
"""
return np.transpose(self._bbox_bounds)[2]
@lazyproperty
def bbox_ymax(self):
"""
The maximum ``y``-pixel index of the bounding box.
Note that this value is inclusive, unlike numpy slice indices.
"""
return np.transpose(self._bbox_bounds)[3]
@lazyproperty
def _center_npixels(self):
"""
The number of unmasked pixels within each aperture using the
"center" mask method.
The result is a `~numpy.ndarray` of per-source pixel counts.
Sources with no unmasked pixels are set to NaN.
"""
gather = self._fast_gather
if gather is not None:
counts = gather[4]
npix = counts.astype(float)
npix[counts == 0] = np.nan
return npix
npix = np.array([np.sum(weight.filled(0.0))
for weight in self._weight_cutout_center])
npix[self._all_masked] = np.nan
return npix
@lazyproperty
def _sem(self):
"""
The standard error of the mean for each aperture.
The result is a `~numpy.ndarray` (without data units). It is
computed as ``s / sqrt(N)``, where ``s`` is the sample (``N -
1``) standard deviation and ``N`` is the number of unmasked
pixels within the aperture. Sources with fewer than two unmasked
pixels are set to NaN.
"""
stats = self._mean_var
if stats is not None:
var = stats[1]
else:
var = np.array([np.var(values)
for values in self._data_values_center])
npix = self._center_npixels
sem = np.full(self.n_apertures, np.nan)
mask = npix >= 2
# var is the population (ddof=0) variance, so var / (N - 1)
# equals the squared standard error of the mean.
sem[mask] = np.sqrt(var[mask] / (npix[mask] - 1.0))
return sem
@lazyproperty
def center_aper_area(self):
"""
The total area of the unmasked pixels within the aperture using
the "center" aperture mask method.
"""
return self._center_npixels * (u.pix**2)
@lazyproperty
def sum_aper_area(self):
"""
The total area of the unmasked pixels within the aperture using
the input ``sum_method`` aperture mask method.
"""
gather = self._fast_sum
if gather is not None:
area, overlap = gather[7].copy(), gather[8]
area[overlap & (area == 0)] = np.nan
return area << (u.pix**2)
areas = np.array([np.sum(weight.filled(0.0))
for weight in self._weight_cutout])
# NaN only when no sum-method pixel survives. The center-method
# ``_all_masked`` flag must not be used here: when ``sum_method``
# is not "center", unmasked boundary pixels can carry a nonzero
# fractional area even if every center-method pixel is masked
# (and their values then also contribute to ``sum``).
areas[areas == 0] = np.nan
return areas << (u.pix**2)
@lazyproperty
def sum(self):
r"""
The sum of the unmasked ``data`` values within the aperture.
.. math::
F = \sum_{i \in A} I_i
where :math:`F` is ``sum``, :math:`I_i` is the
background-subtracted ``data``, and :math:`A` are the unmasked
pixels in the aperture.
Non-finite pixel values (NaN and inf) are excluded
(automatically masked).
"""
gather = self._fast_sum
if gather is not None:
result, area, overlap = gather[5].copy(), gather[7], gather[8]
# No sum-method survivors -> NaN (matches the mask-based
# path, which returns NaN for an all-masked aperture).
result[overlap & (area == 0)] = np.nan
if self._data_unit is not None:
result <<= self._data_unit
return result
if self.sum_method == 'center':
return self._finalize_value_stat(None, np.sum)
data_values = self._get_values(self.data_sum_cutout)
result = np.array([np.sum(arr) for arr in data_values])
if self._data_unit is not None:
result <<= self._data_unit
return result
@lazyproperty
def sum_err(self):
r"""
The uncertainty of `sum`, propagated from the input ``error``
array.
``sum_err`` is the quadrature sum of the total errors over the
unmasked pixels within the aperture:
.. math::
\Delta F = \sqrt{\sum_{i \in A} \sigma_{\mathrm{tot}, i}^2}
where :math:`\Delta F` is the `sum_err`,
:math:`\sigma_{\mathrm{tot, i}}` are the pixel-wise total errors
(``error``), and :math:`A` are the unmasked pixels in the
aperture.
Pixel values that are masked in the input ``data``, including
any non-finite pixel values (NaN and inf) that are automatically
masked, are also masked in the error array.
"""
if self._error is None:
err = self._null_value
else:
gather = self._fast_sum
if gather is not None:
variance = gather[6].copy()
area = gather[7]
overlap = gather[8]
variance[overlap & (area == 0)] = np.nan
err = np.sqrt(variance)
else:
if self.sum_method == 'center':
variance = self._variance_cutout_center
else:
variance = self._variance_cutout
var_values = [arr.compressed() if len(arr.compressed()) > 0
else np.array([np.nan]) for arr in variance]
err = np.sqrt([np.sum(arr) for arr in var_values])
if self._data_unit is not None:
err <<= self._data_unit
return err
@lazyproperty
def min(self):
"""
The minimum of the unmasked pixel values within the aperture.
"""
fast = None if self._order_stats is None else self._order_stats[0]
return self._finalize_value_stat(fast, np.min)
@lazyproperty
def max(self):
"""
The maximum of the unmasked pixel values within the aperture.
"""
fast = None if self._order_stats is None else self._order_stats[1]
return self._finalize_value_stat(fast, np.max)
@lazyproperty
def mean(self):
"""
The mean of the unmasked pixel values within the aperture.
"""
fast = None if self._mean_var is None else self._mean_var[0]
return self._finalize_value_stat(fast, np.mean)
@lazyproperty
def mean_err(self):
r"""
The standard error of the `mean`.
``mean_err`` is the standard deviation of the sampling
distribution of the mean:
.. math::
\sigma_{\bar{x}} = \frac{s}{\sqrt{N}}
where :math:`s` is the sample standard deviation (computed with
``N - 1`` in the denominator) and :math:`N` is the number of
unmasked pixels within the aperture (`center_aper_area`).
Apertures with fewer than two unmasked pixels have an undefined
standard error and are set to NaN.
"""
result = self._sem.copy()
if self._data_unit is not None:
result <<= self._data_unit
return result
@lazyproperty
def median(self):
"""
The median of the unmasked pixel values within the aperture.
"""
fast = None if self._order_stats is None else self._order_stats[2]
return self._finalize_value_stat(fast, np.median)
@lazyproperty
def median_err(self):
r"""
The standard error of the `median`.
``median_err`` is the large-sample approximation of the standard
error of the median:
.. math::
\sigma_{\mathrm{med}} \approx \sqrt{\frac{\pi}{2}}
\ \frac{s}{\sqrt{N}}
where :math:`s` is the sample standard deviation (computed with
``N - 1`` in the denominator) and :math:`N` is the number of
unmasked pixels within the aperture (`center_aper_area`).
This approximation assumes that the pixel values are
approximately normally distributed. Apertures with fewer than
two unmasked pixels have an undefined standard error and are set
to NaN.
"""
result = np.sqrt(np.pi / 2.0) * self._sem
if self._data_unit is not None:
result <<= self._data_unit
return result
@lazyproperty
def mode(self):
"""
The mode of the unmasked pixel values within the aperture.
The mode is estimated as ``(3 * median) - (2 * mean)``.
"""
return 3.0 * self.median - 2.0 * self.mean
@lazyproperty
def _variance(self):
"""
The variance of the unmasked pixel values within each aperture
as a plain `~numpy.ndarray` (without data units).
The population (``ddof=0``) variance is computed using the fast
Cython reduction when available and otherwise the mask-based
per-source code path. When ``self.ddof`` is nonzero, the result
is rescaled by ``N / (N - ddof)``, where ``N`` is the number of
unmasked pixels within the aperture. Apertures with ``N <=
ddof`` unmasked pixels are set to NaN.
"""
fast = None if self._mean_var is None else self._mean_var[1]
var = self._finalize_value_stat(fast, np.var, apply_unit=False)
if self.ddof == 0:
return var
npix = self._center_npixels
result = np.full(self.n_apertures, np.nan)
mask = npix > self.ddof
result[mask] = (var[mask] * npix[mask]
/ (npix[mask] - self.ddof))
return result
@lazyproperty
def std(self):
"""
The standard deviation of the unmasked pixel values within the
aperture.
The divisor used in the calculation is ``N - ddof``, where ``N``
is the number of unmasked pixels within the aperture and ``ddof``
is the value of the ``ddof`` keyword (default 0).
"""
result = np.sqrt(self._variance)
if self._data_unit is not None:
result <<= self._data_unit
return result
@lazyproperty
def mad_std(self):
r"""
The standard deviation calculated using
the `median absolute deviation (MAD)
<https://en.wikipedia.org/wiki/Median_absolute_deviation>`_.
The standard deviation estimator is given by:
.. math::
\sigma \approx \frac{\textrm{MAD}}{\Phi^{-1}(3/4)}
\approx 1.4826 \ \textrm{MAD}
where :math:`\Phi^{-1}(P)` is the normal inverse cumulative
distribution function evaluated at probability :math:`P = 3/4`.
"""
fast = None if self._mad is None else self._mad * _MAD_STD_SCALE
return self._finalize_value_stat(fast, mad_std)
@lazyproperty
def var(self):
"""
The variance of the unmasked pixel values within the aperture.
The divisor used in the calculation is ``N - ddof``, where ``N``
is the number of unmasked pixels within the aperture and ``ddof``
is the value of the ``ddof`` keyword (default 0).
"""
result = self._variance.copy()
if self._data_unit is not None:
result <<= self._data_unit**2
return result
@lazyproperty
def biweight_location(self):
"""
The biweight location of the unmasked pixel values within the
aperture.
The tuning constant is fixed at ``c=6``, the default value
used by `astropy.stats.biweight_location`.
"""
fast = None if self._biweight is None else self._biweight[0]
return self._finalize_value_stat(fast, biweight_location)
@lazyproperty
def biweight_midvariance(self):
"""
The biweight midvariance of the unmasked pixel values within the
aperture.
The tuning constant is fixed at ``c=9``, the default value
used by `astropy.stats.biweight_midvariance`.
"""
fast = None if self._biweight is None else self._biweight[1]
return self._finalize_value_stat(fast, biweight_midvariance,
square_unit=True)
@lazyproperty
def inertia_tensor(self):
"""
The inertia tensor of the source for the rotation around its
center of mass.
"""
moments = self._array('moments_central')
mu_02 = moments[:, 0, 2]
mu_11 = -moments[:, 1, 1]
mu_20 = moments[:, 2, 0]
tensor = np.array([mu_02, mu_11, mu_11, mu_20]).swapaxes(0, 1)
return tensor.reshape((tensor.shape[0], 2, 2)) * u.pix**2
@lazyproperty
def _raw_covariance(self):
"""
The raw ``(N, 2, 2)`` covariance matrix of the 2D Gaussian
function that has the same normalized second-order moments as
the source, before any regularization.
This unregularized matrix is shared by `_covariance` (which
regularizes a copy) and `_singular_covariance_mask` (which tests
it for singularity). Callers that modify the matrix in place
must operate on a copy so the cached value is not corrupted.
"""
moments = self._array('moments_central')
# Ignore divide-by-zero RuntimeWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
mu_norm = moments / moments[:, 0, 0][:, np.newaxis, np.newaxis]
covar = np.array([mu_norm[:, 0, 2], mu_norm[:, 1, 1],
mu_norm[:, 1, 1], mu_norm[:, 2, 0]]).swapaxes(0, 1)
return covar.reshape((covar.shape[0], 2, 2))
@lazyproperty
def _covariance(self):
"""
The covariance matrix of the 2D Gaussian function that has the
same second-order moments as the source, always as an iterable.
"""
# Copy so the regularization below does not mutate the cached
# raw covariance shared with `_singular_covariance_mask`.
covar = self._raw_covariance.copy()
# Regularize the covariance matrix for "infinitely" thin
# detections by incrementally increasing the diagonal elements
# by 1/12, the variance of a uniform distribution across a
# single pixel (the smallest second moment a resolved source can
# have given finite pixel size).
delta = 1.0 / 12
delta2 = delta**2
# Ignore RuntimeWarning from NaN values in covar
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
covar_det = np.linalg.det(covar)
covar_trace = covar[:, 0, 0] + covar[:, 1, 1]
# A valid covariance is positive semidefinite (det >= 0
# and trace >= 0). Any matrix that is not (e.g., from
# net-negative flux weighting) has an undefined shape and is
# set to NaN.
bad = (covar_det < 0) | (covar_trace < 0)
covar[bad] = np.nan
# Regularize "infinitely" thin detections by adding 1/12
# (delta) to each diagonal. A single bump is sufficient.
# For a positive semidefinite matrix the bumped determinant
# exceeds the raw determinant by delta times the trace plus
# delta squared. Since the raw determinant and trace are
# both non-negative, the result is at least delta squared,
# which equals the delta2 threshold, so it clears the
# threshold in a single step.
idx = np.where(covar_det < delta2)[0]
covar[idx, 0, 0] += delta
covar[idx, 1, 1] += delta
return covar
@lazyproperty
def covariance(self):
"""
The covariance matrix of the 2D Gaussian function that has the
same second-order moments as the source.
"""
return self._covariance * (u.pix**2)
@lazyproperty
def covariance_eigvals(self):
"""
The two eigenvalues of the `covariance` matrix in decreasing
order.
"""
eigvals = np.empty((self.n_apertures, 2))
eigvals.fill(np.nan)
# np.linalg.eigvalsh requires finite input values
idx = np.unique(np.where(np.isfinite(self._covariance))[0])
eigvals[idx] = np.linalg.eigvalsh(self._covariance[idx])
# Check for negative variance
# (just in case covariance matrix is not positive semidefinite)
idx2 = np.unique(np.where(eigvals < 0)[0])
eigvals[idx2] = (np.nan, np.nan)
# Sort each eigenvalue pair in descending order
# (eigvalsh returns values in ascending order)
eigvals = np.fliplr(eigvals)
return eigvals * u.pix**2
@lazyproperty
def semimajor_axis(self):
"""
The 1-sigma standard deviation along the semimajor axis of the
2D Gaussian function that has the same second-order central
moments as the source.
"""
eigvals = self._array('covariance_eigvals')
return np.sqrt(eigvals[:, 0])
@lazyproperty
def semiminor_axis(self):
"""
The 1-sigma standard deviation along the semiminor axis of the
2D Gaussian function that has the same second-order central
moments as the source.
"""
eigvals = self._array('covariance_eigvals')
return np.sqrt(eigvals[:, 1])
@lazyproperty
def fwhm(self):
r"""
The circularized full width at half maximum (FWHM) of the 2D
Gaussian function that has the same second-order central moments
as the source.
.. math::
\mathrm{FWHM} & = 2 \sqrt{2 \ln(2)} \sqrt{0.5 (a^2 + b^2)}
\\
& = 2 \sqrt{\ln(2) \ (a^2 + b^2)}
where :math:`a` and :math:`b` are the 1-sigma lengths of the
semimajor (`semimajor_axis`) and semiminor (`semiminor_axis`)
axes, respectively.
"""
return 2.0 * np.sqrt(np.log(2.0) * (self.semimajor_axis**2
+ self.semiminor_axis**2))
@lazyproperty
def orientation(self):
"""
The angle between the ``x`` axis and the major axis of the 2D
Gaussian function that has the same second-order moments as the
source.
The angle increases in the counter-clockwise direction and is
in the range (-90, 90] degrees.
"""
covar = self._covariance
orient_radians = 0.5 * np.arctan2(2.0 * covar[:, 0, 1],
(covar[:, 0, 0] - covar[:, 1, 1]))
return np.rad2deg(orient_radians) * u.deg
@lazyproperty
def eccentricity(self):
r"""
The eccentricity of the 2D Gaussian function that has the same
second-order moments as the source.
The eccentricity is the fraction of the distance along the
semimajor axis at which the focus lies.
.. math::
e = \sqrt{1 - \frac{b^2}{a^2}}
where :math:`a` and :math:`b` are the lengths of the semimajor
and semiminor axes, respectively.
"""
semimajor_var, semiminor_var = np.transpose(self.covariance_eigvals)
return np.sqrt(1.0 - (semiminor_var / semimajor_var))
@lazyproperty
def elongation(self):
r"""
The ratio of the lengths of the semimajor and semiminor axes.
.. math::
\mathrm{elongation} = \frac{a}{b}
where :math:`a` and :math:`b` are the lengths of the semimajor
and semiminor axes, respectively.
"""
return self.semimajor_axis / self.semiminor_axis
@lazyproperty
def ellipticity(self):
r"""
1.0 minus the ratio of the lengths of the semimajor and
semiminor axes (or 1.0 minus the `elongation`).
.. math::
\mathrm{ellipticity} = 1 - \frac{b}{a}
where :math:`a` and :math:`b` are the lengths of the semimajor
and semiminor axes, respectively.
"""
return 1.0 - (self.semiminor_axis / self.semimajor_axis)
@lazyproperty
def covariance_xx(self):
r"""
The ``(0, 0)`` element of the `covariance` matrix, representing
:math:`\sigma_x^2`, in units of pixel**2.
"""
return self._covariance[:, 0, 0] * u.pix**2
@lazyproperty
def covariance_yy(self):
r"""
The ``(1, 1)`` element of the `covariance` matrix, representing
:math:`\sigma_y^2`, in units of pixel**2.
"""
return self._covariance[:, 1, 1] * u.pix**2
@lazyproperty
def covariance_xy(self):
r"""
The ``(0, 1)`` and ``(1, 0)`` elements of the `covariance`
matrix, representing :math:`\sigma_x \sigma_y`, in units of
pixel**2.
"""
return self._covariance[:, 0, 1] * u.pix**2
@lazyproperty
def ellipse_cxx(self):
r"""
Coefficient for ``x**2`` in the generalized ellipse equation in
units of pixel**(-2).
The ellipse is defined as
.. math::
cxx (x - \bar{x})^2 + cxy (x - \bar{x}) (y - \bar{y}) +
cyy (y - \bar{y})^2 = R^2
where :math:`R` is a parameter which scales the ellipse (in
units of the axes lengths).
The isophotal limit of a source is well represented by :math:`R
\approx 3`.
"""
return ((np.cos(self.orientation) / self.semimajor_axis)**2
+ (np.sin(self.orientation) / self.semiminor_axis)**2)
@lazyproperty
def ellipse_cyy(self):
r"""
Coefficient for ``y**2`` in the generalized ellipse equation in
units of pixel**(-2).
The ellipse is defined as
.. math::
cxx (x - \bar{x})^2 + cxy (x - \bar{x}) (y - \bar{y}) +
cyy (y - \bar{y})^2 = R^2
where :math:`R` is a parameter which scales the ellipse (in
units of the axes lengths).
The isophotal limit of a source is well represented by
:math:`R \approx 3`.
"""
return ((np.sin(self.orientation) / self.semimajor_axis)**2
+ (np.cos(self.orientation) / self.semiminor_axis)**2)
@lazyproperty
def ellipse_cxy(self):
r"""
Coefficient for ``x * y`` in the generalized ellipse equation in
units of pixel**(-2).
The ellipse is defined as
.. math::
cxx (x - \bar{x})^2 + cxy (x - \bar{x}) (y - \bar{y}) +
cyy (y - \bar{y})^2 = R^2
where :math:`R` is a parameter which scales the ellipse (in
units of the axes lengths).
The isophotal limit of a source is well represented by :math:`R
\approx 3`.
"""
return (2.0 * np.cos(self.orientation) * np.sin(self.orientation)
* ((1.0 / self.semimajor_axis**2)
- (1.0 / self.semiminor_axis**2)))
@lazyproperty
def gini(self):
r"""
The `Gini coefficient
<https://en.wikipedia.org/wiki/Gini_coefficient>`_ of the
unmasked pixel values within the aperture.
The Gini coefficient of the distribution of absolute flux values
is calculated using the prescription from `Lotz et al. 2004
<https://ui.adsabs.harvard.edu/abs/2004AJ....128..163L/abstract>`_
(Eq. 6) as:
.. math::
G = \frac{1}{\overline{|x|} \, n \, (n - 1)}
\sum^{n}_{i} (2i - n - 1) \left | x_i \right |
where :math:`\overline{|x|}` is the mean of the absolute value
of all pixel values :math:`x_i`. If the sum of all pixel values
is zero, the Gini coefficient is zero.
Negative pixel values are used via their absolute value. Invalid
values (NaN and inf) in the input are automatically excluded
from the calculation. If only a single finite pixel remains
after filtering, the Gini coefficient is 0.0.
"""
return self._finalize_value_stat(self._gini, gini_func,
apply_unit=False)
