fixes for #165, #238, #210, #205

CHANGES.rst

@@ -1,6 +1,8 @@
 version 0.9.2.dev (unreleased)
 ------------------------------
-
+ - updated the gauss1d to set fixed fitting parameters and improve FWHM estimate(#165, #238)
+ - fixed background sky subtraction (#210)
+ - formatting fix for print statement in aper phot plot (#205)
 
 version 0.9.1 (2020-03-30)
 --------------------------

imexam/connect.py

@@ -233,7 +233,7 @@ def _run_imexam(self):
                 # This loop now recognizes the arrow keys
                 # for moving the cursor in the window, it calls
                 # the windows cursor function. However, depending
-                # on how the use has their windowing focus setup
+                # on how the user has their windowing focus setup
                 # they might loose focus on the DS9 window and have to
                 # move the cursor to gain focus again, is there a way
                 # around this? Cursor is not implemented in the ginga

imexam/imexam_defpars.py

@@ -43,6 +43,7 @@
                        "fitplot": [True, "Overplot model fit?"],
                        "func": ["Gaussian1D", "Model form to fit"],
                        "center": [True, "Solve for center using 2d Gaussian? [bool]"],
+                       "delta": [10, "bounding box for centering measurement"],
                        "background": [True, "Subtract background? [bool]"],
                        "skyrad": [10., "Background inner radius in pixels, from center of object"],
                        "width": [5., "Background annulus width in pixels"],

imexam/imexamine.py

@@ -91,16 +91,17 @@ def __init__(self):
         self._define_default_pars()
         # default plot name saved with "s" key
         self.plot_name = "imexam_plot.pdf"
-        # let users have multiple plot windows, the list stores their names
+        # let users have multiple plot windows, the list stores their instance reference
         self._plot_windows = list()
         # this contains the name of the current plotting window
         self._figure_name = "imexam"
-        self._plot_windows.append(self._figure_name)
+        self._plot_windows.append(plt.gcf())
         self._reserved_keys = ['q', '2']  # not to be changed with user funcs
         self._fit_models = ["Gaussian1D",
                             "Moffat1D",
                             "MexicanHat1D",
                             "AiryDisk2D",
+                            "Gaussian2D",
                             "Polynomial1D"]
 
         # see if the package logger was already started
@@ -129,7 +130,7 @@ def setlog(self, filename=None, on=True, level=logging.INFO):
     def _close_plots(self):
         """Make sure to release plot memory at end of exam loop."""
         for plot in self._plot_windows:
-            plt.close()
+            plt.close(plot)
 
     def close(self):
         """For use with the Imexamine object standalone."""
@@ -295,8 +296,9 @@ def new_plot_window(self, x, y, data=None):
         """
         if data is None:
             data = self._data
+        # save the old figure
+        self._plot_windows.append(plt.gcf())
         self._figure_name = "imexam" + str(len(self._plot_windows) + 1)
-        self._plot_windows.append(self._figure_name)
         self.log.info(f"Plots now directed towards {self._figure_name}")
 
     def plot_line(self, x, y, data=None, fig=None):
@@ -352,6 +354,7 @@ def plot_line(self, x, y, data=None, fig=None):
             plt.pause(0.001)
         else:
             fig.canvas.draw_idle()
+        self._plot_windows[-1] = plt.gcf()
 
     def plot_column(self, x, y, data=None, fig=None):
         """column plot of data at point y.
@@ -408,6 +411,7 @@ def plot_column(self, x, y, data=None, fig=None):
             plt.pause(0.001)
         else:
             fig.canvas.draw_idle()
+        self._plot_windows[-1] = plt.gcf()
 
     def show_xy_coords(self, x, y, data=None):
         """print the x,y,value to the screen.
@@ -478,17 +482,10 @@ def save_figure(self, fig=None):
 
         Parameters
         ----------
-        data: numpy array
-            The data array to work on
         fig: figure for redirect
             Used for interaction with the ginga GUI
         """
-        if fig is None:
-            fig = plt.figure(self._figure_name)
-        ax = fig.gca()
-        fig.savefig(self.plot_name)
-        pstr = f"plot saved to {self.plot_name}"
-        self.log.info(pstr)
+        self.save(fig=fig)
 
     def save(self, filename=None, fig=None):
         """Save to file the figure that's currently displayed.
@@ -509,10 +506,10 @@ def save(self, filename=None, fig=None):
             self.set_plot_name(self._figure_name + ".pdf")
 
         if fig is None:
-            fig = plt.figure(self._figure_name)
-        ax = fig.gca()
+            fig = self._plot_windows[-1]
+
         fig.savefig(self.plot_name)
-        pstr = f"plot saved to {self.plot_name}"
+        pstr = f"plot {fig.number} saved to {self.plot_name}"
         self.log.info(pstr)
 
     def aper_phot(self, x, y, data=None,
@@ -633,7 +630,7 @@ def aper_phot(self, x, y, data=None,
 
             # Construct the output strings (header and parameter values)
             pheader = f"x\ty\tradius\tflux\tmag(zpt={magzero:0.2})\t"
-            pstr = f"\n{x:.2f}\t{y:0.2f}\t{radius}\t{total_flux:0.2}\t{mag:0.2}\t"
+            pstr = f"\n{xx:.2f}\t{yy:0.2f}\t{radius}\t{total_flux:0.2}\t{mag:0.2}\t"
 
             if subsky:
                 pheader += "sky/pix\t"
@@ -680,7 +677,6 @@ def aper_phot(self, x, y, data=None,
                                    self.aper_phot_pars['color_max'][0]]
 
                 pad = outer * 1.2  # XXX TODO: Bad magic number
-                print(xx, yy, pad)
                 ax.imshow(data[int(yy - pad):int(yy + pad),
                                int(xx - pad):int(xx + pad)],
                           vmin=color_range[0], vmax=color_range[1],
@@ -697,6 +693,7 @@ def aper_phot(self, x, y, data=None,
                     plt.pause(0.001)
                 else:
                     fig.canvas.draw_idle()
+                self._plot_windows[-1] = plt.gcf()
             else:
                 return (apertures, annulus_apertures, rawflux_table, sky_per_pix)
 
@@ -726,7 +723,7 @@ def line_fit(self, x, y, data=None, form=None, genplot=True, fig=None, col=False
         The background is currently ignored
 
         If centering is True in the parameter set, then the center
-        is fit with a 2d gaussian, not performed for Polynomial1D
+        is fit with a 2d gaussian, center is not performed for Polynomial1D
         """
 
         # Set which parameters to use
@@ -880,6 +877,7 @@ def line_fit(self, x, y, data=None, form=None, genplot=True, fig=None, col=False
                 plt.pause(0.001)
             else:
                 fig.canvas.draw_idle()
+            self._plot_windows[-1] = plt.gcf()
 
         else:
             return fitted
@@ -949,8 +947,6 @@ def com_center(self, x, y, data=None, delta=None, oversampling=1.):
         if delta >= len(data) / 4:
             delta = delta // 2
 
-        if delta is None:
-            delta = int(self.com_center_pars['delta'][0])
         xx = int(x)
         yy = int(y)
 
@@ -1047,17 +1043,30 @@ def radial_profile(self, x, y, data=None, form=None,
         """
         pars = self.radial_profile_pars
 
+        if data is None:
+            data = self._data
+
+        fitplot = bool(pars["fitplot"][0])
+        if fitplot:
+            if form is None:
+                form = getattr(models, pars["func"][0])
+            if form.name not in self._fit_models:
+                raise ValueError(f"Functional {form.name} not in available: {self._fit_models}")
+
+            self.log.info(f"using model: {form}")
+
         subtract_background = bool(pars["background"][0])
         if not photutils_installed and subtract_background:
             self.log.warning("Install photutils to enable "
                              "background subtraction")
             subtract_background = False
 
-        if data is None:
-            data = self._data
-
         getdata = bool(pars["getdata"][0])
-        center = pars["center"][0]
+        delta = int(pars["delta"][0])
+
+        # reset delta for small arrays
+        if delta >= len(data) / 4:
+            delta = delta // 2
 
         # be careful with the clipping since most
         # of the data will be near the low value
@@ -1067,41 +1076,30 @@ def radial_profile(self, x, y, data=None, form=None,
         else:
             sig_factor = 0
 
-        fitplot = bool(pars["fitplot"][0])
-
-        if fitplot:
-            if form is None:
-                fitform = getattr(models, pars["func"][0])
-            else:
-                if form not in self._fit_models:
-                    msg = f"{form} not supported for fitting"
-                    self.log.info(msg)
-                    raise ValueError(msg)
-                else:
-                    fitform = getattr(models, form)
-
         # cut the data down to size and center cutout
         datasize = int(pars["rplot"][0])
         if datasize < 3:
             self.log.info("Insufficient pixels, resetting chunk size to 3.")
             datasize = 3
 
-        if center:
-            # reset delta for small arrays
-            # make it odd if it's even
-            if ((datasize % 2) == 0):
-                datasize = datasize + 1
+        if pars["center"][0]:
             xx = int(x)
             yy = int(y)
-            #  flipped from xpa
             data_chunk = data[yy - datasize:yy + datasize,
                               xx - datasize:xx + datasize]
-            amp, centerx, centery, sigmax, sigmay = self.gauss_center(xx, yy, data, delta=datasize)
-
+            amp, xx, yy, sigma, sigmay = self.gauss_center(xx,
+                                                           yy,
+                                                           data,
+                                                           delta=delta)
+            if (xx < 0 or yy < 0 or xx > data.shape[1] or
+                    yy > data.shape[0]):
+                self.log.warning("Problem with centering, pixel coords")
+            else:
+                centerx = xx
+                centery = yy
         else:
-            centery = y
             centerx = x
-
+            centery = y
         icenterx = int(centerx)
         icentery = int(centery)
 
@@ -1157,7 +1155,7 @@ def radial_profile(self, x, y, data=None, form=None,
                 sky_per_pix = 0.
                 self.log.info("Sky background negative, setting to zero")
             self.log.info(f"Background per pixel: {sky_per_pix}")
-            flux -= sky_per_pix
+            flux = flux - sky_per_pix
 
             if getdata:
                 self.log.info(f"Sky per pixel: {sky_per_pix} using "
@@ -1168,16 +1166,16 @@ def radial_profile(self, x, y, data=None, form=None,
             self.log.info(radius, flux)
 
         # Fit the functional form to the radial profile flux
-        # TODO: Ignore sky subtracted pixels that push flux
-        # below zero?
         if fitplot:
             fline = np.linspace(0, datasize, 100)  # finer sample
             # fit model to data
-            if fitform.name == "Gaussian1D":
+            if form.name == "Gaussian1D":  # make center radii at max ********
+
                 fitted = math_helper.fit_gauss_1d(radius, flux,
                                                   sigma_factor=sig_factor,
-                                                  center_at=0,
+                                                  center_at=radius[0],
                                                   weighted=True)
+                print(f"fitline: {fitted}")
 
                 fwhmx, fwhmy = math_helper.gfwhm(fitted.stddev_0.value)
                 legend = (f"Max. pix. flux = {np.max(flux):9.3f}\n"
@@ -1187,10 +1185,10 @@ def radial_profile(self, x, y, data=None, form=None,
                 legendx = datasize / 2
                 legendy = np.max(flux) / 2
 
-            elif fitform.name == "Moffat1D":
+            elif form.name == "Moffat1D":
                 fitted = math_helper.fit_moffat_1d(flux,
                                                    sigma_factor=sig_factor,
-                                                   center_at=0,
+                                                   center_at=radius[0],
                                                    weighted=True)
                 mfwhm = math_helper.mfwhm(fitted.alpha_0.value,
                                           fitted.gamma_0.value)
@@ -1200,17 +1198,17 @@ def radial_profile(self, x, y, data=None, form=None,
                 legendx = datasize / 2
                 legendy = np.max(flux) / 2
 
-            elif fitform.name == "MexicanHat1D":
+            elif form.name == "MexicanHat1D":
                 fitted = math_helper.fit_mex_hat_1d(flux,
                                                     sigma_factor=sig_factor,
-                                                    center_at=0,
+                                                    center_at=radius[0],
                                                     weighted=True)
                 legend = (f"Max. pix. flux = {np.max(flux):9.3f}\n")
                 legendx = datasize / 2
                 legendy = np.max(flux) / 2
 
             if fitted is None:
-                msg = f"Problem with the {fitform.name} fit"
+                msg = f"Problem with the {form.name} fit"
                 self.log.info(msg)
                 raise ValueError(msg)
 
@@ -1242,14 +1240,14 @@ def radial_profile(self, x, y, data=None, form=None,
 
             if pars["title"][0] is None:
                 if fitplot:
-                    title = f"Radial Profile at ({icenterx},{icentery}) with {fitform.name}"
+                    title = f"Radial Profile at ({icenterx},{icentery}) with {form.name}"
                 else:
                     title = f"Radial Profile for {icenterx} {icentery}"
             else:
                 title = pars["title"][0]
 
             if fitplot:
-                ax.plot(fline, yfit, linestyle='-', c='r', label=fitform.name)
+                ax.plot(fline, yfit, linestyle='-', c='r', label=form.name)
                 ax.set_xlim(0, datasize, 0.5)
                 ax.text(legendx, legendy, legend)
 
@@ -1260,6 +1258,7 @@ def radial_profile(self, x, y, data=None, form=None,
                 plt.pause(0.001)
             else:
                 fig.canvas.draw_idle()
+            self._plot_windows[-1] = plt.gcf()
         else:
             return radius, flux
 
@@ -1356,6 +1355,7 @@ def curve_of_growth(self, x, y, data=None, genplot=True, fig=None):
                     plt.pause(0.001)
                 else:
                     fig.canvas.draw_idle()
+                self._plot_windows[-1] = plt.gcf()
 
             else:
                 return rapert, flux
@@ -1525,6 +1525,7 @@ def histogram(self, x, y, data=None, genplot=True, fig=None):
                 plt.pause(0.001)
             else:
                 fig.canvas.draw_idle()
+            self._plot_windows[-1] = plt.gcf()
         else:
             hist, bin_edges = np.histogram(flat_data,
                                            num_bins,
@@ -1598,6 +1599,7 @@ def contour(self, x, y, data=None, fig=None):
             plt.pause(0.001)
         else:
             fig.canvas.draw_idle()
+        self._plot_windows[-1] = plt.gcf()
 
     def surface(self, x, y, data=None, fig=None):
         """plot a surface around the specified location.
@@ -1622,7 +1624,7 @@ def surface(self, x, y, data=None, fig=None):
         if fig is None:
             fig = plt.figure(self._figure_name)
         fig.clf()
-        fig.add_subplot(111)
+        fig.add_subplot(111, projection='3d')
         ax = fig.gca(projection='3d')
 
         title = self.surface_pars["title"][0]
@@ -1704,6 +1706,7 @@ def surface(self, x, y, data=None, fig=None):
             plt.pause(0.001)
         else:
             fig.canvas.draw_idle()
+        self._plot_windows[-1] = plt.gcf()
 
     def cutout(self, x, y, data=None, size=None, fig=None):
         """Make a fits cutout around the pointer location without wcs.