Gaussian1D for radial profile

Author: pllim

CHANGES.rst

@@ -23,6 +23,9 @@ Cubeviz
 Imviz
 ^^^^^
 
+- Added the ability to fit Gaussian1D model to radial profile in
+  Simple Aperture Photometry plugin. [#1409]
+
 Mosviz
 ^^^^^^
 

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::
@@ -182,7 +184,7 @@ 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,
+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()
@@ -258,6 +260,12 @@ The columns are as follow:
 Once you have the results in a table, you can further manipulated them as
 documented in :ref:`astropy:astropy-table`.
 
+If you opted to fit a `~astropy.modeling.functional_models.Gaussian1D`
+to the radial profile, the last fitted model 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']
+
 .. _imviz-export-plot:
 
 Export Plot

jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.py

@@ -1,8 +1,12 @@
+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.models import Gaussian1D
 from astropy.table import QTable
 from astropy.time import Time
 from ipywidgets import widget_serialization
@@ -41,6 +45,7 @@ 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)
 
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
@@ -329,6 +334,7 @@ def vue_do_aper_phot(self, *args, **kwargs):
                 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'),
@@ -344,13 +350,37 @@ def vue_do_aper_phot(self, *args, **kwargs):
                                                marker_size=32, colors='gray')
                 else:  # Radial Profile (Raw)
                     self._fig.title = 'Raw radial profile from Subset center'
-                    radial_r, radial_img = _radial_profile(
+                    x_data, y_data = _radial_profile(
                         phot_aperstats.data_cutout, phot_aperstats.bbox, aperture, raw=True)
-                    bqplot_line = bqplot.Scatter(x=radial_r, y=radial_img, marker='circle',
+                    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.
+                if self.fit_radial_profile:
+                    fitter = LevMarLSQFitter()
+                    y_max = y_data.max()
+                    x_mean = x_data[np.where(y_data == y_max)].mean()
+                    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=x_mean, stddev=std)
+                    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['phot_radial_profile'] = 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()

jdaviz/configs/imviz/plugins/aper_phot_simple/aper_phot_simple.vue

@@ -112,6 +112,15 @@
           ></v-select>
         </v-row>
 
+        <v-row v-if="current_plot_type != 'Curve of Growth'">
+          <v-switch
+            label="Fit Gaussian"
+            hint="Fit Gaussian1D to radial profile"
+            v-model="fit_radial_profile"
+            persistent-hint>
+          </v-switch>
+        </v-row>
+
         <v-row justify="end">
           <v-btn color="primary" text @click="do_aper_phot">Calculate</v-btn>
         </v-row>

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'

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."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bcb594f2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "my_gaussian = imviz.app.fitted_models['phot_radial_profile']\n",
+    "my_gaussian"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "00435609",
+   "id": "1e546682",
    "metadata": {},
    "outputs": [],
    "source": []