# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
This module provides tools for calculating the properties of sources
defined by an Aperture.
"""
import functools
import inspect
import warnings
from copy import deepcopy
import astropy.units as u
import numpy as np
from astropy.nddata import NDData, StdDevUncertainty
from astropy.stats import (SigmaClip, biweight_location, biweight_midvariance,
mad_std)
from astropy.table import QTable
from astropy.utils import lazyproperty
from astropy.utils.exceptions import AstropyUserWarning
from photutils.aperture import Aperture, SkyAperture
from photutils.utils._misc import _get_meta
from photutils.utils._moments import _moments, _moments_central
from photutils.utils._quantity_helpers import process_quantities
__all__ = ['ApertureStats']
# default table columns for `to_table()` output
DEFAULT_COLUMNS = ['id', 'xcentroid', 'ycentroid', 'sky_centroid',
'sum', 'sum_err', 'sum_aper_area', 'center_aper_area',
'min', 'max', 'mean', 'median', 'mode', 'std',
'mad_std', 'var', 'biweight_location',
'biweight_midvariance', 'fwhm', 'semimajor_sigma',
'semiminor_sigma', 'orientation', 'eccentricity']
def as_scalar(method):
"""
Return a scalar value from a method if the class is scalar.
"""
@functools.wraps(method)
def _decorator(*args, **kwargs):
result = method(*args, **kwargs)
try:
return (result[0] if args[0].isscalar and len(result) == 1
else result)
except TypeError: # if result has no len
return result
return _decorator
[docs]class ApertureStats:
"""
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`
The aperture to apply to the data. The aperture object
may contain more than one position. If ``aperture`` is a
`~photutils.aperture.SkyAperture` object, then a WCS must be
input using the ``wcs`` keyword.
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`. 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 overlap of the aperture on
the pixel grid. This method is used only for calculating the
``sum``, ``sum_error``, ``sum_aper_area``, ``data_sumcutout``,
and ``error_sumcutout`` properties. All other properties use the
"center" aperture mask method. Not all options are available for
all aperture types. The following methods are available:
* ``'exact'`` (default):
The the exact fractional overlap of the aperture and each
pixel is calculated. The aperture weights will contain
values between 0 and 1.
* ``'center'``:
A pixel is considered to be entirely in or out of the
aperture depending on whether its center is in or out of the
aperture. The aperture weights will contain values only of 0
(out) and 1 (in).
* ``'subpixel'``:
A pixel is divided into subpixels (see the ``subpixels``
keyword), each of which are considered to be entirely in
or out of the aperture depending on whether its center is
in or out of the aperture. If ``subpixels=1``, this method
is equivalent to ``'center'``. The aperture weights will
contain values between 0 and 1.
subpixels : int, optional
For the ``'subpixel'`` method, resample pixels by this factor
in each dimension. That is, each pixel is divided into
``subpixels**2`` subpixels. This keyword is ignored unless
``sum_method='subpixel'``.
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 any 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 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.
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.
Most source properties are calculated using the "center"
aperture-mask method, which gives aperture weights of 0 or 1. This
avoids the need to compute weighted statistics --- the ``data``
pixel values are directly used.
The input ``sum_method`` and ``subpixels`` keywords are used
to determine the aperture-mask method when calculating the
sum-related properties: ``sum``, ``sum_error``, ``sum_aper_area``,
``data_sumcutout``, and ``error_sumcutout``. The default is
``sum_method='exact'``, which produces exact aperture-weighted
photometry.
.. _SourceExtractor: https://sextractor.readthedocs.io/en/latest/
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.xcentroid) # doctest: +FLOAT_CMP
149.98737072209013
>>> print(aperstats.ycentroid) # doctest: +FLOAT_CMP
24.99729176183652
>>> print(aperstats.centroid) # doctest: +FLOAT_CMP
[149.98737072 24.99729176]
>>> print(aperstats.mean, aperstats.median, aperstats.std) # doctest: +FLOAT_CMP
46.861845146453526 33.743501730319 38.25291812758177
>>> print(aperstats.sum) # doctest: +FLOAT_CMP
9118.129697119366
>>> print(aperstats.sum_aper_area) # doctest: +FLOAT_CMP
201.0619298297468 pix2
>>> # more than one aperture position
>>> aper2 = CircularAperture(((150, 25), (90, 60)), 10)
>>> aperstats2 = ApertureStats(data, aper2)
>>> print(aperstats2.xcentroid) # doctest: +FLOAT_CMP
[149.97230436 90.00833613]
>>> print(aperstats2.sum) # doctest: +FLOAT_CMP
[ 9863.56195844 36629.52906175]
"""
def __init__(self, data, aperture, *, error=None, mask=None, wcs=None,
sigma_clip=None, sum_method='exact', subpixels=5,
local_bkg=None):
if isinstance(data, NDData):
data, error, mask, wcs = self._unpack_nddata(data, error, mask,
wcs)
(data, error, local_bkg), unit = process_quantities(
(data, error, local_bkg), ('data', 'error', 'local_bkg'))
self._data = self._validate_array(data, 'data', shape=False)
self._data_unit = unit
self.aperture = self._validate_aperture(aperture)
if isinstance(aperture, SkyAperture):
if wcs is None:
raise ValueError('A wcs is required when using a SkyAperture')
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):
raise TypeError('sigma_clip must be a SigmaClip instance')
self.sigma_clip = sigma_clip
self.sum_method = sum_method
self.subpixels = subpixels
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:
raise ValueError('local_bkg must be a 1D array')
n_local_bkg = len(local_bkg)
if n_local_bkg != 1 and n_local_bkg != self.n_apertures:
raise ValueError('local_bkg must be scalar or have the same '
'length as the input aperture')
local_bkg = np.broadcast_to(local_bkg, self.n_apertures)
if np.any(~np.isfinite(local_bkg)):
raise ValueError('local_bkg must not contain any non-finite '
'(e.g., inf or NaN) values')
self._local_bkg = local_bkg # always an iterable
self._ids = np.arange(self.n_apertures) + 1
self.default_columns = DEFAULT_COLUMNS
self.meta = _get_meta()
@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:
warnings.warn(f'The {key!r} keyword is be ignored. Its '
'value is obtained from the input NDData '
'object.', 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):
if not isinstance(aperture, Aperture):
raise TypeError('aperture must be an Aperture object')
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:
raise ValueError(f'{name} must be a {ndim}D array.')
if shape and array.shape != self._data.shape:
raise ValueError(f'data and {name} must have the same shape.')
return array
@property
def _lazyproperties(self):
"""
A list of all class lazyproperties (even in superclasses).
"""
def islazyproperty(obj):
return isinstance(obj, lazyproperty)
return [i[0] for i in inspect.getmembers(self.__class__,
predicate=islazyproperty)]
@property
def properties(self):
"""
A sorted list of built-in source properties.
"""
lazyproperties = [name for name in self._lazyproperties if not
name.startswith('_')]
lazyproperties.sort()
return lazyproperties
def __getitem__(self, index):
if self.isscalar:
raise TypeError(f'A scalar {self.__class__.__name__!r} object '
'cannot be indexed')
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',
'default_columns', 'meta')
for attr in init_attr:
setattr(newcls, attr, getattr(self, attr))
# need to slice _aperture and _ids;
# aperture determines isscalar (needed below)
attrs = ('aperture', '_ids')
for attr in attrs:
setattr(newcls, attr, getattr(self, attr)[index])
# slice evaluated lazyproperty objects
keys = set(self.__dict__.keys()) & set(self._lazyproperties)
keys.add('_local_bkg') # iterable defined in __init__
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:
raise TypeError(f'Scalar {self.__class__.__name__!r} object has '
'no len()')
return self.n_apertures
def __iter__(self):
for item in range(len(self)):
yield self.__getitem__(item)
@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.
"""
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
@as_scalar
def id(self):
"""
The aperture identification number(s).
"""
return self._ids
@property
def ids(self):
"""
The aperture identification number(s), always as an iterable
`~numpy.ndarray`.
"""
_ids = self._ids
if self.isscalar:
_ids = np.array((_ids,))
return _ids
[docs] def get_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.get_ids(id_num)
[docs] def get_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.ids:
raise ValueError(f'{id_num} is not a valid source ID number')
sorter = np.argsort(self.id)
indices = sorter[np.searchsorted(self.id, id_nums, sorter=sorter)]
return self[indices]
[docs] 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.
"""
if columns is None:
table_columns = self.default_columns
else:
table_columns = np.atleast_1d(columns)
tbl = QTable(meta=self.meta)
for column in table_columns:
values = getattr(self, column)
# column assignment requires an object with a length
if self.isscalar:
values = (values,)
tbl[column] = values
return tbl
@lazyproperty
def n_apertures(self):
"""
The number of positions in 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 _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 a iterable.
"""
cutouts = []
for (slices, local_bkg) in zip(self._overlap_slices,
self._local_bkg):
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):
"""
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.
Returns
-------
data, variance, mask, weights : list of `~numpy.ndarray`
A list of cutout arrays for the data, variance, mask and weight
arrays for each source (aperture position).
"""
data_cutouts = []
variance_cutouts = []
mask_cutouts = []
weight_cutouts = []
overlaps = []
for (data_cutout, apermask, slices) in zip(
self._data_cutouts, aperture_masks, self._overlap_slices):
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.
data_mask = ~np.isfinite(data_cutout)
if self._mask is not None:
data_mask |= self._mask[slc_large]
overlap = True
aperweight_cutout = apermask.data[slc_small]
weight_cutout = aperweight_cutout * ~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
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 = self._error[slc_large]**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)
# 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))
@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.
"""
return self._make_aperture_cutouts(self._aperture_masks)
@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))[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))[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)]
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)]
@lazyproperty
@as_scalar
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))[0])
@lazyproperty
@as_scalar
def data_sumcutout(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))[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))[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))[1])
@lazyproperty
@as_scalar
def error_sumcutout(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 of 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))[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 of 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))[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.data_cutout) # self.data_cutout is a list
if self.isscalar:
data = (data,)
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))[4]
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
@as_scalar
def moments(self):
"""
Spatial moments up to 3rd order of the source.
"""
return np.array([_moments(arr, order=3) for arr in
self._moment_data_cutout])
@lazyproperty
@as_scalar
def moments_central(self):
"""
Central moments (translation invariant) of the source up to 3rd
order.
"""
cutout_centroid = self.cutout_centroid
if self.isscalar:
cutout_centroid = cutout_centroid[np.newaxis, :]
return np.array([_moments_central(arr, center=(xcen_, ycen_), order=3)
for arr, xcen_, ycen_ in
zip(self._moment_data_cutout, cutout_centroid[:, 0],
cutout_centroid[:, 1])])
@lazyproperty
@as_scalar
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.moments
if self.isscalar:
moments = moments[np.newaxis, :]
# ignore divide-by-zero RuntimeWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
ycentroid = moments[:, 1, 0] / moments[:, 0, 0]
xcentroid = moments[:, 0, 1] / moments[:, 0, 0]
return np.transpose((xcentroid, ycentroid))
@lazyproperty
@as_scalar
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 _xcentroid(self):
"""
The ``x`` coordinate of the centroid, always as an iterable.
"""
xcentroid = np.transpose(self.centroid)[0]
if self.isscalar:
xcentroid = (xcentroid,)
return xcentroid
@lazyproperty
@as_scalar
def xcentroid(self):
"""
The ``x`` coordinate of the centroid.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
return self._xcentroid
@lazyproperty
def _ycentroid(self):
"""
The ``y`` coordinate of the centroid, always as an iterable.
"""
ycentroid = np.transpose(self.centroid)[1]
if self.isscalar:
ycentroid = (ycentroid,)
return ycentroid
@lazyproperty
@as_scalar
def ycentroid(self):
"""
The ``y`` coordinate of the centroid.
The centroid is computed as the center of mass of the unmasked
pixels within the aperture.
"""
return self._ycentroid
@lazyproperty
@as_scalar
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.xcentroid, self.ycentroid)
@lazyproperty
@as_scalar
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
@as_scalar
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" mode.
"""
return self._bbox
@lazyproperty
@as_scalar
def _bbox_bounds(self):
"""
The bounding box x/y minimum and maximum bounds.
"""
bbox = self.bbox
if self.isscalar:
bbox = (bbox,)
return np.array([(bbox_.ixmin, bbox_.ixmax - 1,
bbox_.iymin, bbox_.iymax - 1)
for bbox_ in bbox])
@lazyproperty
@as_scalar
def bbox_xmin(self):
"""
The minimum ``x`` pixel index of the bounding box.
"""
return np.transpose(self._bbox_bounds)[0]
@lazyproperty
@as_scalar
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
@as_scalar
def bbox_ymin(self):
"""
The minimum ``y`` pixel index of the bounding box.
"""
return np.transpose(self._bbox_bounds)[2]
@lazyproperty
@as_scalar
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]
def _calculate_stats(self, stat_func, unit=None):
"""
Apply the input ``stat_func`` to the 1D array of unmasked data
values in the aperture.
Units are applied if the input ``data`` has units.
Parameters
----------
stat_func : callable
The callable to apply to the 1D `~numpy.ndarray` of unmasked
data values.
unit : `None` or `astropy.unit.Unit`, optional
The unit to apply to the output data. This is used only
if the input ``data`` has units. If `None` then the input
``data`` unit will be used.
"""
result = np.array([stat_func(arr) for arr in self._data_values_center])
if unit is None:
unit = self._data_unit
if unit is not None:
result <<= unit
return result
@lazyproperty
@as_scalar
def center_aper_area(self):
"""
The total area of the unmasked pixels within the aperture using
the "center" aperture mask method.
"""
areas = np.array([np.sum(weight.filled(0.0))
for weight in self._weight_cutout_center])
areas[self._all_masked] = np.nan
return areas << (u.pix**2)
@lazyproperty
@as_scalar
def sum_aper_area(self):
"""
The total area of the unmasked pixels within the aperture using
the input ``sum_method`` aperture mask method.
"""
areas = np.array([np.sum(weight.filled(0.0))
for weight in self._weight_cutout])
areas[self._all_masked] = np.nan
return areas << (u.pix**2)
@lazyproperty
@as_scalar
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).
"""
if self.sum_method == 'center':
return self._calculate_stats(np.sum)
data_values = self._get_values(self.data_sumcutout)
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
@as_scalar
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`, :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:
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
@as_scalar
def min(self):
"""
The minimum of the unmasked pixel values within the aperture.
"""
return self._calculate_stats(np.min)
@lazyproperty
@as_scalar
def max(self):
"""
The maximum of the unmasked pixel values within the aperture.
"""
return self._calculate_stats(np.max)
@lazyproperty
@as_scalar
def mean(self):
"""
The mean of the unmasked pixel values within the aperture.
"""
return self._calculate_stats(np.mean)
@lazyproperty
@as_scalar
def median(self):
"""
The median of the unmasked pixel values within the aperture.
"""
return self._calculate_stats(np.median)
@lazyproperty
@as_scalar
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
@as_scalar
def std(self):
"""
The standard deviation of the unmasked pixel values within the
aperture.
"""
return self._calculate_stats(np.std)
@lazyproperty
@as_scalar
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`.
"""
return self._calculate_stats(mad_std)
@lazyproperty
@as_scalar
def var(self):
"""
The variance of the unmasked pixel values within the aperture.
"""
unit = self._data_unit
if unit is not None:
unit **= 2
return self._calculate_stats(np.var, unit=unit)
@lazyproperty
@as_scalar
def biweight_location(self):
"""
The biweight location of the unmasked pixel values within the
aperture.
See `astropy.stats.biweight_location`.
"""
return self._calculate_stats(biweight_location)
@lazyproperty
@as_scalar
def biweight_midvariance(self):
"""
The biweight midvariance of the unmasked pixel values within the
aperture.
See `astropy.stats.biweight_midvariance`
"""
unit = self._data_unit
if unit is not None:
unit **= 2
return self._calculate_stats(biweight_midvariance, unit=unit)
@lazyproperty
@as_scalar
def inertia_tensor(self):
"""
The inertia tensor of the source for the rotation around its
center of mass.
"""
moments = self.moments_central
if self.isscalar:
moments = moments[np.newaxis, :]
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 _covariance(self):
"""
The covariance matrix of the 2D Gaussian function that has the
same second-order moments as the source, always as an iterable.
"""
moments = self.moments_central
if self.isscalar:
moments = moments[np.newaxis, :]
# 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)
covar = covar.reshape((covar.shape[0], 2, 2))
# Modify the covariance matrix in the case of "infinitely" thin
# detections. This follows SourceExtractor's prescription of
# incrementally increasing the diagonal elements by 1/12.
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)
idx = np.where(covar_det < delta2)[0]
while idx.size > 0: # pragma: no cover
covar[idx, 0, 0] += delta
covar[idx, 1, 1] += delta
covar_det = np.linalg.det(covar)
idx = np.where(covar_det < delta2)[0]
return covar
@lazyproperty
@as_scalar
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
@as_scalar
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.eivals requires finite input values
idx = np.unique(np.where(np.isfinite(self._covariance))[0])
eigvals[idx] = np.linalg.eigvals(self._covariance[idx])
# check for negative variance
# (just in case covariance matrix is not positive (semi)definite)
idx2 = np.unique(np.where(eigvals < 0)[0]) # pragma: no cover
eigvals[idx2] = (np.nan, np.nan) # pragma: no cover
# sort each eigenvalue pair in descending order
eigvals.sort(axis=1)
eigvals = np.fliplr(eigvals)
return eigvals * u.pix**2
@lazyproperty
@as_scalar
def semimajor_sigma(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.covariance_eigvals
if self.isscalar:
eigvals = eigvals[np.newaxis, :]
# this matches SourceExtractor's A parameter
return np.sqrt(eigvals[:, 0])
@lazyproperty
@as_scalar
def semiminor_sigma(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.covariance_eigvals
if self.isscalar:
eigvals = eigvals[np.newaxis, :]
# this matches SourceExtractor's B parameter
return np.sqrt(eigvals[:, 1])
@lazyproperty
@as_scalar
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_sigma`) and semiminor (`semiminor_sigma`)
axes, respectively.
"""
return 2.0 * np.sqrt(np.log(2.0) * (self.semimajor_sigma**2
+ self.semiminor_sigma**2))
@lazyproperty
@as_scalar
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.
"""
covar = self._covariance
orient_radians = 0.5 * np.arctan2(2.0 * covar[:, 0, 1],
(covar[:, 0, 0] - covar[:, 1, 1]))
return orient_radians * 180.0 / np.pi * u.deg
@lazyproperty
@as_scalar
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
@as_scalar
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_sigma / self.semiminor_sigma
@lazyproperty
@as_scalar
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_sigma / self.semimajor_sigma)
@lazyproperty
@as_scalar
def covar_sigx2(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
@as_scalar
def covar_sigy2(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
@as_scalar
def covar_sigxy(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
@as_scalar
def cxx(self):
r"""
`SourceExtractor`_'s CXX ellipse parameter 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). `SourceExtractor`_ reports that the
isophotal limit of a source is well represented by :math:`R
\approx 3`.
"""
return ((np.cos(self.orientation) / self.semimajor_sigma)**2
+ (np.sin(self.orientation) / self.semiminor_sigma)**2)
@lazyproperty
@as_scalar
def cyy(self):
r"""
`SourceExtractor`_'s CYY ellipse parameter 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). `SourceExtractor`_ reports that the
isophotal limit of a source is well represented by :math:`R
\approx 3`.
"""
return ((np.sin(self.orientation) / self.semimajor_sigma)**2
+ (np.cos(self.orientation) / self.semiminor_sigma)**2)
@lazyproperty
@as_scalar
def cxy(self):
r"""
`SourceExtractor`_'s CXY ellipse parameter 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). `SourceExtractor`_ reports that 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_sigma**2)
- (1.0 / self.semiminor_sigma**2)))
@lazyproperty
@as_scalar
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 is calculated using the prescription from
`Lotz et al. 2004
<https://ui.adsabs.harvard.edu/abs/2004AJ....128..163L/abstract>`_
as:
.. math::
G = \frac{1}{\left | \bar{x} \right | n (n - 1)}
\sum^{n}_{i} (2i - n - 1) \left | x_i \right |
where :math:`\bar{x}` is the mean over pixel values :math:`x_i`
within the aperture.
The Gini coefficient is a way of measuring the inequality in a
given set of values. In the context of galaxy morphology, it
measures how the light of a galaxy image is distributed among
its pixels. A Gini coefficient value of 0 corresponds to a
galaxy image with the light evenly distributed over all pixels
while a Gini coefficient value of 1 represents a galaxy image
with all its light concentrated in just one pixel.
"""
gini = []
for arr in self._data_values_center:
if np.all(np.isnan(arr)):
gini.append(np.nan)
continue
npix = np.size(arr)
normalization = np.abs(np.mean(arr)) * npix * (npix - 1)
kernel = ((2.0 * np.arange(1, npix + 1) - npix - 1)
* np.abs(np.sort(arr)))
gini.append(np.sum(kernel) / normalization)
return np.array(gini)