Add masking and limits to pstat

wfc3tools/pstat.help

@@ -68,6 +68,7 @@ Parameters:
 
 * extname = "sci" [string, allowed values: sci | err | dq ]
     Extension name (EXTNAME keyword value) of data to plot.
+
 * units = "counts" [string, allowed values: counts | rate]
     Plot "sci" or  "err"  data  in  units  of  counts  or  countrate
     ("rate").   Input data can be in either unit; conversion will be
@@ -77,6 +78,20 @@ Parameters:
 * stat = "midpt" [string, allowed values: mean|midpt|mode|stddev|min|max]
     Type of statistic to compute.
 
+* mask = None [`numpy.ndarray` (bool), optional]
+    A boolean mask with the same shape as ``data``, where a `True`
+    value indicates the corresponding element of ``data`` is masked.
+    Masked pixels are excluded when computing the statistics.
+
+* low_limit = None [float, (optional)]
+       If set, the statistics will be calculated using only pixels that are
+       above this value.
+
+*  high_limit = None [float, (optional)]
+       If set, the statistics will be calculated using only pixels that are
+       below this value.
+
+
 * title = "" [string]
    Title  for  the  plot.   If  left  blank,  the name of the input
    image, appended with the extname and image section, is used.
@@ -92,6 +107,10 @@ Parameters:
 
 * plot = True [bool]  set plot to false if you only want the data returned
 
+* overplot = False [bool] 
+   If True, the results will be overplotted on the previous plot
+
+
 Usage:
 
 ::

wfc3tools/pstat.py

@@ -17,11 +17,77 @@
 
 plt.ion()
 
-def pstat(filename, extname="sci", units="counts", stat="midpt", title=None,
-          xlabel=None, ylabel=None, plot=True):
+
+def pstat(filename, extname="sci", units="counts", stat="midpt", mask=None,
+          low_limit = None, high_limit = None, title=None, xlabel=None,
+          ylabel=None, plot=True, overplot=False):
     """A function to plot the statistics of one or more pixels up an IR ramp.
 
-    Pixel values here are 0 based, not 1 based """
+    Parameters
+    ----------
+    filename: string
+       Input   MultiAccum   image  name  with  optional  image  section
+       specification.  If no image section  is  specified,  the  entire image
+       is  used.   This  should  be  either a _raw or _ima file, containing all
+       the  data  from  multiple  readouts.   You  must specify  just  the
+       file name and image section, with no extname designation.
+
+    extname:  {"sci", "err", "dq"}
+       Extension name (EXTNAME keyword value) of data to plot.
+
+    units: {"counts", "rate"}
+       Plot "sci" or  "err"  data  in  units  of  counts  or  countrate
+       ("rate").   Input data can be in either unit; conversion will be
+       performed automatically.  Ignored when  plotting  "dq",  "samp", or
+       "time" data.
+
+    stat: { "mean", "midpt", "mode", "stddev", "min", "max"}
+       Type of statistic to compute.
+
+    mask: `numpy.ndarray` (bool), optional
+       A boolean mask with the same shape as ``data``, where a `True`
+       value indicates the corresponding element of ``data`` is masked.
+       Masked pixels are excluded when computing the statistics.
+
+    low_limit: float, (optional)
+       If set, the statistics will be calculated using only pixels that are
+       above this value.
+
+    high_limit: float, (optional)
+       If set, the statistics will be calculated using only pixels that are
+       below this value.
+
+    title: str
+       Title  for  the  plot.   If  left  blank,  the name of the input image,
+       appended with the extname and image section, is used.
+
+    xlabel: str
+       Label for the X-axis of the plot.  If  left  blank,  a  suitable default
+       is generated.
+
+    ylabel: str
+       Label  for  the  Y-axis  of  the plot. If left blank, a suitable default
+       based on the plot units and the extname of the  data  is generated.
+
+    plot: Bool
+       Set plot to false if you only want the data returned
+
+    overplot: Bool
+       If True, the results will be overplotted on the previous plot
+
+    Returns
+    -------
+    xaxis: `numpy.ndarray`
+       Array of x-axis values that will be plotted
+
+    yaxis: `numpuy.ndarray`
+       Array of y-axis values that will be plotted as specified by 'units'
+
+
+    Notes
+    -----
+    Pixel values here are 0 based, not 1 based
+    """
 
     # pull the image extension from the filename string
     section_start = filename.find("[")
@@ -98,25 +164,41 @@ def pstat(filename, extname="sci", units="counts", stat="midpt", title=None,
             yend = myfile[1].header["NAXIS2"]  # full y size
 
         for i in range(1, nsamp, 1):
+            data = myfile[extname.upper(), i].data[xstart:xend, ystart:yend]
+
+            # mask the data and remove outlyier values
+            if mask is not None:
+                data = data[~mask[xstart:xend, ystart:yend]]
+
+            if high_limit is not None:
+                data = data[data < high_limit]
+
+            if low_limit is not None:
+                data = data[data > low_limit]
+
+            if data.size == 0:
+                print("No valid pixels in ext {} of {}".format(i, imagename))
+                return xaxis, yaxis
+
             if "midpt" in stat:
-                yaxis[i-1] = np.median(myfile[extname.upper(), i].data[xstart:xend, ystart:yend])
+                yaxis[i-1] = np.median(data)
 
             if "mean" in stat:
-                yaxis[i-1] = np.mean(myfile[extname.upper(), i].data[xstart:xend, ystart:yend])
+                yaxis[i-1] = np.mean(data)
 
             if "mode" in stat:
-                yaxis[i-1] = mode(myfile[extname.upper(), i].data[xstart:xend, ystart:yend], axis=None)[0]
+                yaxis[i-1] = mode(data, axis=None)[0]
 
             if "min" in stat:
-                yaxis[i-1] = np.min(myfile[extname.upper(), i].data[xstart:xend, ystart:yend])
+                yaxis[i-1] = np.min(data)
 
             if "max" in stat:
-                yaxis[i-1] = np.max(myfile[extname.upper(), i].data[xstart:xend, ystart:yend])
+                yaxis[i-1] = np.max(data)
 
             if "stddev" in stat:
-                yaxis[i-1] = np.std(myfile[extname.upper(), i].data[xstart:xend, ystart:yend])
+                yaxis[i-1] = np.std(data)
 
-            exptime=myfile["SCI",i ].header['SAMPTIME']
+            exptime = myfile["SCI", i].header['SAMPTIME']
             xaxis[i-1] = exptime
 
             # convert to countrate
@@ -127,13 +209,14 @@ def pstat(filename, extname="sci", units="counts", stat="midpt", title=None,
                 yaxis[i-1] *= exptime
 
     if plot:
-        plt.clf()  # clear out any current plot
+        if not overplot:
+            plt.clf()  # clear out any current plot
         if not ylabel:
             if "rate" in units.lower():
                 if "/" in bunit.lower():
                     ylabel = bunit
                 else:
-                    ylabel = bunit+" per second"
+                    ylabel = bunit + " per second"
             else:
                 if "/" in bunit:
                     stop_index = bunit.find("/")
@@ -190,4 +273,5 @@ def help(file=None):
         f.write(helpstr)
         f.close()
 
+
 pstat.__doc__ = getHelpAsString(docstring=True)