Merge pull request #433 from kevin218/flexo

Fig 3304 + BG bug fix + v0.4

docs/source/installation.rst

@@ -6,8 +6,8 @@ Installation methods
 --------------------
 
 In order to have consistent, repeatable results across the ``Eureka!`` user community, we recommend that all general users install
-the most recent stable release of ``Eureka!``, v0.3. The following installation instructions are written with this in mind,
-and the most recent stable release is also available as a zipped archive `here <https://github.com/kevin218/Eureka/releases/tag/v0.3>`_.
+the most recent stable release of ``Eureka!``, v0.4. The following installation instructions are written with this in mind,
+and the most recent stable release is also available as a zipped archive `here <https://github.com/kevin218/Eureka/releases/tag/v0.4>`_.
 Also note that if you are using a macOS device with an M1 processor, you will need to use the ``conda`` environment.yml file
 installation instructions below as the pip dependencies fail to build on the M1 processor.
 
@@ -31,7 +31,7 @@ Once in your new conda environment, you can install ``Eureka!`` directly from so
 
 .. code-block:: bash
 
-	git clone -b v0.3 https://github.com/kevin218/Eureka.git
+	git clone -b v0.4 https://github.com/kevin218/Eureka.git
 	cd Eureka
 	pip install -e '.[jwst]'
 
@@ -49,7 +49,7 @@ Once in your new conda environment, you can install the ``Eureka!`` package with
 
 .. code-block:: bash
 
-	pip install -e git+https://github.com/kevin218/Eureka.git@v0.3#egg=eureka[jwst]
+	pip install -e git+https://github.com/kevin218/Eureka.git@v0.4#egg=eureka[jwst]
 
 Other specific branches can be installed using:
 
@@ -68,7 +68,7 @@ you'd prefer not to use ``pip`` to install dependencies. This can be done follow
 
 .. code-block:: bash
 
-	git clone -b v0.3 https://github.com/kevin218/Eureka.git
+	git clone -b v0.4 https://github.com/kevin218/Eureka.git
 	cd Eureka
 	conda env create --file environment.yml --force
 	conda activate eureka

src/eureka/S3_data_reduction/background.py

@@ -116,7 +116,7 @@ def writeBG_WFC3(arg):
                                      args=(data.flux[n].values,
                                            data.mask[n].values,
                                            data.v0[n].values,
-                                           data.variance[n].values, 
+                                           data.variance[n].values,
                                            data.guess[n].values,
                                            n, meta, isplots,),
                                      callback=writeBG_WFC3)
@@ -246,7 +246,7 @@ def fitbg(dataim, meta, mask, x1, x2, deg=1, threshold=5, isrotate=False,
                     # stdres = np.std(residuals)
                     # Median Absolute Deviation (slower but more robust)
                     # stdres  = np.median(np.abs(np.ediff1d(residuals)))
-                    # Mean Absolute Deviation (good comprimise)
+                    # Mean Absolute Deviation (good compromise)
                     stdres = np.mean(np.abs(np.ediff1d(residuals)))
                     if stdres == 0:
                         stdres = np.inf
@@ -269,7 +269,7 @@ def fitbg(dataim, meta, mask, x1, x2, deg=1, threshold=5, isrotate=False,
                     plt.title(str(j))
                     plt.plot(goodxvals, dataslice, 'bo')
                     plt.plot(range(nx), bg[j], 'g-')
-                    fname = 'figs'+os.sep+'Fig6_BG_'+str(j)+figure_filetype
+                    fname = 'figs'+os.sep+'Fig3601_BG_'+str(j)+figure_filetype
                     plt.savefig(meta.outputdir + fname, dpi=300)
                     plt.pause(0.01)
 

src/eureka/S3_data_reduction/optspex.py

@@ -1,5 +1,6 @@
 import numpy as np
 import matplotlib.pyplot as plt
+import scipy.interpolate as spi
 from tqdm import tqdm
 from ..lib import gaussian as g
 from ..lib import smooth
@@ -508,7 +509,7 @@ def optimize_wrapper(data, meta, log, apdata, apmask, apbg, apv0, gain=1,
         average smoothing. Defaults to 'hanning'.
     m : int; optional
         File number. Defaults to 0.
-    
+
     Returns
     -------
     data : Xarray Dataset
@@ -521,7 +522,7 @@ def optimize_wrapper(data, meta, log, apdata, apmask, apbg, apv0, gain=1,
     Notes
     -----
     History:
-    
+
     - 2022-07-18, Taylor J Bell
         Added optimize_wrapper to iterate over each frame.
     '''
@@ -543,7 +544,15 @@ def optimize_wrapper(data, meta, log, apdata, apmask, apbg, apv0, gain=1,
 
     # Compute median frame
     data_ma = np.ma.masked_where(apmask == 0, apdata)
-    medflux = np.ma.median(data_ma, axis=0).data
+    medflux = np.ma.median(data_ma, axis=0)
+    # Replace masked pixels through interpolation
+    ny, nx = medflux.shape
+    xx, yy = np.meshgrid(np.arange(nx), np.arange(ny))
+    x1 = xx[~medflux.mask].ravel()
+    y1 = yy[~medflux.mask].ravel()
+    goodmed = medflux[~medflux.mask].ravel()
+    interpmed = spi.griddata((x1, y1), goodmed, (xx, yy),
+                             method='cubic', fill_value=0)
 
     # Perform optimal extraction on each of the frames
     iterfn = range(meta.int_start, meta.n_int)
@@ -558,7 +567,7 @@ def optimize_wrapper(data, meta, log, apdata, apmask, apbg, apv0, gain=1,
                      fittype=meta.fittype,
                      window_len=meta.window_len,
                      deg=meta.prof_deg, windowtype=windowtype,
-                     n=n, m=m, meddata=medflux)
+                     n=n, m=m, meddata=interpmed)
 
     # Mask out NaNs and Infs
     optspec_ma = np.ma.masked_invalid(data.optspec.values)

src/eureka/S3_data_reduction/plots_s3.py

@@ -1,7 +1,9 @@
 import numpy as np
 import os
 from tqdm import tqdm
+import matplotlib as mpl
 import matplotlib.pyplot as plt
+import scipy.interpolate as spi
 from .source_pos import gauss
 from ..lib import util
 from ..lib.plots import figure_filetype
@@ -93,6 +95,15 @@ def image_and_background(data, meta, log, m):
     xmin, xmax = data.flux.x.min().values, data.flux.x.max().values
     ymin, ymax = data.flux.y.min().values, data.flux.y.max().values
 
+    # Commented out vmax calculation is sensitive to unflagged hot pixels
+    # vmax = np.ma.max(np.ma.masked_invalid(subdata))/40
+    vmin = -200
+    vmax = 1000
+    median = np.ma.median(subbg)
+    std = np.ma.std(subbg)
+    # Set bad pixels to plot as black
+    cmap = mpl.cm.get_cmap("plasma").copy()
+    cmap.set_bad('k', 1.)
     iterfn = range(meta.int_end-meta.int_start)
     if meta.verbose:
         iterfn = tqdm(iterfn)
@@ -101,18 +112,16 @@ def image_and_background(data, meta, log, m):
         plt.clf()
         plt.suptitle(f'Integration {intstart + n}')
         plt.subplot(211)
-        plt.title('Background-Subtracted Flux')
-        max = np.ma.max(subdata[n])
-        plt.imshow(subdata[n], origin='lower', aspect='auto',
-                   vmin=0, vmax=max/10, extent=[xmin, xmax, ymin, ymax])
+        plt.title('Background-Subtracted Frame')
+        plt.imshow(subdata[n], origin='lower', aspect='auto', cmap=cmap,
+                   vmin=vmin, vmax=vmax, extent=[xmin, xmax, ymin, ymax])
         plt.colorbar()
         plt.ylabel('Detector Pixel Position')
         plt.subplot(212)
         plt.title('Subtracted Background')
-        median = np.ma.median(subbg[n])
-        std = np.ma.std(subbg[n])
-        plt.imshow(subbg[n], origin='lower', aspect='auto', vmin=median-3*std,
-                   vmax=median+3*std, extent=[xmin, xmax, ymin, ymax])
+        plt.imshow(subbg[n], origin='lower', aspect='auto', cmap=cmap,
+                   vmin=median-3*std, vmax=median+3*std,
+                   extent=[xmin, xmax, ymin, ymax])
         plt.colorbar()
         plt.ylabel('Detector Pixel Position')
         plt.xlabel('Detector Pixel Position')
@@ -426,3 +435,87 @@ def driftywidth(data, meta):
     plt.savefig(meta.outputdir+fname, bbox_inches='tight', dpi=300)
     if not meta.hide_plots:
         plt.pause(0.2)
+
+
+def residualBackground(data, meta, m, vmin=-200, vmax=1000):
+    '''Plot the median, BG-subtracted frame to study the residual BG region and
+    aperture/BG sizes. (Fig 3304)
+
+    Parameters
+    ----------
+    data : Xarray Dataset
+        The Dataset object.
+    meta : eureka.lib.readECF.MetaClass
+        The metadata object.
+    m : int
+        The file number.
+    vmin : int; optional
+        Minimum value of colormap. Default is -200.
+    vmax : int; optional
+        Maximum value of colormap. Default is 1000.
+
+    Returns
+    -------
+    None
+
+    Notes
+    -----
+    History:
+
+    - 2022-07-29 KBS
+        First version
+    '''
+    # Median flux of segment
+    subdata = np.ma.masked_where(~data.mask.values, data.flux.values)
+    flux = np.ma.median(subdata, axis=0)
+    # Compute vertical slice of with 10 columns
+    slice = np.nanmedian(flux[:, meta.subnx//2-5:meta.subnx//2+5], axis=1)
+    # Interpolate to 0.01-pixel resolution
+    f = spi.interp1d(np.arange(meta.subny), slice, 'cubic')
+    ny_hr = np.arange(0, meta.subny-1, 0.01)
+    flux_hr = f(ny_hr)
+    # Set bad pixels to plot as black
+    cmap = mpl.cm.get_cmap("plasma").copy()
+    cmap.set_bad('k', 1.)
+
+    plt.figure(3304, figsize=(8, 3.5))
+    plt.clf()
+    fig, (a0, a1) = plt.subplots(1, 2, gridspec_kw={'width_ratios': [3, 1]},
+                                 num=3304, figsize=(8, 3.5))
+    a0.imshow(flux, origin='lower', aspect='auto', vmax=vmax, vmin=vmin,
+              cmap=cmap)
+    a0.hlines([meta.bg_y1, meta.bg_y2], 0, meta.subnx, color='orange')
+    a0.hlines([meta.src_ypos+meta.spec_hw, meta.src_ypos-meta.spec_hw], 0,
+              meta.subnx, color='mediumseagreen', linestyle='dashed')
+    a0.axes.set_xlim(0, meta.subnx)
+    a0.axes.set_ylabel("Pixel Position")
+    a0.axes.set_xlabel("Pixel Position")
+    a1.scatter(flux_hr, ny_hr, 5, flux_hr, cmap='plasma',
+               norm=plt.Normalize(vmin, vmax))
+    a1.vlines([0], 0, meta.subny, color='0.5', linestyle='dotted')
+    a1.hlines([meta.bg_y1, meta.bg_y2], vmin, vmax, color='orange',
+              linestyle='solid', label='bg'+str(meta.bg_hw))
+    a1.hlines([meta.src_ypos+meta.spec_hw, meta.src_ypos-meta.spec_hw], vmin,
+              vmax, color='mediumseagreen', linestyle='dashed',
+              label='ap'+str(meta.spec_hw))
+    a1.legend(loc='upper right', fontsize=8)
+    a1.axes.set_xlabel("Flux [e-]")
+    a1.axes.set_xlim(vmin, vmax)
+    a1.axes.set_ylim(0, meta.subny)
+    a1.axes.set_yticklabels([])
+    # a1.yaxis.set_visible(False)
+    a1.axes.set_xticks(np.linspace(vmin, vmax, 3))
+    fig.colorbar(plt.cm.ScalarMappable(norm=plt.Normalize(vmin, vmax),
+                 cmap='plasma'), ax=a1)
+    fig.subplots_adjust(top=0.97,
+                        bottom=0.155,
+                        left=0.08,
+                        right=0.925,
+                        hspace=0.2,
+                        wspace=0.08)
+    file_number = str(m).zfill(int(np.floor(np.log10(meta.num_data_files))+1))
+    fname = (f'figs{os.sep}fig3304_file{file_number}' +
+             '_ResidualBG'+figure_filetype)
+    plt.savefig(meta.outputdir+fname, dpi=300)
+    if not meta.hide_plots:
+        plt.pause(0.1)

src/eureka/S3_data_reduction/s3_reduce.py

@@ -258,8 +258,8 @@ def reduce(eventlabel, ecf_path=None, s2_meta=None):
                 data['mask'] = (['time', 'y', 'x'],
                                 np.ones(data.flux.shape, dtype=bool))
 
-                # Check if arrays have NaNs
-                log.writelog('  Masking NaNs in data arrays...',
+                # Check if arrays have NaNs/infs
+                log.writelog('  Masking NaNs/infs in data arrays...',
                              mute=(not meta.verbose))
                 data.mask.values = util.check_nans(data.flux.values,
                                                    data.mask.values,
@@ -271,6 +271,14 @@ def reduce(eventlabel, ecf_path=None, s2_meta=None):
                                                    data.mask.values,
                                                    log, name='V0')
 
+                # Start masking pixels based on DQ flags
+                # https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html
+                # Odd numbers in DQ array are bad pixels. Do not use.
+                if hasattr(meta, 'dqmask') and meta.dqmask:
+                    # dqmask = np.where(data['dq'] > 0)
+                    dqmask = np.where(data.dq % 2 == 1)
+                    data['mask'].values[dqmask] = 0
+
                 # Manually mask regions [colstart, colend, rowstart, rowend]
                 if hasattr(meta, 'manmask'):
                     data = util.manmask(data, meta, log)
@@ -340,6 +348,8 @@ def reduce(eventlabel, ecf_path=None, s2_meta=None):
                     for n in iterfn:
                         # make optimal spectrum plot
                         plots_s3.optimal_spectrum(data, meta, n, m)
+                    if meta.inst != 'wfc3':
+                        plots_s3.residualBackground(data, meta, m)
 
                 if meta.save_output:
                     # Save flux data from current segment

src/eureka/lib/medstddev.py

@@ -118,10 +118,10 @@ def medstddev(data, mask=None, medi=False, axis=0):
 
     # critical case fixes:
     if np.any(ngood == 0):
-        std[np.where(ngood == 0)[0]] = np.nan
-        median[np.where(ngood == 0)[0]] = np.nan
+        std[np.where(ngood == 0)] = np.nan
+        median[np.where(ngood == 0)] = np.nan
     if np.any(ngood == 1):
-        std[np.where(ngood == 1)[0]] = 0.
+        std[np.where(ngood == 1)] = 0.
 
     if len(std) == 1:
         std = std[0]

src/eureka/lib/util.py

@@ -159,15 +159,18 @@ def check_nans(data, mask, log, name=''):
     """
     data = np.ma.masked_where(mask == 0, np.copy(data))
     num_nans = np.sum(np.ma.masked_invalid(data).mask)
+    num_pixels = np.size(data)
+    perc_nans = 100*num_nans/num_pixels
     if num_nans > 0:
-        log.writelog(f"  WARNING: {name} has {num_nans} NaNs/infs. Your "
-                     "subregion may be off the edge of the detector "
-                     "subarray.\n    Masking NaN region and continuing, "
-                     "but you should really stop and reconsider your "
-                     "choices.")
+        log.writelog(f"  {name} has {num_nans} NaNs/infs, which is "
+                     f"{perc_nans:.2f}% of all pixels.")
         inan = np.where(np.ma.masked_invalid(data).mask)
-        # subdata[inan]  = 0
         mask[inan] = 0
+    if perc_nans > 10:
+        log.writelog("  WARNING: Your region of interest may be off the edge "
+                     "of the detector subarray.  Masking NaN/inf regions and "
+                     "continuing, but you should really stop and reconsider "
+                     "your choices.")
     return mask
 
 

tests/NIRCam_ecfs/S3_NIRCam.ecf

@@ -12,13 +12,13 @@ src_pos_type gaussian   # Determine source position when not given in header (Op
 record_ypos True        # Option to record the y position and width for each integration (only records if src_pos_type is gaussian)
 
 # Background parameters
-bg_hw       20          # Half-width of exclusion region for BG subtraction (relative to source position)
+bg_hw       12          # Half-width of exclusion region for BG subtraction (relative to source position)
 bg_thresh   [5,5]       # Double-iteration X-sigma threshold for outlier rejection along time axis
 bg_deg      1           # Polynomial order for column-by-column background subtraction, -1 for median of entire frame
 p3thresh    5           # X-sigma threshold for outlier rejection during background subtraction
 
 # Spectral extraction parameters
-spec_hw     20			# Half-width of aperture region for spectral extraction (relative to source position)
+spec_hw     8			# Half-width of aperture region for spectral extraction (relative to source position)
 fittype     meddata     # Method for constructing spatial profile (Options: smooth, meddata, poly, gauss, wavelet, or wavelet2D)
 window_len  31          # Smoothing window length, when fittype = smooth
 prof_deg    3           # Polynomial degree, when fittype = poly

tests/test_NIRCam.py

@@ -53,28 +53,28 @@ def test_NIRCam(capsys):
     # run assertions for S3
     meta.outputdir_raw = (f'data{os.sep}JWST-Sim{os.sep}NIRCam{os.sep}'
                           f'Stage3{os.sep}')
-    name = pathdirectory(meta, 'S3', 1, ap=20, bg=20)
+    name = pathdirectory(meta, 'S3', 1, ap=8, bg=12)
     assert os.path.exists(name)
     assert os.path.exists(name+os.sep+'figs')
 
     # run assertions for S4
     meta.outputdir_raw = (f'data{os.sep}JWST-Sim{os.sep}NIRCam{os.sep}'
                           f'Stage4{os.sep}')
-    name = pathdirectory(meta, 'S4', 1, ap=20, bg=20)
+    name = pathdirectory(meta, 'S4', 1, ap=8, bg=12)
     assert os.path.exists(name)
     assert os.path.exists(name+os.sep+'figs')
 
     # run assertions for S5
     meta.outputdir_raw = (f'data{os.sep}JWST-Sim{os.sep}NIRCam{os.sep}'
                           f'Stage5{os.sep}')
-    name = pathdirectory(meta, 'S5', 1, ap=20, bg=20)
+    name = pathdirectory(meta, 'S5', 1, ap=8, bg=12)
     assert os.path.exists(name)
     assert os.path.exists(name+os.sep+'figs')
 
     # run assertions for S6
     meta.outputdir_raw = (f'data{os.sep}JWST-Sim{os.sep}NIRCam{os.sep}'
                           f'Stage6{os.sep}')
-    name = pathdirectory(meta, 'S6', 1, ap=20, bg=20)
+    name = pathdirectory(meta, 'S6', 1, ap=8, bg=12)
     assert os.path.exists(name)
     assert os.path.exists(name+os.sep+'figs')