JP-3730: Re-enable saving of blot models during outlier detection (sp…

…acetelescope#8758)

CHANGES.rst

@@ -496,6 +496,8 @@ outlier_detection
 
 - Fix errors in documentation describing arguments. [#8603]
 
+- Re-enabled saving of blot models when `save_intermediate_results` is True. [#8758]
+
 pathloss
 --------
 

docs/jwst/outlier_detection/outlier_detection_imaging.rst

@@ -75,7 +75,8 @@ Specifically, this routine performs the following operations:
    * The median image is written out to disk as `_<asn_id>_median.fits` by default.
 
 #. By default, the median image is blotted back (inverse of resampling) to
-   match each original input image.
+   match each original input image.  Resampled/blotted images are written out to disk if 
+   the ``save_intermediate_results`` parameter is set to `True`.
 
    * **If resampling is turned off**, the median image is compared directly to
      each input image.

docs/jwst/outlier_detection/outlier_detection_spec.rst

@@ -39,6 +39,8 @@ Specifically, this routine performs the following operations (modified from the
      if the ``save_intermediate_results`` parameter is set to `True`
 #. Blot median image to match each original input image
 
+   - Resampled/blotted images are written out to disk if the ``save_intermediate_results``
+     parameter is set to `True`
    - **If resampling is turned off**, the median image is used for comparison
      with the original input models for detecting outliers
 #. Perform statistical comparison between blotted image and original image to identify outliers

jwst/outlier_detection/imaging.py

@@ -127,7 +127,16 @@ def detect_outliers(
     with input_models:
         for image in input_models:
             if resample_data:
-                flag_resampled_model_crs(image, median_data, median_wcs, snr1, snr2, scale1, scale2, backg)
+                flag_resampled_model_crs(image,
+                                         median_data,
+                                         median_wcs,
+                                         snr1,
+                                         snr2,
+                                         scale1,
+                                         scale2,
+                                         backg,
+                                         save_blot=save_intermediate_results,
+                                         make_output_path=make_output_path)
             else:
                 flag_model_crs(image, median_data, snr1)
             input_models.shelve(image, modify=True)

jwst/outlier_detection/spec.py

@@ -130,7 +130,11 @@ def detect_outliers(
             scale1,
             scale2,
             backg,
+            save_blot=save_intermediate_results,
+            make_output_path=make_output_path
         )
     else:
-        flag_crs_in_models(input_models, median_data, snr1)
+        flag_crs_in_models(input_models,
+                           median_data,
+                           snr1)
     return input_models

jwst/outlier_detection/tests/test_outlier_detection.py

@@ -288,11 +288,13 @@ def test_outlier_step_spec(tmp_cwd, tmp_path):
         i2d_files = glob(os.path.join(dirname, '*i2d*.fits'))
         s2d_files = glob(os.path.join(dirname, '*outlier_s2d.fits'))
         median_files = glob(os.path.join(dirname, '*median.fits'))
+        blot_files = glob(os.path.join(dirname, '*blot.fits'))
 
         # intermediate files are removed
         assert len(i2d_files) == 0
         assert len(s2d_files) == 0
         assert len(median_files) == 0
+        assert len(blot_files) == 0
 
         # result files are written to the output directory
         if dirname == output_dir:
@@ -320,25 +322,31 @@ def test_outlier_step_spec(tmp_cwd, tmp_path):
         i2d_files = glob(os.path.join(dirname, '*i2d*.fits'))
         s2d_files = glob(os.path.join(dirname, '*outlier_s2d.fits'))
         median_files = glob(os.path.join(dirname, '*median.fits'))
+        blot_files = glob(os.path.join(dirname, '*blot.fits'))
         if dirname == output_dir:
             # result files are written to the output directory
             assert len(result_files) == len(container)
 
-            # s2d and median files are written to the output directory
+            # s2d, median, and blot files are written to the output directory
             assert len(s2d_files) == len(container)
+            assert len(blot_files) == len(container)
             assert len(median_files) == 1
 
             # i2d files not written
             assert len(i2d_files) == 0
 
             # nothing else was written
-            assert len(all_files) == len(s2d_files) + len(median_files) + len(result_files)
+            assert len(all_files) == len(s2d_files) + \
+                                     len(median_files) + \
+                                     len(result_files) + \
+                                     len(blot_files)
         else:
             # nothing should be written to the current directory
             assert len(result_files) == 0
             assert len(s2d_files) == 0
             assert len(median_files) == 0
             assert len(i2d_files) == 0
+            assert len(blot_files) == 0
             assert len(all_files) == 0
 
     miri_rate.close()
@@ -405,6 +413,30 @@ def test_outlier_step_on_disk(three_sci_as_asn, tmp_cwd):
         assert zeroth.dq[12, 12] == OUTLIER_DO_NOT_USE
         result.shelve(zeroth, modify=False)
 
+    # Verify intermediate results were written to disk
+    dirname = tmp_cwd
+    all_files = glob(os.path.join(dirname, '*.fits'))
+    input_files = glob(os.path.join(dirname, '*_cal.fits'))
+    result_files = glob(os.path.join(dirname, '*outlierdetectionstep.fits'))
+    i2d_files = glob(os.path.join(dirname, '*i2d*.fits'))
+    s2d_files = glob(os.path.join(dirname, '*outlier_s2d.fits'))
+    median_files = glob(os.path.join(dirname, '*median.fits'))
+    blot_files = glob(os.path.join(dirname, '*blot.fits'))
+
+    assert len(result_files) == len(container)
+
+    # i2d, median, blot files are written to the output directory
+    assert len(i2d_files) == len(container)
+    assert len(blot_files) == len(container)
+    assert len(median_files) == 1
+
+    # s2d files not written
+    assert len(s2d_files) == 0
+
+    # nothing else was written
+    assert len(all_files) == len(input_files) + len(i2d_files) + len(median_files) + len(result_files) + len(blot_files)
+
+
 
 def test_outlier_step_square_source_no_outliers(we_three_sci, tmp_cwd):
     """Test whole step with square source with sharp edges, no outliers"""

jwst/outlier_detection/utils.py

@@ -206,6 +206,8 @@ def flag_resampled_model_crs(
     scale1,
     scale2,
     backg,
+    save_blot=False,
+    make_output_path=None,
 ):
     if 'SPECTRAL' not in input_model.meta.wcs.output_frame.axes_type:
         input_pixflux_area = input_model.meta.photometry.pixelarea_steradians
@@ -217,6 +219,15 @@ def flag_resampled_model_crs(
         pix_ratio = 1.0
 
     blot = gwcs_blot(median_data, median_wcs, input_model.data.shape, input_model.meta.wcs, pix_ratio)
+    if save_blot:
+        if make_output_path is None:
+            raise ValueError("make_output_path must be provided if save_blot is True")
+        model_path = make_output_path(input_model.meta.filename, suffix='blot')
+        blot_model = _make_blot_model(input_model, blot)
+        blot_model.meta.filename = model_path
+        blot_model.save(model_path)
+        log.info(f"Saved model in {model_path}")
+        del blot_model
     # dq flags will be updated in-place
     _flag_resampled_model_crs(input_model, blot, snr1, snr2, scale1, scale2, backg)
 
@@ -260,9 +271,20 @@ def flag_crs_in_models_with_resampling(
     scale1,
     scale2,
     backg,
+    save_blot=False,
+    make_output_path=None,
 ):
     for image in input_models:
-        flag_resampled_model_crs(image, median_data, median_wcs, snr1, snr2, scale1, scale2, backg)
+        flag_resampled_model_crs(image,
+                                 median_data,
+                                 median_wcs,
+                                 snr1,
+                                 snr2,
+                                 scale1,
+                                 scale2,
+                                 backg,
+                                 save_blot=save_blot,
+                                 make_output_path=make_output_path)
 
 
 def flag_model_crs(image, blot, snr):
@@ -276,3 +298,10 @@ def flag_model_crs(image, blot, snr):
     match_nans_and_flags(image)
 
     log.info(f"{np.count_nonzero(cr_mask)} pixels marked as outliers")
+
+
+def _make_blot_model(input_model, blot):
+    blot_model = type(input_model)()
+    blot_model.data = blot
+    blot_model.update(input_model)
+    return blot_model