diff --git a/CHANGES.rst b/CHANGES.rst
index e27820312e..ed8cff47fe 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -27,6 +27,10 @@ Cubeviz
Imviz
^^^^^
+- Added the ability to fit Gaussian1D model to radial profile in
+ Simple Aperture Photometry plugin. Radial profile and curve of growth now center
+ on source centroid, not Subset center. [#1409]
+
Mosviz
^^^^^^
diff --git a/docs/imviz/plugins.rst b/docs/imviz/plugins.rst
index b683640810..0571d4a92d 100644
--- a/docs/imviz/plugins.rst
+++ b/docs/imviz/plugins.rst
@@ -160,7 +160,9 @@ an interactively selected region. A typical workflow is as follows:
Caution: having too many data points may cause performance issues with this feature.
The exact limitations depend on your hardware.
-10. Once all inputs are populated correctly, click on the :guilabel:`CALCULATE`
+10. Toggle :guilabel:`Fit Gaussian` on to fit a `~astropy.modeling.functional_models.Gaussian1D`
+ model to the radial profile data. This is disabled for curve-of-growth.
+11. Once all inputs are populated correctly, click on the :guilabel:`CALCULATE`
button to perform simple aperture photometry.
.. note::
@@ -169,8 +171,8 @@ an interactively selected region. A typical workflow is as follows:
However, if NaN exists in data, it will be treated as 0.
When calculation is complete, a plot would show the radial profile
-of the background subtracted data and the photometry results are displayed under the
-:guilabel:`CALCULATE` button.
+of the background subtracted data and the photometry and model fitting (if requested)
+results are displayed under the :guilabel:`CALCULATE` button.
.. figure:: img/imviz_radial_profile.png
:alt: Imviz radial profile plot.
@@ -182,7 +184,14 @@ of the background subtracted data and the photometry results are displayed under
Radial profile (raw).
-You can also retrieve the results as `~astropy.table.QTable` as follows,
+If you opted to fit a `~astropy.modeling.functional_models.Gaussian1D`
+to the radial profile, the last fitted model parameters will be displayed
+under the radial profile plot. The model itself can be obtained by as follows.
+See :ref:`astropy:astropy-modeling` on how to manipulate the model::
+
+ my_gaussian1d = imviz.app.fitted_models['phot_radial_profile']
+
+You can also retrieve the photometry results as `~astropy.table.QTable` as follows,
assuming ``imviz`` is the instance of your Imviz application::
results = imviz.get_aperture_photometry_results()
diff --git a/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py b/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py
index 2aca9f8c69..a49eba7569 100644
--- a/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py
+++ b/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py
@@ -1,8 +1,13 @@
+import os
+import warnings
from datetime import datetime
import bqplot
import numpy as np
from astropy import units as u
+from astropy.modeling.fitting import LevMarLSQFitter
+from astropy.modeling import Parameter
+from astropy.modeling.models import Gaussian1D
from astropy.table import QTable
from astropy.time import Time
from ipywidgets import widget_serialization
@@ -41,6 +46,8 @@ class SimpleAperturePhotometry(TemplateMixin, DatasetSelectMixin):
current_plot_type = Unicode().tag(sync=True)
plot_available = Bool(False).tag(sync=True)
radial_plot = Any('').tag(sync=True, **widget_serialization)
+ fit_radial_profile = Bool(False).tag(sync=True)
+ fit_results = List().tag(sync=True)
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
@@ -63,6 +70,7 @@ def __init__(self, *args, **kwargs):
self._fig = bqplot.Figure()
self.plot_types = ["Curve of Growth", "Radial Profile", "Radial Profile (Raw)"]
self.current_plot_type = self.plot_types[0]
+ self._fitted_model_name = 'phot_radial_profile'
def reset_results(self):
self.result_available = False
@@ -218,6 +226,11 @@ def vue_do_aper_phot(self, *args, **kwargs):
data = self._selected_data
reg = self._selected_subset
+ # Reset last fitted model
+ fit_model = None
+ if self._fitted_model_name in self.app.fitted_models:
+ del self.app.fitted_models[self._fitted_model_name]
+
try:
comp = data.get_component(data.main_components[0])
try:
@@ -319,16 +332,17 @@ def vue_do_aper_phot(self, *args, **kwargs):
line_y_sc = bqplot.LinearScale()
if self.current_plot_type == "Curve of Growth":
- self._fig.title = 'Curve of growth from Subset center'
+ self._fig.title = 'Curve of growth from source centroid'
x_arr, sum_arr, x_label, y_label = _curve_of_growth(
- comp_data, aperture, phot_table['sum'][0], wcs=data.coords,
- background=bg, pixarea_fac=pixarea_fac)
+ comp_data, phot_aperstats.centroid, aperture, phot_table['sum'][0],
+ wcs=data.coords, background=bg, pixarea_fac=pixarea_fac)
self._fig.axes = [bqplot.Axis(scale=line_x_sc, label=x_label),
bqplot.Axis(scale=line_y_sc, orientation='vertical',
label=y_label)]
bqplot_line = bqplot.Lines(x=x_arr, y=sum_arr, marker='circle',
scales={'x': line_x_sc, 'y': line_y_sc},
marker_size=32, colors='gray')
+ bqplot_marks = [bqplot_line]
else: # Radial profile
self._fig.axes = [bqplot.Axis(scale=line_x_sc, label='pix'),
@@ -336,21 +350,50 @@ def vue_do_aper_phot(self, *args, **kwargs):
label=comp.units or 'Value')]
if self.current_plot_type == "Radial Profile":
- self._fig.title = 'Radial profile from Subset center'
+ self._fig.title = 'Radial profile from source centroid'
x_data, y_data = _radial_profile(
- phot_aperstats.data_cutout, phot_aperstats.bbox, aperture, raw=False)
+ phot_aperstats.data_cutout, phot_aperstats.bbox, phot_aperstats.centroid,
+ raw=False)
bqplot_line = bqplot.Lines(x=x_data, y=y_data, marker='circle',
scales={'x': line_x_sc, 'y': line_y_sc},
marker_size=32, colors='gray')
else: # Radial Profile (Raw)
- self._fig.title = 'Raw radial profile from Subset center'
- radial_r, radial_img = _radial_profile(
- phot_aperstats.data_cutout, phot_aperstats.bbox, aperture, raw=True)
- bqplot_line = bqplot.Scatter(x=radial_r, y=radial_img, marker='circle',
+ self._fig.title = 'Raw radial profile from source centroid'
+ x_data, y_data = _radial_profile(
+ phot_aperstats.data_cutout, phot_aperstats.bbox, phot_aperstats.centroid,
+ raw=True)
+ bqplot_line = bqplot.Scatter(x=x_data, y=y_data, marker='circle',
scales={'x': line_x_sc, 'y': line_y_sc},
default_size=1, colors='gray')
- self._fig.marks = [bqplot_line]
+ # Fit Gaussian1D to radial profile data.
+ # mean is fixed at 0 because we recentered to centroid.
+ if self.fit_radial_profile:
+ fitter = LevMarLSQFitter()
+ y_max = y_data.max()
+ std = 0.5 * (phot_table['semimajor_sigma'][0] +
+ phot_table['semiminor_sigma'][0])
+ if isinstance(std, u.Quantity):
+ std = std.value
+ gs = Gaussian1D(amplitude=y_max, mean=0, stddev=std,
+ fixed={'mean': True, 'amplitude': True},
+ bounds={'amplitude': (y_max * 0.5, y_max)})
+ with warnings.catch_warnings(record=True) as warns:
+ fit_model = fitter(gs, x_data, y_data)
+ if len(warns) > 0:
+ msg = os.linesep.join([str(w.message) for w in warns])
+ self.hub.broadcast(SnackbarMessage(
+ f"Radial profile fitting: {msg}", color='warning', sender=self))
+ y_fit = fit_model(x_data)
+ self.app.fitted_models[self._fitted_model_name] = fit_model
+ bqplot_fit = bqplot.Lines(x=x_data, y=y_fit, marker=None,
+ scales={'x': line_x_sc, 'y': line_y_sc},
+ colors='magenta', line_style='dashed')
+ bqplot_marks = [bqplot_line, bqplot_fit]
+ else:
+ bqplot_marks = [bqplot_line]
+
+ self._fig.marks = bqplot_marks
except Exception as e: # pragma: no cover
self.reset_results()
@@ -379,7 +422,18 @@ def vue_do_aper_phot(self, *args, **kwargs):
f'{x:.4e} ({phot_table["aperture_sum_counts_err"][0]:.4e})'})
else:
tmp.append({'function': key, 'result': str(x)})
+
+ # Also display fit results
+ fit_tmp = []
+ if fit_model is not None and isinstance(fit_model, Gaussian1D):
+ for param in ('fwhm', 'amplitude'): # mean is fixed at 0
+ p_val = getattr(fit_model, param)
+ if isinstance(p_val, Parameter):
+ p_val = p_val.value
+ fit_tmp.append({'function': param, 'result': f'{p_val:.4e}'})
+
self.results = tmp
+ self.fit_results = fit_tmp
self.result_available = True
self.radial_plot = self._fig
self.bqplot_figs_resize = [self._fig]
@@ -389,7 +443,7 @@ def vue_do_aper_phot(self, *args, **kwargs):
# NOTE: These are hidden because the APIs are for internal use only
# but we need them as a separate functions for unit testing.
-def _radial_profile(radial_cutout, reg_bb, aperture, raw=False):
+def _radial_profile(radial_cutout, reg_bb, centroid, raw=False):
"""Calculate radial profile.
Parameters
@@ -400,8 +454,8 @@ def _radial_profile(radial_cutout, reg_bb, aperture, raw=False):
reg_bb : obj
Bounding box from ``ApertureStats``.
- aperture : obj
- ``photutils`` aperture object.
+ centroid : tuple of int
+ ``ApertureStats`` centroid or desired center in ``(x, y)``.
raw : bool
If `True`, returns raw data points for scatter plot.
@@ -409,14 +463,15 @@ def _radial_profile(radial_cutout, reg_bb, aperture, raw=False):
"""
reg_ogrid = np.ogrid[reg_bb.iymin:reg_bb.iymax, reg_bb.ixmin:reg_bb.ixmax]
- radial_dx = reg_ogrid[1] - aperture.positions[0]
- radial_dy = reg_ogrid[0] - aperture.positions[1]
+ radial_dx = reg_ogrid[1] - centroid[0]
+ radial_dy = reg_ogrid[0] - centroid[1]
radial_r = np.hypot(radial_dx, radial_dy)[~radial_cutout.mask].ravel() # pix
radial_img = radial_cutout.compressed() # data unit
if raw:
- x_arr = radial_r
- y_arr = radial_img
+ i_arr = np.argsort(radial_r)
+ x_arr = radial_r[i_arr]
+ y_arr = radial_img[i_arr]
else:
# This algorithm is from the imexam package,
# see licenses/IMEXAM_LICENSE.txt for more details
@@ -427,7 +482,7 @@ def _radial_profile(radial_cutout, reg_bb, aperture, raw=False):
return x_arr, y_arr
-def _curve_of_growth(data, aperture, final_sum, wcs=None, background=0, n_datapoints=10,
+def _curve_of_growth(data, centroid, aperture, final_sum, wcs=None, background=0, n_datapoints=10,
pixarea_fac=None):
"""Calculate curve of growth for aperture photometry.
@@ -436,8 +491,14 @@ def _curve_of_growth(data, aperture, final_sum, wcs=None, background=0, n_datapo
data : ndarray or `~astropy.units.Quantity`
Data for the calculation.
+ centroid : tuple of int
+ ``ApertureStats`` centroid or desired center in ``(x, y)``.
+
aperture : obj
- ``photutils`` aperture object.
+ ``photutils`` aperture to use, except its center will be
+ changed to the given ``centroid``. This is because the aperture
+ might be hand-drawn and a more accurate centroid has been
+ recalculated separately.
final_sum : float or `~astropy.units.Quantity`
Aperture sum that is already calculated in the
@@ -477,20 +538,20 @@ def _curve_of_growth(data, aperture, final_sum, wcs=None, background=0, n_datapo
if isinstance(aperture, CircularAperture):
x_label = 'Radius (pix)'
x_arr = np.linspace(0, aperture.r, num=n_datapoints)[1:]
- aper_list = [CircularAperture(aperture.positions, cur_r) for cur_r in x_arr[:-1]]
+ aper_list = [CircularAperture(centroid, cur_r) for cur_r in x_arr[:-1]]
elif isinstance(aperture, EllipticalAperture):
x_label = 'Semimajor axis (pix)'
x_arr = np.linspace(0, aperture.a, num=n_datapoints)[1:]
a_arr = x_arr[:-1]
b_arr = aperture.b * a_arr / aperture.a
- aper_list = [EllipticalAperture(aperture.positions, cur_a, cur_b, theta=aperture.theta)
+ aper_list = [EllipticalAperture(centroid, cur_a, cur_b, theta=aperture.theta)
for (cur_a, cur_b) in zip(a_arr, b_arr)]
elif isinstance(aperture, RectangularAperture):
x_label = 'Width (pix)'
x_arr = np.linspace(0, aperture.w, num=n_datapoints)[1:]
w_arr = x_arr[:-1]
h_arr = aperture.h * w_arr / aperture.w
- aper_list = [RectangularAperture(aperture.positions, cur_w, cur_h, theta=aperture.theta)
+ aper_list = [RectangularAperture(centroid, cur_w, cur_h, theta=aperture.theta)
for (cur_w, cur_h) in zip(w_arr, h_arr)]
else:
raise TypeError(f'Unsupported aperture: {aperture}')
diff --git a/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue b/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue
index f913d891b8..a4d285f24b 100644
--- a/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue
+++ b/jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue
@@ -110,6 +110,16 @@
>
+
+
+
+
+
+
Calculate
@@ -123,8 +133,25 @@
+
+ Gaussian Fit Results
+
+ Result
+ Value
+
+
+
+ {{ item.function }}
+
+ {{ item.result }}
+
+
+
- Results
+ Photometry Results
Result
Value
diff --git a/jdaviz/configs/imviz/tests/test_simple_aper_phot.py b/jdaviz/configs/imviz/tests/test_simple_aper_phot.py
index 668efbcdb4..ac602550b4 100644
--- a/jdaviz/configs/imviz/tests/test_simple_aper_phot.py
+++ b/jdaviz/configs/imviz/tests/test_simple_aper_phot.py
@@ -18,6 +18,10 @@ def test_plugin_wcs_dithered(self):
phot_plugin = self.imviz.app.get_tray_item_from_name('imviz-aper-phot-simple')
+ # Model fitting is already tested in astropy.
+ # Here, we enable it just to make sure it does not crash.
+ phot_plugin.fit_radial_profile = True
+
# Make sure invalid Data/Subset selection does not crash plugin.
with pytest.raises(ValueError):
phot_plugin.dataset_selected = 'no_such_data'
@@ -164,7 +168,7 @@ def test_plugin_wcs_dithered(self):
# Curve of growth
phot_plugin.current_plot_type = 'Curve of Growth'
phot_plugin.vue_do_aper_phot()
- assert phot_plugin._fig.title == 'Curve of growth from Subset center'
+ assert phot_plugin._fig.title == 'Curve of growth from source centroid'
class TestSimpleAperPhot_NoWCS(BaseImviz_WCS_NoWCS):
@@ -255,13 +259,14 @@ def test_annulus_background(imviz_helper):
class TestRadialProfile():
def setup_class(self):
data = np.ones((51, 51)) * u.nJy
- self.aperture = EllipticalAperture((25, 25), 20, 15)
- phot_aperstats = ApertureStats(data, self.aperture)
+ aperture = EllipticalAperture((25, 25), 20, 15)
+ phot_aperstats = ApertureStats(data, aperture)
self.data_cutout = phot_aperstats.data_cutout
self.bbox = phot_aperstats.bbox
+ self.centroid = phot_aperstats.centroid
def test_profile_raw(self):
- x_arr, y_arr = _radial_profile(self.data_cutout, self.bbox, self.aperture, raw=True)
+ x_arr, y_arr = _radial_profile(self.data_cutout, self.bbox, self.centroid, raw=True)
# Too many data points to compare each one for X.
assert x_arr.shape == y_arr.shape == (923, )
assert_allclose(x_arr.min(), 0)
@@ -269,7 +274,7 @@ def test_profile_raw(self):
assert_allclose(y_arr, 1)
def test_profile_imexam(self):
- x_arr, y_arr = _radial_profile(self.data_cutout, self.bbox, self.aperture, raw=False)
+ x_arr, y_arr = _radial_profile(self.data_cutout, self.bbox, self.centroid, raw=False)
assert_allclose(x_arr, range(20))
assert_allclose(y_arr, 1)
@@ -293,11 +298,12 @@ def test_curve_of_growth(with_unit):
RectangularAperture(cen, 20, 15))
for aperture in apertures:
- final_sum = ApertureStats(data, aperture).sum
+ astat = ApertureStats(data, aperture)
+ final_sum = astat.sum
if pixarea_fac is not None:
final_sum = final_sum * pixarea_fac
x_arr, sum_arr, x_label, y_label = _curve_of_growth(
- data, aperture, final_sum, background=bg, pixarea_fac=pixarea_fac)
+ data, astat.centroid, aperture, final_sum, background=bg, pixarea_fac=pixarea_fac)
assert_allclose(x_arr, [2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
assert y_label == expected_ylabel
@@ -316,5 +322,5 @@ def test_curve_of_growth(with_unit):
assert_allclose(sum_arr, [3, 12, 27, 48, 75, 108, 147, 192, 243, 300])
with pytest.raises(TypeError, match='Unsupported aperture'):
- _curve_of_growth(data, EllipticalAnnulus(cen, 3, 8, 5), 100,
+ _curve_of_growth(data, cen, EllipticalAnnulus(cen, 3, 8, 5), 100,
pixarea_fac=pixarea_fac)
diff --git a/notebooks/concepts/imviz_simple_aper_phot.ipynb b/notebooks/concepts/imviz_simple_aper_phot.ipynb
index e0549ed983..55de8717be 100644
--- a/notebooks/concepts/imviz_simple_aper_phot.ipynb
+++ b/notebooks/concepts/imviz_simple_aper_phot.ipynb
@@ -190,10 +190,29 @@
"results"
]
},
+ {
+ "cell_type": "markdown",
+ "id": "b45fcdfe",
+ "metadata": {},
+ "source": [
+ "If you fitted Gaussian to radial profile, you can get it back out like this. If it does not exist, you will get `None`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bcb594f2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "my_gaussian = imviz.app.fitted_models.get('phot_radial_profile', None)\n",
+ "my_gaussian"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
- "id": "00435609",
+ "id": "b7c59fef",
"metadata": {},
"outputs": [],
"source": []