# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Tools for generating radial profiles.
"""
import warnings
import numpy as np
from astropy.modeling.fitting import TRFLSQFitter
from astropy.modeling.models import Gaussian1D, Moffat1D
from astropy.stats import gaussian_sigma_to_fwhm
from astropy.utils import lazyproperty
from astropy.utils.exceptions import AstropyUserWarning
from photutils.aperture.core import _update_method_subpixels_docstring
from photutils.profiles.core import ProfileBase
__all__ = ['RadialProfile']
[docs]
@_update_method_subpixels_docstring
class RadialProfile(ProfileBase):
"""
Class to create a radial profile using concentric circular annulus
apertures.
The radial profile represents the azimuthally-averaged flux in
circular annuli apertures as a function of radius.
For this class, the input radii represent the edges of the radial
bins. This differs from the `CurveOfGrowth` class, where the input
radii are the radii of the circular apertures used to compute the
cumulative flux.
The output `radius` attribute represents the bin centers.
Parameters
----------
data : 2D `~numpy.ndarray`
The 2D data array. The data should be background-subtracted.
Non-finite values (e.g., NaN or inf) in the ``data`` or
``error`` array are automatically masked.
xycen : tuple of 2 floats
The ``(x, y)`` pixel coordinate of the source center.
radii : 1D float `~numpy.ndarray`
An array of radii defining the *edges* of the radial bins.
``radii`` must be strictly increasing with a minimum value
greater than or equal to zero, and contain at least 2 values.
The radial spacing does not need to be constant. The output
`radius` attribute will be defined at the bin centers.
error : 2D `~numpy.ndarray`, optional
The 1-sigma errors of the input ``data``. ``error`` is assumed
to include all sources of error, including the Poisson error of
the sources (see `~photutils.utils.calc_total_error`). ``error``
must have the same shape as the input ``data``. Non-finite
values (e.g., NaN or inf) in the ``data`` or ``error`` array are
automatically masked.
mask : 2D bool `~numpy.ndarray`, optional
A boolean mask with the same shape as ``data`` where a `True`
value indicates the corresponding element of ``data`` is masked.
Masked data are excluded from all calculations.
<method_subpixels_descriptions>
See Also
--------
CurveOfGrowth
Notes
-----
If the minimum of ``radii`` is zero, then a circular aperture
with radius equal to ``radii[1]`` will be used for the innermost
aperture.
Examples
--------
>>> import numpy as np
>>> from astropy.modeling.models import Gaussian2D
>>> from photutils.centroids import centroid_2dg
>>> from photutils.datasets import make_noise_image
>>> from photutils.profiles import RadialProfile
Create an artificial single source. Note that this image does not
have any background.
>>> gmodel = Gaussian2D(42.1, 47.8, 52.4, 4.7, 4.7, 0)
>>> yy, xx = np.mgrid[0:100, 0:100]
>>> data = gmodel(xx, yy)
>>> bkg_sig = 2.1
>>> noise = make_noise_image(data.shape, mean=0., stddev=bkg_sig, seed=0)
>>> data += noise
>>> error = np.zeros_like(data) + bkg_sig
Create the radial profile.
>>> xycen = centroid_2dg(data)
>>> edge_radii = np.arange(25)
>>> rp = RadialProfile(data, xycen, edge_radii, error=error)
>>> print(rp.radius) # doctest: +FLOAT_CMP
[ 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5
14.5 15.5 16.5 17.5 18.5 19.5 20.5 21.5 22.5 23.5]
>>> print(rp.profile) # doctest: +FLOAT_CMP
[ 4.30187860e+01 4.02502046e+01 3.57758011e+01 3.16071235e+01
2.61511082e+01 2.10539746e+01 1.63701300e+01 1.16674718e+01
8.12828014e+00 5.78962699e+00 3.59342666e+00 2.35353336e+00
1.20355937e+00 7.67093923e-01 4.24650784e-01 8.67989701e-02
5.11484374e-02 -9.82041768e-02 2.37482124e-02 -3.66602855e-02
6.84802299e-02 1.72239596e-01 -3.86056497e-02 7.30423743e-02]
>>> print(rp.profile_error) # doctest: +FLOAT_CMP
[1.18479813 0.68404352 0.52985783 0.4478116 0.39493271 0.35723008
0.32860388 0.30591356 0.28735575 0.27181133 0.25854415 0.24704749
0.23695963 0.22801451 0.22001149 0.21279603 0.20624688 0.20026744
0.19477961 0.18971954 0.18503438 0.18068002 0.17661928 0.17282057]
Plot the radial profile, including the raw data profile.
.. plot::
import matplotlib.pyplot as plt
import numpy as np
from astropy.modeling.models import Gaussian2D
from photutils.centroids import centroid_2dg
from photutils.datasets import make_noise_image
from photutils.profiles import RadialProfile
# Create an artificial single source
gmodel = Gaussian2D(42.1, 47.8, 52.4, 4.7, 4.7, 0)
yy, xx = np.mgrid[0:100, 0:100]
data = gmodel(xx, yy)
bkg_sig = 2.1
noise = make_noise_image(data.shape, mean=0., stddev=bkg_sig, seed=0)
data += noise
error = np.zeros_like(data) + bkg_sig
# Find the source centroid
xycen = centroid_2dg(data)
# Create the radial profile
edge_radii = np.arange(26)
rp = RadialProfile(data, xycen, edge_radii, error=error)
# Plot the radial profile
fig, ax = plt.subplots(figsize=(8, 6))
rp.plot(ax=ax, color='C0')
rp.plot_error(ax=ax)
ax.scatter(rp.data_radius, rp.data_profile, s=1, color='C1')
Normalize the profile and plot the normalized radial profile.
.. plot::
import matplotlib.pyplot as plt
import numpy as np
from astropy.modeling.models import Gaussian2D
from photutils.centroids import centroid_2dg
from photutils.datasets import make_noise_image
from photutils.profiles import RadialProfile
# Create an artificial single source
gmodel = Gaussian2D(42.1, 47.8, 52.4, 4.7, 4.7, 0)
yy, xx = np.mgrid[0:100, 0:100]
data = gmodel(xx, yy)
bkg_sig = 2.1
noise = make_noise_image(data.shape, mean=0., stddev=bkg_sig, seed=0)
data += noise
error = np.zeros_like(data) + bkg_sig
# Find the source centroid
xycen = centroid_2dg(data)
# Create the radial profile
edge_radii = np.arange(26)
rp = RadialProfile(data, xycen, edge_radii, error=error)
# Plot the radial profile
rp.normalize()
fig, ax = plt.subplots(figsize=(8, 6))
rp.plot(ax=ax)
rp.plot_error(ax=ax)
Plot three of the annulus apertures on the data.
.. plot::
import matplotlib.pyplot as plt
import numpy as np
from astropy.modeling.models import Gaussian2D
from astropy.visualization import simple_norm
from photutils.centroids import centroid_2dg
from photutils.datasets import make_noise_image
from photutils.profiles import RadialProfile
# Create an artificial single source
gmodel = Gaussian2D(42.1, 47.8, 52.4, 4.7, 4.7, 0)
yy, xx = np.mgrid[0:100, 0:100]
data = gmodel(xx, yy)
bkg_sig = 2.1
noise = make_noise_image(data.shape, mean=0., stddev=bkg_sig, seed=0)
data += noise
error = np.zeros_like(data) + bkg_sig
# Find the source centroid
xycen = centroid_2dg(data)
# Create the radial profile
edge_radii = np.arange(26)
rp = RadialProfile(data, xycen, edge_radii, error=error)
norm = simple_norm(data, 'sqrt')
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(data, norm=norm, origin='lower')
rp.apertures[5].plot(ax=ax, color='C0', lw=2)
rp.apertures[10].plot(ax=ax, color='C1', lw=2)
rp.apertures[15].plot(ax=ax, color='C3', lw=2)
Fit 1D Gaussian and Moffat models to the normalized radial profile.
>>> rp.normalize()
>>> rp.gaussian_fit # doctest: +FLOAT_CMP
<Gaussian1D(amplitude=0.98231088, mean=0., stddev=4.67512776)>
>>> rp.moffat_fit # doctest: +ELLIPSIS
<Moffat1D(amplitude=..., x_0=0., gamma=..., alpha=...)>
>>> print(rp.gaussian_fwhm) # doctest: +FLOAT_CMP
11.009084813327846
>>> print(rp.moffat_fwhm) # doctest: +FLOAT_CMP
10.868426507928344
Plot the fitted 1D Gaussian and Moffat models on the radial profile.
.. plot::
import matplotlib.pyplot as plt
import numpy as np
from astropy.modeling.models import Gaussian2D
from photutils.centroids import centroid_2dg
from photutils.datasets import make_noise_image
from photutils.profiles import RadialProfile
# Create an artificial single source
gmodel = Gaussian2D(42.1, 47.8, 52.4, 4.7, 4.7, 0)
yy, xx = np.mgrid[0:100, 0:100]
data = gmodel(xx, yy)
bkg_sig = 2.1
noise = make_noise_image(data.shape, mean=0., stddev=bkg_sig, seed=0)
data += noise
error = np.zeros_like(data) + bkg_sig
# Find the source centroid
xycen = centroid_2dg(data)
# Create the radial profile
edge_radii = np.arange(26)
rp = RadialProfile(data, xycen, edge_radii, error=error)
# Plot the radial profile
rp.normalize()
fig, ax = plt.subplots(figsize=(8, 6))
rp.plot(ax=ax, label='Radial Profile')
rp.plot_error(ax=ax)
ax.plot(rp.radius, rp.gaussian_profile, label='Gaussian Fit')
ax.plot(rp.radius, rp.moffat_profile, label='Moffat Fit')
ax.legend()
"""
# Define y-axis label used by `~photutils.profiles.ProfileBase.plot`
_ylabel = 'Radial Profile'
# Define the fit properties that should be invalidated when the
# profile normalization is changed, so they are always consistent
# with the current profile.
_fit_properties = ('_profile_nanmask', 'gaussian_fit',
'gaussian_profile', 'gaussian_fwhm', 'moffat_fit',
'moffat_profile', 'moffat_fwhm')
@lazyproperty
def radius(self):
"""
The profile radius (bin centers) in pixels as a 1D
`~numpy.ndarray`.
The returned radius values are defined as the arithmetic means
of the input radial-bins edges (``radii``).
For logarithmically-spaced input ``radii``, one could instead
use a radius array defined using the geometric mean of the bin
edges, i.e. ``np.sqrt(radii[:-1] * radii[1:])``.
"""
# Define the radial bin centers from the radial bin edges
return (self.radii[:-1] + self.radii[1:]) / 2
@lazyproperty
def apertures(self):
"""
A list of the circular annulus apertures used to measure the
radial profile, as `~photutils.aperture.CircularAnnulus`
objects.
If the minimum of ``radii`` is zero, then the innermost element
will be a `~photutils.aperture.CircularAperture` with radius
equal to ``radii[1]``.
"""
# Prevent circular imports
from photutils.aperture import CircularAnnulus, CircularAperture
apertures = []
for i in range(len(self.radii) - 1):
try:
aperture = CircularAnnulus(self.xycen, self.radii[i],
self.radii[i + 1])
except ValueError: # zero radius
aperture = CircularAperture(self.xycen,
self.radii[i + 1])
apertures.append(aperture)
return apertures
@lazyproperty
def _flux(self):
"""
The flux in a circular annulus.
"""
return np.diff(self._photometry[0])
@lazyproperty
def _flux_err(self):
"""
The flux error in a circular annulus.
"""
return np.sqrt(np.diff(self._photometry[1] ** 2))
@lazyproperty
def area(self):
"""
The unmasked area in each circular annulus (or aperture) as a
function of radius as a 1D `~numpy.ndarray`.
"""
return np.diff(self._photometry[2])
@lazyproperty
def profile(self):
"""
The radial profile as a 1D `~numpy.ndarray`.
"""
# Ignore divide-by-zero RuntimeWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
return self._flux / self.area
@lazyproperty
def profile_error(self):
"""
The radial profile errors as a 1D `~numpy.ndarray`.
If no ``error`` array was provided, an empty array with shape
``(0,)`` is returned.
"""
if self.error is None:
return self._flux_err
# Ignore divide-by-zero RuntimeWarning
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
return self._flux_err / self.area
@lazyproperty
def _profile_nanmask(self):
return np.isfinite(self.profile)
@lazyproperty
def gaussian_fit(self):
"""
The fitted 1D Gaussian to the radial profile as a
`~astropy.modeling.functional_models.Gaussian1D` model.
The cached fit is automatically invalidated when the profile
normalization is changed, so the fit is always consistent with
the current profile.
"""
profile = self.profile[self._profile_nanmask]
radius = self.radius[self._profile_nanmask]
if len(profile) == 0:
msg = ('The radial profile is entirely non-finite or masked; '
'cannot fit a Gaussian.')
warnings.warn(msg, AstropyUserWarning)
return None
amplitude = np.max(profile)
sum_profile = np.sum(profile)
if sum_profile == 0:
msg = ('The profile sum is zero; the Gaussian fit initial '
'guess may be inaccurate.')
warnings.warn(msg, AstropyUserWarning)
std = 1.0 # fallback to avoid zero initial guess
else:
std = np.sqrt(abs(np.sum(profile * radius**2) / sum_profile))
std = max(std, 1.0) # guard against near-zero initial guess
g_init = Gaussian1D(amplitude=amplitude, mean=0.0, stddev=std)
g_init.mean.fixed = True
fitter = TRFLSQFitter()
gaussian_fit = fitter(g_init, radius, profile)
if radius.min() > 0.3 * gaussian_fit.stddev.value:
msg = ('Gaussian fit may be unreliable because the input '
'radii do not extend close to the source center.')
warnings.warn(msg, AstropyUserWarning)
return gaussian_fit
@lazyproperty
def gaussian_profile(self):
"""
The fitted 1D Gaussian profile to the radial profile as a 1D
`~numpy.ndarray`.
The cached profile is automatically invalidated when the profile
normalization is changed.
Returns `None` if the fit failed (e.g., the profile is entirely
non-finite or masked).
"""
if self.gaussian_fit is None:
return None
return self.gaussian_fit(self.radius)
@lazyproperty
def gaussian_fwhm(self):
"""
The full-width at half-maximum (FWHM) in pixels of the 1D
Gaussian fitted to the radial profile.
The cached value is automatically invalidated when the profile
normalization is changed.
Returns `None` if the fit failed (e.g., the profile is entirely
non-finite or masked).
"""
if self.gaussian_fit is None:
return None
return self.gaussian_fit.stddev.value * gaussian_sigma_to_fwhm
@lazyproperty
def moffat_fit(self):
"""
The fitted 1D Moffat to the radial profile as a
`~astropy.modeling.functional_models.Moffat1D` model.
The cached fit is automatically invalidated when the profile
normalization is changed, so the fit is always consistent with
the current profile.
"""
profile = self.profile[self._profile_nanmask]
radius = self.radius[self._profile_nanmask]
if len(profile) == 0:
msg = ('The radial profile is entirely non-finite or masked; '
'cannot fit a Moffat.')
warnings.warn(msg, AstropyUserWarning)
return None
amplitude = np.max(profile)
sum_profile = np.sum(profile)
if sum_profile == 0:
msg = ('The profile sum is zero; the Moffat fit initial '
'guess may be inaccurate.')
warnings.warn(msg, AstropyUserWarning)
gamma = 1.0 # fallback to avoid zero initial guess
else:
# Estimate gamma from the half-max radius
half_max = amplitude / 2.0
above = profile >= half_max
gamma = (max(np.max(radius[above]), 1.0)
if np.any(above) else 1.0)
m_init = Moffat1D(amplitude=amplitude, x_0=0.0, gamma=gamma,
alpha=2.5)
m_init.x_0.fixed = True
m_init.gamma.bounds = (0, None)
m_init.alpha.bounds = (1, 25)
fitter = TRFLSQFitter()
return fitter(m_init, radius, profile)
@lazyproperty
def moffat_profile(self):
"""
The fitted 1D Moffat profile to the radial profile as a 1D
`~numpy.ndarray`.
The cached profile is automatically invalidated when the profile
normalization is changed.
Returns `None` if the fit failed (e.g., the profile is entirely
non-finite or masked).
"""
if self.moffat_fit is None:
return None
return self.moffat_fit(self.radius)
@lazyproperty
def moffat_fwhm(self):
"""
The full-width at half-maximum (FWHM) in pixels of the 1D Moffat
fitted to the radial profile.
The cached value is automatically invalidated when the profile
normalization is changed.
Returns `None` if the fit failed (e.g., the profile is entirely
non-finite or masked).
"""
if self.moffat_fit is None:
return None
return self.moffat_fit.fwhm
@lazyproperty
def _data_profile(self):
"""
The raw data profile returned as 1D arrays (`~numpy.ndarray`) of
radii and data values.
Returns the radii and values of the unmasked data points within
the maximum radius defined by the input radii. Pixels flagged
in ``self.mask`` (including auto-masked non-finite values) are
excluded.
"""
shape = self.data.shape
max_radius = np.max(self.radii)
x_min = int(max(np.floor(self.xycen[0] - max_radius), 0))
x_max = int(min(np.ceil(self.xycen[0] + max_radius), shape[1]))
y_min = int(max(np.floor(self.xycen[1] - max_radius), 0))
y_max = int(min(np.ceil(self.xycen[1] + max_radius), shape[0]))
yidx, xidx = np.indices((y_max - y_min, x_max - x_min))
xidx += x_min
yidx += y_min
# Calculate the radii of the pixels from the center and select
# those within the maximum radius defined by the input radii
radii = np.hypot(xidx - self.xycen[0], yidx - self.xycen[1])
within = radii <= max_radius
radii = radii[within]
yidx_sub = yidx[within]
xidx_sub = xidx[within]
# Exclude masked pixels (user mask and auto-masked non-finite
# values)
valid = ~self.mask[yidx_sub, xidx_sub]
radii = radii[valid]
data_values = self.data[yidx_sub[valid], xidx_sub[valid]]
return radii, data_values
@lazyproperty
def data_radius(self):
"""
The radii of the raw data profile as a 1D `~numpy.ndarray`.
"""
return self._data_profile[0]
@lazyproperty
def data_profile(self):
"""
The raw data profile as a 1D `~numpy.ndarray`.
"""
return self._data_profile[1]
def _invalidate_fit_cache(self):
"""
Remove cached Gaussian and Moffat fit lazy properties so they
are recomputed on next access using the current profile.
"""
for key in self._fit_properties:
self.__dict__.pop(key, None)
def _normalize_hook(self, normalization):
"""
Also normalize ``data_profile`` if it has been computed, and
invalidate fit caches so they are recomputed on next access.
"""
if 'data_profile' in self.__dict__:
self.__dict__['data_profile'] = self.data_profile / normalization
self._invalidate_fit_cache()
def _unnormalize_hook(self):
"""
Also unnormalize ``data_profile`` if it has been computed, and
invalidate fit caches so they are recomputed on next access.
"""
if 'data_profile' in self.__dict__:
self.__dict__['data_profile'] = (self.data_profile
* self.normalization_value)
self._invalidate_fit_cache()
