JP-824: Add option to save combined background image (#5954)

* added option to save combined background

* fix hanging call from previous method

* TODO: removed commented code if review is okay

* TODO: remove commented code if review is okay

* TODO: remove commented code // flake8

* removed commented code, pull data from masked array

* remove comment, add regtest update for saved bkgd

* docs and changelog

* turn on resample in regtest

CHANGES.rst

@@ -43,6 +43,9 @@ background
 
 - Remove unused ``SubtractImagesStep`` [#5919]
 
+- Added new step parameter to optionally save the combined, average
+  background image: ``save_combined_background``. [#5954]
+
 calwebb_spec2
 -------------
 

docs/jwst/background_step/arguments.rst

@@ -16,3 +16,7 @@ control the sigma clipping, and are passed as arguments to the astropy
   The number of clipping iterations to perform, or ``None`` to clip until
   convergence is achieved. Defaults to ``None``.
 
+``--save_combined_background``
+  Saves the combined, averaged background image used for background
+  subtraction. Defaults to ``False``.
+

docs/jwst/background_step/description.rst

@@ -53,6 +53,9 @@ subtracted from it. If the target exposure is in the form of a 3-D CubeModel
 (e.g. the result of a time series exposure), the background image
 is subtracted from each plane of the CubeModel.
 
+The combined, averaged background image can be saved using the step parameter
+``save_combined_background``.
+
 WFSS Mode
 ---------
 For Wide-Field Slitless Spectroscopy expsoures (NIS_WFSS and NRC_WFSS),

jwst/background/background_step.py

@@ -13,13 +13,17 @@ class BackgroundStep(Step):
     """
 
     spec = """
+        save_combined_background = boolean(default=False)  # Save combined background image
         sigma = float(default=3.0)  # Clipping threshold
         maxiters = integer(default=None)  # Number of clipping iterations
     """
 
     # These reference files are only used for WFSS/GRISM data.
     reference_file_types = ["wfssbkg", "wavelengthrange"]
 
+    # Define a suffix for optional saved output of the combined background
+    bkg_suffix = 'combinedbackground'
+
     def process(self, input, bkg_list):
         """
         Subtract the background signal from target exposures by subtracting
@@ -77,11 +81,15 @@ def process(self, input, bkg_list):
                                 break
                 # Do the background subtraction
                 if do_sub:
-                    result = background_sub.background_sub(input_model,
-                                                           bkg_list,
-                                                           self.sigma,
-                                                           self.maxiters)
+                    bkg_model, result = background_sub.background_sub(input_model,
+                                                                      bkg_list,
+                                                                      self.sigma,
+                                                                      self.maxiters)
                     result.meta.cal_step.back_sub = 'COMPLETE'
+                    if self.save_combined_background:
+                        comb_bkg_path = self.save_model(bkg_model, suffix=self.bkg_suffix, force=True)
+                        self.log.info(f'Combined background written to "{comb_bkg_path}".')
+
                 else:
                     result = input_model.copy()
                     result.meta.cal_step.back_sub = 'SKIPPED'

jwst/background/background_sub.py

@@ -41,11 +41,8 @@ def background_sub(input_model, bkg_list, sigma, maxiters):
     log.debug(' subtracting avg bkg from {}'.format(input_model.meta.filename))
     result = subtract_images.subtract(input_model, bkg_model)
 
-    # Close the average background image and update the step status
-    bkg_model.close()
-
     # We're done. Return the result.
-    return result
+    return bkg_model, result
 
 
 def average_background(bkg_list, sigma, maxiters):
@@ -101,7 +98,7 @@ def average_background(bkg_list, sigma, maxiters):
     merr = np.ma.masked_array(cerr, mask=mdata.mask)
 
     # Compute the combined ERR as the uncertainty in the mean
-    avg_bkg.err = np.sqrt(merr.sum(axis=0)) / (num_bkg - merr.mask.sum(axis=0))
+    avg_bkg.err = (np.sqrt(merr.sum(axis=0)) / (num_bkg - merr.mask.sum(axis=0))).data
 
     return avg_bkg
 

jwst/regtest/test_miri_lrs_slit_spec2.py

@@ -26,19 +26,19 @@ def run_pipeline(jail, rtdata_module):
     # NOTE: THE RESAMPLE_SPEC STEP IS SKIPPED FOR NOW, BECAUSE IT HAS A BUG
     # (the s2d image is all zeros)
     args = ["config/calwebb_spec2.cfg", rtdata.input,
-            "--steps.resample_spec.skip=true",       # remove when bug fixed
             "--save_bsub=true",
             "--steps.assign_wcs.save_results=true",
             "--steps.flat_field.save_results=true",
-            "--steps.srctype.save_results=true"]
+            "--steps.srctype.save_results=true",
+            "--steps.bkg_subtract.save_combined_background=true"]
     Step.from_cmdline(args)
 
     return rtdata
 
 
 @pytest.mark.bigdata
 @pytest.mark.parametrize("output", [
-    "bsub", "flat_field", "assign_wcs", "srctype", "cal", "x1d"])
+    "bsub", "flat_field", "assign_wcs", "srctype", "cal", "x1d", "combinedbackground"])
 def test_miri_lrs_slit_spec2(run_pipeline, fitsdiff_default_kwargs, output):
     """Regression test of the calwebb_spec2 pipeline on MIRI
        LRS fixedslit data using along-slit-nod pattern for