-
-
Notifications
You must be signed in to change notification settings - Fork 131
/
fwhm.py
184 lines (151 loc) · 7.82 KB
/
fwhm.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
This module defines tools to perform aperture photometry.
"""
import warnings
import numpy as np
from astropy.nddata import NDData, StdDevUncertainty
import astropy.units as u
from astropy.utils.exceptions import AstropyUserWarning
from photutils.aperture import Aperture, CircularAperture, SkyCircularAperture, SkyAperture
from photutils.aperture._photometry_utils import (_handle_units, _prepare_photometry_data,
_validate_inputs)
__all__ = ['fwhm_cog']
def fwhm_cog(data, aperture, uncertainty=None, mask=None,
aperture_method='center', subpixels=5, wcs=None):
"""
Compute the full width at half max of a source using the "curve of growth"
method.
Parameters
----------
data : array_like, `~astropy.units.Quantity`, `~astropy.nddata.NDData`
The 2D array on which to perform photometry. ``data`` should be
background-subtracted. If ``data`` is a
`~astropy.units.Quantity` array, then ``uncertainty`` (if input) must
also be a `~astropy.units.Quantity` array with the same units.
See the Notes section below for more information about
`~astropy.nddata.NDData` input.
aperture : `~photutils.aperture.Aperture` or a tupople
The aperture to use for calculating the FWHM. For convenience, this can
also be a length-3 tuple giving the (xcen, ycen, radius), or a length-2
tuple (skycoord, radius) which are equivalent to passing in
`~photutils.aperture.CircularAperture` or
`~photutils.aperture.SkyCircularAperture`. Note that if passing an
aperture in it must be a single scalar position.
uncertainty : array_like or `~astropy.units.Quantity` or bool, optional
The pixel-wise Gaussian 1-sigma uncertainty of the input ``data``.
These are used to weight pixels using inverse-variance weighting. If a
`~astropy.units.Quantity` array, then ``data`` must also be a
`~astropy.units.Quantity` array with the same units.
If False, the weighting step will be skipped even if uncertainties are
present from a passed-in `NDData` object
mask : array_like (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.
aperture_method : {'exact', 'center', 'subpixel'}, optional
The method used to determine the overlap of the aperture on the
pixel grid. Not all options are available for all aperture
types. Note that the more precise methods are generally slower.
The following methods are available:
* ``'exact'`` (default):
The the exact fractional overlap of the aperture and
each pixel is calculated. The returned mask 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 returned mask 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 returned mask
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.
wcs : WCS object, 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`). Used only if the input
``aperture`` is expressed in sky coordinates.
Returns
-------
result
"""
# this parsing is from aperture_photometry. TODO: break out into a utility function shared between these functions
if isinstance(data, NDData):
nddata_attr = {'uncertainty': uncertainty, 'mask': mask, 'wcs': wcs}
for key, value in nddata_attr.items():
if value is not None:
warnings.warn(f'The {key!r} keyword will 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:
uncertainty = data.uncertainty.array
else:
uncertainty = 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 fwhm_cog(data, aperture, uncertainty=uncertainty, mask=mask, wcs=wcs, aperture_method=aperture_method, subpixels=subpixels)
# validate inputs
data, uncertainty = _validate_inputs(data, uncertainty)
# handle data, error, and unit inputs
# output data and error are ndarray without units
data, uncertainty, unit = _handle_units(data, uncertainty)
# compute variance and apply input mask
data, variance = _prepare_photometry_data(data, uncertainty, mask)
# now the fwhm-specific bits
if not isinstance(aperture, Aperture):
try:
nap = len(aperture)
except TypeError:
raise TypeError('aperture is not an Aperture object or a sequence')
if nap == 3:
aperture = CircularAperture(aperture[:2], aperture[2])
#TODO: allow a unitful radius that uses the platescale at the given location as a best-guess
elif nap == 2:
aperture = SkyCircularAperture(*aperture)
else:
raise ValueError('aperture is not length 2 or 3 nor an Aperture')
if not aperture.isscalar:
raise ValueError('aperture has a list of positions, not a single position')
if isinstance(aperture, SkyAperture):
pixel_aperture = aperture.to_pixel(wcs)
assert pixel_aperture is not aperture, 'These should always be different now' #TODO: can remove this if someone else is certain this is
else:
pixel_aperture = aperture
apermask = aperture.to_mask(method=aperture_method, subpixels=subpixels)
masked_values = apermask.get_values(data)
# coordinates masking cannot be "exact" because that's fractional, but subpixel still makes sense
coomask = aperture.to_mask(method='center' if aperture_method=='exact' else aperture_method,
subpixels=subpixels)
ygrid, xgrid = np.meshgrid(*[np.arange(sh) for sh in data.shape], copy=False)
xs = coomask.get_values(xgrid.T)
ys = coomask.get_values(ygrid.T)
if aperture is pixel_aperture:
# a pixel-space COG
dx = xs - pixel_aperture.positions[0]
dy = ys - pixel_aperture.positions[1]
px_distance = np.hypot(dx, dy)
dsortidx = np.argsort(px_distance)
dsort = px_distance[dsortidx]
else:
# an on-sky COG
skycoords = wcs.to_world(xs, ys)
sky_distance = aperture.positions.separation(skycoords)
dsortidx = np.argsort(sky_distance)
dsort = sky_distance[dsortidx]
curve_of_growth = np.cumsum(masked_values[dsortidx])
hwhm = np.interp(curve_of_growth[-1]/2., curve_of_growth, dsort) # interp works correctly for unitful `fp`!
return hwhm*2, dsort, curve_of_growth