JP-1649: Add option for user-defined sky levels in skymatch (#9053)

changes/9053.skymatch.rst

@@ -0,0 +1 @@
+Add option to pass in user-defined sky levels

docs/jwst/skymatch/arguments.rst

@@ -8,7 +8,7 @@ The ``skymatch`` step uses the following optional arguments:
 
 ``skymethod`` (str, default='match')
   The sky computation algorithm to be used.
-  Allowed values: `local`, `global`, `match`, `global+match`
+  Allowed values: `local`, `global`, `match`, `global+match`, `user`
 
 ``match_down`` (boolean, default=True)
   Specifies whether the sky *differences* should be subtracted from images with
@@ -23,6 +23,12 @@ The ``skymatch`` step uses the following optional arguments:
   the images. The BKGSUB keyword (boolean) will be set in each output image to
   record whether or not the background was subtracted.
 
+``skylist`` (string, default=None)
+  A filename pointing to a two-column whitespace-delimited list of user-defined
+  (filename, skyval) pairs to be used for sky subtraction. The list
+  must have the same length as the input images and contain exactly one line per
+  image. This argument is used only when ``skymethod`` is set to `user`.
+
 **Image bounding polygon parameters:**
 
 ``stepsize`` (int, default=None)

docs/jwst/skymatch/description.rst

@@ -117,6 +117,17 @@ of sky in the images. This method cannot measure the true sky level, but
 instead provides additive corrections that can be used to equalize the signal
 between overlapping images.
 
+User-Supplied Sky Values
+-------------------------
+The ``skymatch`` step can also accept user-supplied sky values for each image.
+This is useful when sky values have been determined based on a custom workflow
+outside the pipeline. To use this feature, the user must provide a list of sky
+values matching the number of images (``skylist`` parameter) and set the
+``skymethod`` parameter to "user". The ``skylist`` must be a two-column 
+whitespace-delimited file with the first column containing the image filenames
+and the second column containing the sky values. There must be exactly one line
+per image in the input list.
+
 Examples
 --------
 To get a better idea of the behavior of these different methods, the tables below

jwst/skymatch/skymatch.py

@@ -12,7 +12,7 @@
 from . skyimage import SkyImage, SkyGroup
 
 
-__all__ = ['match']
+__all__ = ['skymatch']
 
 
 __author__ = 'Mihai Cara'
@@ -25,7 +25,7 @@
 log.setLevel(logging.DEBUG)
 
 
-def match(images, skymethod='global+match', match_down=True, subtract=False):
+def skymatch(images, skymethod='global+match', match_down=True, subtract=False):
     """
     A function to compute and/or "equalize" sky background in input images.
 
@@ -94,7 +94,6 @@ def match(images, skymethod='global+match', match_down=True, subtract=False):
     subtract : bool (Default = False)
         Subtract computed sky value from image data.
 
-
     Raises
     ------
 
@@ -233,7 +232,7 @@ def match(images, skymethod='global+match', match_down=True, subtract=False):
         in sky levels.
 
     """
-    function_name = match.__name__
+    function_name = skymatch.__name__
 
     # Time it
     runtime_begin = datetime.now()

jwst/skymatch/skymatch_step.py

@@ -20,14 +20,19 @@
 from stdatamodels.jwst.datamodels.dqflags import pixel
 
 from jwst.datamodels import ModelLibrary
+from jwst.lib.suffix import remove_suffix
+from pathlib import Path
 
-from ..stpipe import Step
+from jwst.stpipe import Step
 
 # LOCAL:
-from .skymatch import match
+from .skymatch import skymatch
 from .skyimage import SkyImage, SkyGroup
 from .skystatistics import SkyStats
 
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
 
 __all__ = ['SkyMatchStep']
 
@@ -42,9 +47,10 @@ class SkyMatchStep(Step):
 
     spec = """
         # General sky matching parameters:
-        skymethod = option('local', 'global', 'match', 'global+match', default='match') # sky computation method
+        skymethod = option('local', 'global', 'match', 'global+match', 'user', default='match') # sky computation method
         match_down = boolean(default=True) # adjust sky to lowest measured value?
         subtract = boolean(default=False) # subtract computed sky from image data?
+        skylist = string(default=None) # Filename pointing to list of (imagename skyval) pairs
 
         # Image's bounding polygon parameters:
         stepsize = integer(default=None) # Max vertex separation
@@ -68,13 +74,17 @@ class SkyMatchStep(Step):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
-    def process(self, input):
+    def process(self, input_models):
         self.log.setLevel(logging.DEBUG)
 
-        if isinstance(input, ModelLibrary):
-            library = input
+        if isinstance(input_models, ModelLibrary):
+            library = input_models
         else:
-            library = ModelLibrary(input, on_disk=not self.in_memory)
+            library = ModelLibrary(input_models, on_disk=not self.in_memory)
+
+        # Method: "user". Use user-provided sky values, and bypass skymatch() altogether.
+        if self.skymethod == 'user':
+            return self._user_sky(library)
 
         self._dqbits = interpret_bit_flags(self.dqbits, flag_name_map=pixel)
 
@@ -105,7 +115,7 @@ def process(self, input):
                     images.append(SkyGroup(sky_images, id=group_index))
 
         # match/compute sky values:
-        match(images, skymethod=self.skymethod, match_down=self.match_down,
+        skymatch(images, skymethod=self.skymethod, match_down=self.match_down,
               subtract=self.subtract)
 
         # set sky background value in each image's meta:
@@ -216,3 +226,50 @@ def _set_sky_background(self, sky_image, library, step_status):
 
         dm.meta.cal_step.skymatch = step_status
         library.shelve(dm, index)
+
+
+    def _user_sky(self, library):
+        """Handle user-provided sky values for each image.
+        """
+
+        if self.skylist is None:
+            raise ValueError('skymethod set to "user", but no sky value file provided.')
+
+        log.info(" ")
+        log.info("Setting sky background of input images to user-provided values "
+                 f"from `skylist` ({self.skylist}).")
+
+        # read the comma separated file and get just the stem of the filename
+        skylist = np.genfromtxt(
+            self.skylist,
+            dtype=[("fname", "<S128"),  ("sky", "f")],
+        )
+        skyfnames, skyvals = skylist['fname'], skylist['sky']
+        skyfnames = skyfnames.astype(str)
+        skyfnames = [remove_suffix(Path(fname).stem)[0] for fname in skyfnames]
+        skyfnames = np.array(skyfnames)
+
+        if len(skyvals) != len(library):
+            raise ValueError(f"Number of entries in skylist ({len(self.skylist)}) does not match "
+                             f"number of input images ({len(library)}).")
+
+        with library:
+            for model in library:
+                fname, _ = remove_suffix(Path(model.meta.filename).stem)
+                sky = skyvals[np.where(skyfnames == fname)]
+                if len(sky) == 0:
+                    raise ValueError(f"Image with stem '{fname}' not found in the skylist.")
+                if len(sky) > 1:
+                    raise ValueError(f"Image with stem '{fname}' found multiple times in the skylist.")
+
+                log.debug(f"Setting sky background of image '{model.meta.filename}' to {float(sky)}.")
+
+                model.meta.background.level = float(sky)
+                model.meta.background.subtracted = self.subtract
+                model.meta.background.method = self.skymethod
+                if self.subtract:
+                    model.data -= sky
+                model.meta.cal_step.skymatch = "COMPLETE"
+                library.shelve(model)
+
+        return library

jwst/skymatch/skystatistics.py

@@ -1,6 +1,6 @@
 """
 `skystatistics` module provides statistics computation class used by
-:py:func:`~jwst.skymatch.skymatch.match`
+:py:func:`~jwst.skymatch.skymatch.skymatch`
 and :py:class:`~jwst.skymatch.skyimage.SkyImage`.
 
 :Authors: Mihai Cara (contact: [email protected])

jwst/skymatch/tests/test_skymatch.py

@@ -547,3 +547,96 @@ def test_skymatch_2x(tmp_cwd, nircam_rate, tmp_path, skymethod, subtract):
             else:
                 assert abs(np.mean(im2.data[dq_mask]) - lev) < 0.01
             result2.shelve(im2)
+
+
+@pytest.mark.parametrize("subtract", (False, True))
+def test_user_skyfile(tmp_cwd, nircam_rate, subtract):
+
+    # give them all different suffixes to ensure they get stripped properly
+    im1 = nircam_rate.copy()
+    im1.meta.filename = "one_tweakregstep.fits"
+    im2 = im1.copy()
+    im2.meta.filename = "two_unknown.fits"
+    im3 = im1.copy()
+    im3.meta.filename = "dir/three.fits"
+
+    # give filenames in skyfile same stems but different suffix
+    fnames_skyfile = ["other_dir/one_cal.fits", "two_unknown_cal.fits", "three_cal.fits"]
+
+    container = [im1, im2, im3]
+
+    # put levels into the skylist file for when skylist='user'
+    levels = [9.12, 8.28, 2.56]
+    fnames = [model.meta.filename for model in container]
+
+    for im, lev in zip(container, levels):
+        im.data += lev
+
+    skyfile = "skylist.txt"
+    with open(skyfile, "w") as f:
+        for fname, lev in zip(fnames_skyfile, levels):
+            f.write(f"{fname} {lev}\n")
+
+    #test good inputs
+    result = SkyMatchStep.call(
+        container,
+        subtract=subtract,
+        skymethod="user",
+        skylist=skyfile,
+    )
+
+    ref_levels = levels
+    sub_levels = np.subtract(levels, ref_levels)
+
+    with result:
+        for im, lev, rlev, slev in zip(result, levels, ref_levels, sub_levels):
+            # check that meta was set correctly:
+            assert im.meta.background.method == "user"
+            assert im.meta.background.subtracted == subtract
+
+            # test computed/measured sky values:
+            assert abs(im.meta.background.level - rlev) < 0.01
+
+            # test
+            if subtract:
+                assert abs(np.mean(im.data) - slev) < 0.01
+            else:
+                assert abs(np.mean(im.data) - lev) < 0.01
+            result.shelve(im, modify=False)
+
+
+    # test failures
+    # no skylist file
+    with pytest.raises(ValueError):
+        # skylist must be provided
+        SkyMatchStep.call(
+            container,
+            skymethod='user',
+        )
+
+    # skylist file doesn't have right number of lines
+    skyfile = "skylist_short.txt"
+    with open(skyfile, "w") as f:
+        for fname, lev in zip(fnames[1:], levels[1:]):
+            f.write(f"{fname} {lev}\n")
+
+    with pytest.raises(ValueError):
+        SkyMatchStep.call(
+            container,
+            skymethod='user',
+            skylist=skyfile
+        )
+
+    # skylist file does not contain all filenames
+    skyfile = "skylist_missing.txt"
+    fnames_wrong = ["two.fits"] + fnames[1:]
+    with open(skyfile, "w") as f:
+        for fname, lev in zip(fnames_wrong, levels):
+            f.write(f"{fname} {lev}\n")
+
+    with pytest.raises(ValueError):
+        SkyMatchStep.call(
+            container,
+            skymethod='user',
+            skylist=skyfile
+        )