Merge pull request #40 from radicamc/dev

Dev

README.md

@@ -1,6 +1,8 @@
 # supreme-SPOON
 supreme-**S**teps to **P**rocess S**O**SS **O**bservatio**N**s
 
+⭐⭐ **Please note that supreme-SPOON will be undergoing a rebrand and future versions (>=1.4.0) will use a new name, exoTEDRF.** ⭐⭐
+
 **supreme-SPOON** is an end-to-end pipeline for NIRISS/SOSS time series observations (TSOs).
 The pipeline is divided into four stages:
  - Stage 1: Detector Level Processing 

changelog.md

@@ -1,6 +1,19 @@
 # Changelog
 All notable changes to this project will be documented in this file.
 
+### [1.3.2] -- 2024-03-08
+#### Added
+- Allow the trace mask width to be specified directly during 1/f correction.
+- More support for SUBSTRIP96, particularly in plotting routines. 
+- Add possibility for flux offset in background correction.
+- Best-fitting models are now automatically saved in Stage 4. 
+- Patch of ATOCA implementation to circumvent occasional error in estimation of wavelength grid.
+- supreme-SPOON end-of-life notices.
+- Misc. bug fixes. 
+
+#### Removed
+- Trace masking functionality in TracingStep.
+
 ### [1.3.1] -- 2024-01-19
 #### Added
 - Fixed implementation of STScI up-the-ramp JumpStep which would not run consistently.

setup.py

@@ -4,7 +4,7 @@
 from setuptools import setup
 
 setup(name='supreme_spoon',
-      version='1.3.1',
+      version='1.3.2',
       license='MIT',
       author='Michael Radica',
       author_email='[email protected]',

supreme_spoon/fit_lightcurves.py

@@ -267,7 +267,7 @@
     # make summary plots if necessary.
     fancyprint('Summarizing fit results.')
     data = np.ones((nints, nbins)) * np.nan
-    models = np.ones((3, nints, nbins)) * np.nan
+    models = np.ones((4, nints, nbins)) * np.nan
     residuals = np.ones((nints, nbins)) * np.nan
     order_results = {'dppm': [], 'dppm_err': [], 'wave': wave,
                      'wave_err': np.mean([wave - wave_low, wave_up - wave],
@@ -366,6 +366,7 @@
                 models[1, :, i] = systematics
             if gp_model is not None:
                 models[2, :, i] = gp_model
+            models[3, :, i] = norm_flux[:, i]
             residuals[:, i] = norm_flux[:, i] - transit_model
 
     results_dict['order {}'.format(order)] = order_results

supreme_spoon/plotting.py

@@ -47,6 +47,14 @@ def make_background_row_plot(before, after, background_model, row_start=230,
     else:
         bkg = background_model
 
+    # If SUBSTRIP96, change rows to use.
+    if np.shape(af)[-2] == 96:
+        sub96 = True
+        row_start = 5
+        row_end = 21
+    else:
+        sub96 = False
+
     # Create medians.
     if np.ndim(af) == 4:
         before = bn.nanmedian(bf[:, -1], axis=0)
@@ -70,8 +78,12 @@ def make_background_row_plot(before, after, background_model, row_start=230,
 
     plt.axvline(700, ls=':', c='grey')
     plt.axhline(0, ls=':', c='grey')
-    plt.ylim(np.min([np.nanmin(aafter), np.nanmin(bbefore[700:] - bkg_scale)]),
-             np.nanpercentile(bbefore, 95))
+    if sub96 is True:
+        plt.ylim(np.min([np.nanmin(aafter), np.nanmin(bbefore[700:] - bkg_scale)]),
+                 np.nanpercentile(bbefore, 75))
+    else:
+        plt.ylim(np.min([np.nanmin(aafter), np.nanmin(bbefore[700:] - bkg_scale)]),
+                 np.nanpercentile(bbefore, 95))
     plt.xlabel('Spectral Pixel', fontsize=12)
     plt.ylabel('Counts', fontsize=12)
 
@@ -399,7 +411,12 @@ def make_jump_location_plot(results, outfile=None, show_plot=True):
                                      fill=False)
                     ax.add_patch(marker)
 
-            ax.text(30, 230, '({0}, {1})'.format(i, g), c='white', fontsize=12)
+            if dimy == 96:
+                ax.text(30, 80, '({0}, {1})'.format(i, g), c='white',
+                        fontsize=12)
+            else:
+                ax.text(30, 230, '({0}, {1})'.format(i, g), c='white',
+                        fontsize=12)
             if j != 0:
                 ax.yaxis.set_major_formatter(plt.NullFormatter())
             else:
@@ -569,7 +586,7 @@ def make_linearity_plot(results, old_results, outfile=None, show_plot=True):
     """Plot group differences before and after linearity correction.
     """
 
-    fancyprint('Doing diagnostic plot.')
+    fancyprint('Doing diagnostic plot 1.')
     results = np.atleast_1d(results)
     old_results = np.atleast_1d(old_results)
     for i, file in enumerate(results):
@@ -635,6 +652,67 @@ def make_linearity_plot(results, old_results, outfile=None, show_plot=True):
         plt.show()
 
 
+def make_linearity_plot2(results, old_results, outfile=None, show_plot=True):
+    """Plot residuals to flat line before and after linearity correction.
+    """
+
+    fancyprint('Doing diagnostic plot 2.')
+    results = np.atleast_1d(results)
+    old_results = np.atleast_1d(old_results)
+    for i, file in enumerate(results):
+        with utils.open_filetype(file) as datamodel:
+            if i == 0:
+                cube = datamodel.data
+            else:
+                cube = np.concatenate([cube, datamodel.data])
+    for i, file in enumerate(old_results):
+        with utils.open_filetype(file) as datamodel:
+            if i == 0:
+                old_cube = datamodel.data
+            else:
+                old_cube = np.concatenate([old_cube, datamodel.data])
+
+    nint, ngroup, dimy, dimx = np.shape(cube)
+    # Get bright pixels in the trace.
+    stack = bn.nanmedian(cube[:, -1], axis=0)
+    ii = np.where((stack >= np.nanpercentile(stack, 80)) &
+                  (stack < np.nanpercentile(stack, 99)))
+    jj = np.random.randint(0, len(ii[0]), 1000)
+    y = ii[0][jj]
+    x = ii[1][jj]
+    i = np.random.randint(0, nint, 1000)
+
+    oold = np.zeros((1000, ngroup))
+    nnew = np.zeros((1000, ngroup))
+    for j in range(1000):
+        o = old_cube[i[j], :, y[j], x[j]]
+        ol = np.linspace(np.min(o), np.max(o), ngroup)
+        oold[j] = (o - ol) / np.max(o) * 100
+        n = cube[i[j], :, y[j], x[j]]
+        nl = np.linspace(np.min(n), np.max(n), ngroup)
+        nnew[j] = (n - nl) / np.max(n) * 100
+
+    # Plot up mean group differences before and after linearity correction.
+    plt.figure(figsize=(5, 3))
+    plt.plot(np.arange(ngroup)+1, np.nanmedian(oold, axis=0),
+             label='Before Correction', c='blue')
+    plt.plot(np.arange(ngroup)+1, np.nanmedian(nnew, axis=0),
+             label='After Correction', c='red')
+    plt.axhline(0, ls='--', c='black')
+    plt.xticks(np.arange(ngroup)+1, (np.arange(ngroup)+1).astype(str))
+    plt.xlabel('Group Number', fontsize=12)
+    plt.ylabel('Residual [%]', fontsize=12)
+    plt.legend()
+
+    if outfile is not None:
+        plt.savefig(outfile, bbox_inches='tight')
+        fancyprint('Plot saved to {}'.format(outfile))
+    if show_plot is False:
+        plt.close()
+    else:
+        plt.show()
+
+
 def make_oneoverf_chromatic_plot(m_e, m_o, b_e, b_o, ngroup, outfile=None,
                                  show_plot=True):
     """Make plot of chromatic 1/f slope and intercept values.
@@ -742,7 +820,7 @@ def make_oneoverf_plot(results, baseline_ints, timeseries=None,
 
 def make_oneoverf_psd(results, old_results, timeseries, baseline_ints,
                       nsample=25,  pixel_masks=None, tframe=5.494, tpix=1e-5,
-                      tgap=1.2e-4, outfile=None, show_plot=True, window=False):
+                      tgap=1.2e-4, outfile=None, show_plot=True):
     """Make a PSD plot to see PSD of background before and after 1/f removal.
     """
 
@@ -765,22 +843,12 @@ def make_oneoverf_psd(results, old_results, timeseries, baseline_ints,
                 old_cube = np.concatenate([old_cube, datamodel.data])
     old_cube = np.where(np.isnan(old_cube), np.nanmedian(old_cube), old_cube)
     if pixel_masks is not None:
-        if window is True:
-            for i, file in enumerate(pixel_masks):
-                mask_in = fits.getdata(file, 3)
-                mask_out = fits.getdata(file, 5)
-                window = ~(mask_out - mask_in).astype(bool)
-                if i == 0:
-                    mask_cube = window
-                else:
-                    mask_cube = np.concatenate([mask_cube, window])
-        else:
-            for i, file in enumerate(pixel_masks):
-                data = fits.getdata(file, 1)
-                if i == 0:
-                    mask_cube = data
-                else:
-                    mask_cube = np.concatenate([mask_cube, data])
+        for i, file in enumerate(pixel_masks):
+            data = fits.getdata(file)
+            if i == 0:
+                mask_cube = data
+            else:
+                mask_cube = np.concatenate([mask_cube, data])
 
     else:
         mask_cube = None
@@ -937,7 +1005,7 @@ def make_photon_noise_plot(spectrum_files, ngroup, baseline_ints, order=1,
             if order == 1:
                 wave = np.mean([spectrum[1].data[0], spectrum[2].data[0]],
                                axis=0)
-                ii = np.ones(2040)
+                ii = np.ones_like(wave)
             else:
                 wave = np.mean([spectrum[5].data[0], spectrum[6].data[0]],
                                axis=0)

supreme_spoon/run_DMS.py

@@ -11,7 +11,6 @@
 from datetime import datetime
 import numpy as np
 import os
-import pandas as pd
 import shutil
 import sys
 
@@ -67,14 +66,10 @@
 
 # Open some of the input files.
 background_model = np.load(config['background_file'])
-if np.ndim(background_model) == 3:
-    background_model = background_model[0]
 if config['timeseries'] is not None:
     config['timeseries'] = np.load(config['timeseries'])
 if config['timeseries_o2'] is not None:
     config['timeseries_o2'] = np.load(config['timeseries_o2'])
-if config['centroids'] is not None:
-    config['centroids'] = pd.read_csv(config['centroids'], comment='#')
 if config['f277w'] is not None:
     config['f277w'] = np.load(config['f277w'])
 
@@ -95,6 +90,8 @@
     stage1_results = run_stage1(input_files, background_model=background_model,
                                 baseline_ints=config['baseline_ints'],
                                 oof_method=config['oof_method'],
+                                inner_mask_width=config['inner_mask_width'],
+                                outer_mask_width=config['outer_mask_width'],
                                 save_results=config['save_results'],
                                 pixel_masks=config['outlier_maps'],
                                 force_redo=config['force_redo'],
@@ -107,6 +104,7 @@
                                 do_plot=config['do_plots'],
                                 timeseries=config['timeseries'],
                                 timeseries_o2=config['timeseries_o2'],
+                                centroids=config['centroids'],
                                 **config['stage1_kwargs'])
 else:
     stage1_results = input_files
@@ -141,9 +139,9 @@
                                 smoothing_scale=config['smoothing_scale'],
                                 skip_steps=stage2_skip,
                                 generate_lc=config['generate_lc'],
-                                generate_tracemask=config['generate_tracemask'],
                                 inner_mask_width=config['inner_mask_width'],
                                 outer_mask_width=config['outer_mask_width'],
+                                centroids=config['centroids'],
                                 pixel_masks=config['outlier_maps'],
                                 generate_order0_mask=config['generate_order0_mask'],
                                 f277w=config['f277w'],

supreme_spoon/run_DMS.yaml

@@ -33,6 +33,10 @@ timeseries : None
 timeseries_o2 : None
 # For 1/f correction; outlier pixel maps (optional).
 outlier_maps : None
+# Width, in pixels, of inner trace mask.
+inner_mask_width : 40
+# Width, in pixels, of outer trace mask.
+outer_mask_width : 70
 # For jump detection; whether to flag outliers up-the-ramp (jwst default).
 flag_up_ramp : False
 # For jump detection; sigma threshold for a pixel to be flagged as an outlier.
@@ -57,12 +61,6 @@ TracingStep : 'run'
 space_outlier_threshold : 15
 # For outlier detection; sigma threshold for a pixel to be flagged as an outlier, temporally.
 time_outlier_threshold : 10
-# If True, create a mask of the target trace orders.
-generate_tracemask : True
-# Width, in pixels, of inner trace mask. Only necessary if generate_tracemask is True.
-inner_mask_width : 40
-# Width, in pixels, of outer trace mask. Only necessary if generate_tracemask is True.
-outer_mask_width : 70
 # If True, calculate the stability of the SOSS trace over the course of the TSO using principle component analysis.
 calculate_stability : True
 # Number of principle components to extract. Only necessary if calculate_stability is True.
@@ -84,11 +82,9 @@ stage2_kwargs : {}
 extract_method : 'box'
 # Box width to extract around the trace center. Passing 'optimize' will find and use the aperture width that minimizes
 # the white light curve scatter.
-soss_width : 40
+soss_width : 35
 # Specprofile reference file for ATOCA (optional).
 specprofile : None
-# File with centroids for all three orders for box extraction (optional).
-centroids : None
 # Stellar effective temperature in K for wavelength calibration (optional).
 st_teff : None
 # Stellar log surface gravity for wavelength calibration (optional).
@@ -113,6 +109,8 @@ save_results : True
 force_redo : False
 # Integrations of ingress and/or egress.
 baseline_ints : [50, -50]
+# File with centroids for all three orders for box extraction and 1/f correction (optional).
+centroids : None
 # If True, make diagnostic plots.
 do_plots : True
 # ========================= END DMS CONFIG FILE =============================