From 55fc57a4a30b20e85ca40c947006f53efcac4b23 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Fri, 6 Dec 2024 11:54:03 -0500
Subject: [PATCH 01/21] Rename ResampleData class weight type parameter.

---
 romancal/resample/resample.py            | 4 ++--
 romancal/resample/tests/test_resample.py | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py
index 0e9ad543c..4def81aff 100644
--- a/romancal/resample/resample.py
+++ b/romancal/resample/resample.py
@@ -52,7 +52,7 @@ def __init__(
         pixfrac=1.0,
         kernel="square",
         fillval="INDEF",
-        wht_type="ivm",
+        weight_type="ivm",
         good_bits="0",
         pscale_ratio=1.0,
         pscale=None,
@@ -94,7 +94,7 @@ def __init__(
         self.pixfrac = pixfrac
         self.kernel = kernel
         self.fillval = fillval
-        self.weight_type = wht_type
+        self.weight_type = weight_type
         self.good_bits = good_bits
         self.in_memory = kwargs.get("in_memory", True)
         if "target" in input_models.asn:
diff --git a/romancal/resample/tests/test_resample.py b/romancal/resample/tests/test_resample.py
index 2923a2625..2e65e8dcc 100644
--- a/romancal/resample/tests/test_resample.py
+++ b/romancal/resample/tests/test_resample.py
@@ -309,7 +309,7 @@ def test_resampledata_init(exposure_1):
         pixfrac=pixfrac,
         kernel=kernel,
         fillval=fillval,
-        wht_type=wht_type,
+        weight_type=wht_type,
         good_bits=good_bits,
         pscale_ratio=pscale_ratio,
         pscale=pscale,
@@ -710,7 +710,7 @@ def test_resampledata_do_drizzle_default_single_exposure_weight_array(
     """Test that resample methods return non-empty weight arrays."""
 
     input_models = ModelLibrary(exposure_1)
-    resample_data = ResampleData(input_models, wht_type=weight_type)
+    resample_data = ResampleData(input_models, weight_type=weight_type)
 
     output_models_many_to_one = resample_data.resample_many_to_one()
     output_models_many_to_many = resample_data.resample_many_to_many()

From 11595621e2e3c7d01661e9df811dc2167f33d889 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Mon, 9 Dec 2024 15:29:32 -0500
Subject: [PATCH 02/21] Add var_sky array to mosaic and enable usage in drizzle
 weight image.

---
 romancal/regtest/test_resample.py        |  6 +++--
 romancal/resample/resample.py            | 11 ++++-----
 romancal/resample/resample_utils.py      | 31 ++++++++++++++++++++++++
 romancal/resample/tests/test_resample.py |  2 +-
 4 files changed, 41 insertions(+), 9 deletions(-)

diff --git a/romancal/regtest/test_resample.py b/romancal/regtest/test_resample.py
index a1c715925..446102e51 100644
--- a/romancal/regtest/test_resample.py
+++ b/romancal/regtest/test_resample.py
@@ -60,6 +60,7 @@ def test_resample_single_file(rtdata, ignore_asdf_paths):
                         "var_poisson",
                         "var_rnoise",
                         "var_flat",
+                        "var_sky",
                     ]
                 )
             }"""
@@ -76,6 +77,7 @@ def test_resample_single_file(rtdata, ignore_asdf_paths):
                     np.sum(~np.isnan(getattr(resample_out, x))) for x in [
                         "var_poisson",
                         "var_rnoise",
+                        "var_sky",
                     ]
                 )
             }"""
@@ -94,14 +96,14 @@ def test_resample_single_file(rtdata, ignore_asdf_paths):
                             np.isnan(getattr(resample_out, x)),
                             np.equal(getattr(resample_out, x), 0)
                         )
-                    ) > 0 for x in ["var_poisson", "var_rnoise", "var_flat"]
+                    ) > 0 for x in ["var_poisson", "var_rnoise", "var_flat", "var_sky"]
                 )
 
             }"""
     )
     assert all(
         np.sum(np.isnan(getattr(resample_out, x)))
-        for x in ["var_poisson", "var_rnoise", "var_flat"]
+        for x in ["var_poisson", "var_rnoise", "var_flat", "var_sky"]
     )
 
     step.log.info(
diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py
index 4def81aff..c01192491 100644
--- a/romancal/resample/resample.py
+++ b/romancal/resample/resample.py
@@ -87,6 +87,8 @@ def __init__(
             )
 
         self.input_models = input_models
+        # add sky variance array
+        resample_utils.add_var_sky_array(self.input_models)
         self.output_filename = output
         self.pscale_ratio = pscale_ratio
         self.single = single
@@ -157,9 +159,6 @@ def __init__(
         #         f'Model cannot be instantiated.'
         #     )
 
-        # NOTE: wait for William to fix bug in datamodels' init and then
-        # use datamodels.ImageModel(shape=(nx, ny)) instead of mk_datamodel()
-
         # n_images sets the number of context image planes.
         # This should be 1 to start (not the default of 2).
         self.blank_output = maker_utils.mk_datamodel(
@@ -197,6 +196,8 @@ def __init__(
             self.blank_output["individual_image_cal_logs"] = [
                 model.meta.cal_logs for model in models
             ]
+            # add sky variance array
+            self.blank_output["var_sky"] = np.zeros_like(self.blank_output.var_flat)
             for i, m in enumerate(models):
                 self.input_models.shelve(m, i, modify=False)
 
@@ -396,11 +397,11 @@ def resample_many_to_one(self):
         self.resample_variance_array("var_rnoise", output_model)
         self.resample_variance_array("var_poisson", output_model)
         self.resample_variance_array("var_flat", output_model)
+        self.resample_variance_array("var_sky", output_model)
 
         # Make exposure time image
         exptime_tot = self.resample_exposure_time(output_model)
 
-        # TODO: fix unit here
         output_model.err = np.sqrt(
             np.nansum(
                 [
@@ -414,7 +415,6 @@ def resample_many_to_one(self):
 
         self.update_exposure_times(output_model, exptime_tot)
 
-        # TODO: fix RAD to expect a context image datatype of int32
         output_model.context = output_model.context.astype(np.uint32)
 
         return ModelLibrary([output_model])
@@ -493,7 +493,6 @@ def resample_variance_array(self, name, output_model):
 
         # We now have a sum of the inverse resampled variances.  We need the
         # inverse of that to get back to units of variance.
-        # TODO: fix unit here
         output_variance = np.reciprocal(inverse_variance_sum)
 
         setattr(output_model, name, output_variance)
diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index 4b8c78aa0..e8c6ce60f 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -10,6 +10,7 @@
 from stcal.alignment.util import wcs_from_footprints
 
 from romancal.assign_wcs.utils import wcs_bbox_from_shape
+from romancal.datamodels.library import ModelLibrary
 
 log = logging.getLogger(__name__)
 log.setLevel(logging.DEBUG)
@@ -163,6 +164,20 @@ def build_driz_weight(
     elif weight_type == "exptime":
         exptime = model.meta.exposure.exposure_time
         result = exptime * dqmask
+    elif weight_type == "ivsky":
+        if (
+            hasattr(model, "var_sky")
+            and model.var_sky is not None
+            and model.var_sky.shape == model.data.shape
+        ):
+            with np.errstate(divide="ignore", invalid="ignore"):
+                inv_sky_variance = model.var_sky**-1
+            inv_sky_variance[~np.isfinite(inv_sky_variance)] = 1
+        else:
+            warnings.warn(
+                "var_rnoise and var_poisson arrays are not available. Setting drizzle weight map to 1"
+            )
+        result = inv_sky_variance * dqmask
     elif weight_type is None:
         result = np.ones(model.data.shape, dtype=model.data.dtype) * dqmask
     else:
@@ -402,3 +417,19 @@ def resample_range(data_shape, bbox=None):
         ymax = min(data_shape[0] - 1, int(y2 + 0.5))
 
     return xmin, xmax, ymin, ymax
+
+
+def add_var_sky_array(input_models: ModelLibrary):
+    input_models = (
+        input_models
+        if isinstance(input_models, ModelLibrary)
+        else ModelLibrary([input_models])
+    )
+    with input_models:
+        ref_img = input_models.borrow(index=0)
+        var_sky = np.zeros_like(ref_img.data)
+        input_models.shelve(model=ref_img, index=0)
+        for i, img in enumerate(input_models):
+            var_sky += img.var_rnoise + img.var_poisson / img.data * np.median(img.data)
+            img["var_sky"] = var_sky
+            input_models.shelve(img, i, modify=True)
diff --git a/romancal/resample/tests/test_resample.py b/romancal/resample/tests/test_resample.py
index 2e65e8dcc..044f280b5 100644
--- a/romancal/resample/tests/test_resample.py
+++ b/romancal/resample/tests/test_resample.py
@@ -702,7 +702,7 @@ def get_footprint(model, index):
     )
 
 
-@pytest.mark.parametrize("weight_type", ["ivm", "exptime"])
+@pytest.mark.parametrize("weight_type", ["ivm", "exptime", "ivsky"])
 def test_resampledata_do_drizzle_default_single_exposure_weight_array(
     exposure_1,
     weight_type,

From 6d9f4373aad5c6f9f4f2d866581adc5211b1ab22 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Wed, 11 Dec 2024 09:20:01 -0500
Subject: [PATCH 03/21] Add changelog entry.

---
 changes/1556.resample.rst | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 changes/1556.resample.rst

diff --git a/changes/1556.resample.rst b/changes/1556.resample.rst
new file mode 100644
index 000000000..5cce86d3c
--- /dev/null
+++ b/changes/1556.resample.rst
@@ -0,0 +1 @@
+Allow resample to use inverse sky variance when building the drizzle weight map.

From d7df70968da7e10e75039e3f1eb5cd5b22eed03b Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Wed, 11 Dec 2024 09:26:49 -0500
Subject: [PATCH 04/21] Update docs and docstring.

---
 docs/roman/resample/arguments.rst   | 4 ++++
 romancal/resample/resample_utils.py | 2 +-
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/docs/roman/resample/arguments.rst b/docs/roman/resample/arguments.rst
index b72bb630a..d07d70bdc 100644
--- a/docs/roman/resample/arguments.rst
+++ b/docs/roman/resample/arguments.rst
@@ -81,6 +81,10 @@ image.
     not exist, the variance is set to 1 for all pixels (equal weighting).
     If `weight_type=exptime`, the scaling value will be set equal to the
     exposure time found in the image header.
+    If `weight_type=ivsky`, the scaling value will be determined per-pixel 
+    using the inverse of the sky variance (VAR_SKY) array calculated in the 
+    resample step for each input image. If the VAR_SKY array does
+    not exist, the variance is set to 1 for all pixels (equal weighting).
 
 ``--single`` (bool, default=False)
     If set to `True`, resample each input image into a separate output.  If
diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index e8c6ce60f..73df1d397 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -112,7 +112,7 @@ def build_driz_weight(
     model : object
         The input model.
     weight_type : str, optional
-        The type of weight to use. Allowed values are 'ivm' or 'exptime'.
+        The type of weight to use. Allowed values are 'ivm', 'exptime', or 'ivsky'.
         Defaults to None.
     good_bits : str, optional
         The good bits to use for building the mask. Defaults to None.

From a3e8634a2f6229dde90d3c5e24a220035e9696a7 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 11 Dec 2024 14:27:10 +0000
Subject: [PATCH 05/21] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 docs/roman/resample/arguments.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/roman/resample/arguments.rst b/docs/roman/resample/arguments.rst
index d07d70bdc..c0078739b 100644
--- a/docs/roman/resample/arguments.rst
+++ b/docs/roman/resample/arguments.rst
@@ -81,8 +81,8 @@ image.
     not exist, the variance is set to 1 for all pixels (equal weighting).
     If `weight_type=exptime`, the scaling value will be set equal to the
     exposure time found in the image header.
-    If `weight_type=ivsky`, the scaling value will be determined per-pixel 
-    using the inverse of the sky variance (VAR_SKY) array calculated in the 
+    If `weight_type=ivsky`, the scaling value will be determined per-pixel
+    using the inverse of the sky variance (VAR_SKY) array calculated in the
     resample step for each input image. If the VAR_SKY array does
     not exist, the variance is set to 1 for all pixels (equal weighting).
 

From 25e4a16aabb9f780dc3426853aec8375fd561621 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Wed, 11 Dec 2024 16:47:53 -0500
Subject: [PATCH 06/21] Remove unnecessary code.

---
 romancal/resample/resample_utils.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index 73df1d397..964f5a47c 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -420,6 +420,18 @@ def resample_range(data_shape, bbox=None):
 
 
 def add_var_sky_array(input_models: ModelLibrary):
+    """
+    Add sky variance array to each input model.
+
+    Parameters
+    ----------
+    input_models : ModelLibrary
+        A library of models or a single model to which the sky variance array will be added.
+
+    Returns
+    -------
+    None
+    """
     input_models = (
         input_models
         if isinstance(input_models, ModelLibrary)
@@ -427,9 +439,9 @@ def add_var_sky_array(input_models: ModelLibrary):
     )
     with input_models:
         ref_img = input_models.borrow(index=0)
-        var_sky = np.zeros_like(ref_img.data)
         input_models.shelve(model=ref_img, index=0)
         for i, img in enumerate(input_models):
-            var_sky += img.var_rnoise + img.var_poisson / img.data * np.median(img.data)
-            img["var_sky"] = var_sky
+            img["var_sky"] = img.var_rnoise + img.var_poisson / img.data * np.median(
+                img.data
+            )
             input_models.shelve(img, i, modify=True)

From 9f6df6f7b98446b799d8889cb8545a47ad28c61c Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Thu, 12 Dec 2024 10:16:17 -0500
Subject: [PATCH 07/21] Fix terminology in resample docs.

---
 docs/roman/resample/arguments.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/roman/resample/arguments.rst b/docs/roman/resample/arguments.rst
index c0078739b..78bcf3b49 100644
--- a/docs/roman/resample/arguments.rst
+++ b/docs/roman/resample/arguments.rst
@@ -78,13 +78,13 @@ image.
     If `weight_type=ivm` (the default), the scaling value
     will be determined per-pixel using the inverse of the read noise
     (VAR_RNOISE) array stored in each input image. If the VAR_RNOISE array does
-    not exist, the variance is set to 1 for all pixels (equal weighting).
+    not exist, the weight is set to 1 for all pixels (equal weighting).
     If `weight_type=exptime`, the scaling value will be set equal to the
     exposure time found in the image header.
     If `weight_type=ivsky`, the scaling value will be determined per-pixel
     using the inverse of the sky variance (VAR_SKY) array calculated in the
     resample step for each input image. If the VAR_SKY array does
-    not exist, the variance is set to 1 for all pixels (equal weighting).
+    not exist, the weight is set to 1 for all pixels (equal weighting).
 
 ``--single`` (bool, default=False)
     If set to `True`, resample each input image into a separate output.  If

From d277a66086ccd25a9984c993b7a06b8643a41427 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Mon, 16 Dec 2024 15:33:19 -0500
Subject: [PATCH 08/21] Update step spec.

---
 romancal/resample/resample_step.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/romancal/resample/resample_step.py b/romancal/resample/resample_step.py
index 8796a52ff..114979538 100644
--- a/romancal/resample/resample_step.py
+++ b/romancal/resample/resample_step.py
@@ -58,7 +58,7 @@ class ResampleStep(RomanStep):
         pixfrac = float(default=1.0) # change back to None when drizpar reference files are updated
         kernel = string(default='square') # change back to None when drizpar reference files are updated
         fillval = string(default='INDEF' ) # change back to None when drizpar reference files are updated
-        weight_type = option('ivm', 'exptime', None, default='ivm')  # change back to None when drizpar ref update
+        weight_type = option('ivm', 'exptime', 'ivsky', None, default='ivm')  # change back to None when drizpar ref update
         output_shape = int_list(min=2, max=2, default=None)  # [x, y] order
         crpix = float_list(min=2, max=2, default=None)
         crval = float_list(min=2, max=2, default=None)

From 3ba72ff3a456b64fd8e49ef24ec56ff29d24806a Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Mon, 16 Dec 2024 15:39:51 -0500
Subject: [PATCH 09/21] Add default value for inverse of variance sky if
 var_sky array is not present.

---
 romancal/resample/resample_utils.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index 964f5a47c..50efb90d8 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -177,6 +177,7 @@ def build_driz_weight(
             warnings.warn(
                 "var_rnoise and var_poisson arrays are not available. Setting drizzle weight map to 1"
             )
+            inv_sky_variance = 1.0
         result = inv_sky_variance * dqmask
     elif weight_type is None:
         result = np.ones(model.data.shape, dtype=model.data.dtype) * dqmask

From 5900cf05adeacc90b5bf66cb42efe4092d45420a Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Mon, 16 Dec 2024 15:45:19 -0500
Subject: [PATCH 10/21] Remove unnecessary check.

---
 romancal/resample/resample_utils.py | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index 50efb90d8..e1e2b8606 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -422,22 +422,17 @@ def resample_range(data_shape, bbox=None):
 
 def add_var_sky_array(input_models: ModelLibrary):
     """
-    Add sky variance array to each input model.
+    Add sky variance array to each model of a ModelLibrary.
 
     Parameters
     ----------
     input_models : ModelLibrary
-        A library of models or a single model to which the sky variance array will be added.
+        A library of models to which the sky variance array will be added.
 
     Returns
     -------
     None
     """
-    input_models = (
-        input_models
-        if isinstance(input_models, ModelLibrary)
-        else ModelLibrary([input_models])
-    )
     with input_models:
         ref_img = input_models.borrow(index=0)
         input_models.shelve(model=ref_img, index=0)

From bc62c2b0328e7cb21630281e378d4392810794e5 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Mon, 16 Dec 2024 21:56:39 -0500
Subject: [PATCH 11/21] Prevent zero division issues.

---
 romancal/resample/resample_utils.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index e1e2b8606..34a682ffa 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -437,7 +437,10 @@ def add_var_sky_array(input_models: ModelLibrary):
         ref_img = input_models.borrow(index=0)
         input_models.shelve(model=ref_img, index=0)
         for i, img in enumerate(input_models):
-            img["var_sky"] = img.var_rnoise + img.var_poisson / img.data * np.median(
-                img.data
-            )
+            if np.all(img.data != 0):
+                img["var_sky"] = (
+                    img.var_rnoise + img.var_poisson / img.data * np.median(img.data)
+                )
+            else:
+                raise ValueError("Input model contains invalid data array.")
             input_models.shelve(img, i, modify=True)

From 5717c3bf454f3ca2131ee9ec3dcbd5af69d812cd Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Mon, 30 Dec 2024 15:04:30 -0500
Subject: [PATCH 12/21] Add description of calculation of VAR_SKY.

---
 docs/roman/resample/arguments.rst | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/docs/roman/resample/arguments.rst b/docs/roman/resample/arguments.rst
index 78bcf3b49..668ef208d 100644
--- a/docs/roman/resample/arguments.rst
+++ b/docs/roman/resample/arguments.rst
@@ -77,14 +77,20 @@ image.
     The weighting type for each input image.
     If `weight_type=ivm` (the default), the scaling value
     will be determined per-pixel using the inverse of the read noise
-    (VAR_RNOISE) array stored in each input image. If the VAR_RNOISE array does
+    (``VAR_RNOISE``) array stored in each input image. If the ``VAR_RNOISE`` array does
     not exist, the weight is set to 1 for all pixels (equal weighting).
     If `weight_type=exptime`, the scaling value will be set equal to the
     exposure time found in the image header.
     If `weight_type=ivsky`, the scaling value will be determined per-pixel
-    using the inverse of the sky variance (VAR_SKY) array calculated in the
-    resample step for each input image. If the VAR_SKY array does
-    not exist, the weight is set to 1 for all pixels (equal weighting).
+    using the inverse of the sky variance (``VAR_SKY``) array calculated in the
+    resample step for each input image. ``VAR_SKY`` is given by the following equation:
+
+    .. math::
+
+        \text{VAR\_SKY} = \text{VAR\_RNOISE} + \text{VAR\_POISSON} \, \frac{ med(\text{DATA}) }{ \text{DATA} },
+
+    where :math:`\text{DATA}` and :math:`med(\text{DATA})` correspond to the data array and its median, respectively.
+    If the ``VAR_SKY`` array does not exist, the weight is set to 1 for all pixels (equal weighting).
 
 ``--single`` (bool, default=False)
     If set to `True`, resample each input image into a separate output.  If

From e8925c61de470b1e15ea22ede962d2c43b793930 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Thu, 9 Jan 2025 20:22:38 -0500
Subject: [PATCH 13/21] Fix creating base image.

---
 .../outlier_detection/tests/test_outlier_detection.py    | 9 +++++++--
 romancal/resample/tests/test_resample.py                 | 1 +
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/romancal/outlier_detection/tests/test_outlier_detection.py b/romancal/outlier_detection/tests/test_outlier_detection.py
index 7779dadc9..cf7b1ae64 100644
--- a/romancal/outlier_detection/tests/test_outlier_detection.py
+++ b/romancal/outlier_detection/tests/test_outlier_detection.py
@@ -103,9 +103,11 @@ def test_outlier_do_detection_write_files_to_custom_location(tmp_path, base_imag
     img_1 = base_image()
     img_1.meta.filename = "img1_cal.asdf"
     img_1.meta.background.level = 0
+    img_1.data[:] = 0.01
     img_2 = base_image()
     img_2.meta.filename = "img2_cal.asdf"
     img_2.meta.background.level = 0
+    img_2.data[:] = 0.01
     input_models = ModelLibrary([img_1, img_2])
 
     outlier_step = OutlierDetectionStep(
@@ -135,11 +137,13 @@ def test_find_outliers(tmp_path, base_image, on_disk):
     cr_value = 100
     source_value = 10
     err_value = 10  # snr=1
+    sky_value = 0.5
 
     imgs = []
     for i in range(3):
         img = base_image()
-        img.data[42, 72] = source_value
+        img.data[:] = sky_value
+        img.data[42, 72] += source_value
         img.err[:] = err_value
         img.meta.filename = str(tmp_path / f"img{i}_suffix.asdf")
         img.meta.observation.observation_id = str(i)
@@ -207,7 +211,8 @@ def test_identical_images(tmp_path, base_image, caplog):
     img_1_input_coords = np.array(
         [(5, 45), (25, 25), (45, 85), (65, 65), (85, 5)], dtype=[("x", int), ("y", int)]
     )
-    img_1.data[img_1_input_coords["x"], img_1_input_coords["y"]] = 100000
+    img_1.data[:] = 0.01
+    img_1.data[img_1_input_coords["x"], img_1_input_coords["y"]] += 100
 
     img_2 = img_1.copy()
     img_2.meta.filename = "img2_suffix.asdf"
diff --git a/romancal/resample/tests/test_resample.py b/romancal/resample/tests/test_resample.py
index 044f280b5..98bd9bc86 100644
--- a/romancal/resample/tests/test_resample.py
+++ b/romancal/resample/tests/test_resample.py
@@ -100,6 +100,7 @@ def create_image(self):
             },
         )
         # data from WFISim simulation of SCA #01
+        l2.data[:] = 0.01
         l2.meta.filename = self.filename
         l2.meta["wcs"] = create_wcs_object_without_distortion(
             fiducial_world=self.fiducial_world,

From 7173bf3438bd16ce90ac11b1770c8a013df52828 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Fri, 10 Jan 2025 17:16:36 -0500
Subject: [PATCH 14/21] Fix missing stacklevel keyword in call to
 warnings.warn().

---
 romancal/resample/resample_utils.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index 34a682ffa..199aec7d4 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -175,7 +175,8 @@ def build_driz_weight(
             inv_sky_variance[~np.isfinite(inv_sky_variance)] = 1
         else:
             warnings.warn(
-                "var_rnoise and var_poisson arrays are not available. Setting drizzle weight map to 1"
+                "var_rnoise and var_poisson arrays are not available. Setting drizzle weight map to 1",
+                stacklevel=2,
             )
             inv_sky_variance = 1.0
         result = inv_sky_variance * dqmask

From f850013829609c38fc2d864c7ded9ea0d5e67c77 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Thu, 16 Jan 2025 12:00:30 -0500
Subject: [PATCH 15/21] Fix unit test by adding non-zero background level.

---
 romancal/outlier_detection/tests/test_outlier_detection.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/romancal/outlier_detection/tests/test_outlier_detection.py b/romancal/outlier_detection/tests/test_outlier_detection.py
index cf7b1ae64..4233b8c4d 100644
--- a/romancal/outlier_detection/tests/test_outlier_detection.py
+++ b/romancal/outlier_detection/tests/test_outlier_detection.py
@@ -204,14 +204,15 @@ def test_identical_images(tmp_path, base_image, caplog):
     """
     Test that OutlierDetection does not flag any outliers in the DQ array if images are identical.
     """
+    background_level = 0.01
     img_1 = base_image()
     img_1.meta.filename = "img1_suffix.asdf"
-    img_1.meta.background.level = 0
+    img_1.meta.background.level = background_level
     # add outliers
     img_1_input_coords = np.array(
         [(5, 45), (25, 25), (45, 85), (65, 65), (85, 5)], dtype=[("x", int), ("y", int)]
     )
-    img_1.data[:] = 0.01
+    img_1.data[:] = background_level
     img_1.data[img_1_input_coords["x"], img_1_input_coords["y"]] += 100
 
     img_2 = img_1.copy()

From fadc18d1d175b856d2927928218b9e13afeaa928 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Thu, 16 Jan 2025 13:34:54 -0500
Subject: [PATCH 16/21] Add var_sky to model error calculation.

---
 romancal/resample/resample.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py
index c01192491..f814c57f6 100644
--- a/romancal/resample/resample.py
+++ b/romancal/resample/resample.py
@@ -408,6 +408,7 @@ def resample_many_to_one(self):
                     output_model.var_rnoise,
                     output_model.var_poisson,
                     output_model.var_flat,
+                    output_model.var_sky,
                 ],
                 axis=0,
             )

From c448778299f53c23c5482fb6ea4815f28e747b6a Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Wed, 29 Jan 2025 15:50:44 -0500
Subject: [PATCH 17/21] Address review comments.

---
 romancal/resample/resample.py            |  1 -
 romancal/resample/resample_utils.py      | 15 +++++++++------
 romancal/resample/tests/test_resample.py |  6 ++++++
 3 files changed, 15 insertions(+), 7 deletions(-)

diff --git a/romancal/resample/resample.py b/romancal/resample/resample.py
index 7f6f91a1f..376fb38bc 100644
--- a/romancal/resample/resample.py
+++ b/romancal/resample/resample.py
@@ -389,7 +389,6 @@ def resample_many_to_one(self):
                     output_model.var_rnoise,
                     output_model.var_poisson,
                     output_model.var_flat,
-                    output_model.var_sky,
                 ],
                 axis=0,
             )
diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index 199aec7d4..c47de1bdd 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -172,7 +172,7 @@ def build_driz_weight(
         ):
             with np.errstate(divide="ignore", invalid="ignore"):
                 inv_sky_variance = model.var_sky**-1
-            inv_sky_variance[~np.isfinite(inv_sky_variance)] = 1
+            inv_sky_variance[~np.isfinite(inv_sky_variance)] = 0
         else:
             warnings.warn(
                 "var_rnoise and var_poisson arrays are not available. Setting drizzle weight map to 1",
@@ -438,10 +438,13 @@ def add_var_sky_array(input_models: ModelLibrary):
         ref_img = input_models.borrow(index=0)
         input_models.shelve(model=ref_img, index=0)
         for i, img in enumerate(input_models):
-            if np.all(img.data != 0):
-                img["var_sky"] = (
-                    img.var_rnoise + img.var_poisson / img.data * np.median(img.data)
-                )
-            else:
+            try:
+                ok_data = img.data != 0
+                img["var_sky"] = np.empty_like(img.data)
+                img["var_sky"][ok_data] = img.var_rnoise[ok_data] + img.var_poisson[
+                    ok_data
+                ] / img.data[ok_data] * np.median(img.data)
+                img["var_sky"][~ok_data] = img.var_rnoise[~ok_data]
+            except:
                 raise ValueError("Input model contains invalid data array.")
             input_models.shelve(img, i, modify=True)
diff --git a/romancal/resample/tests/test_resample.py b/romancal/resample/tests/test_resample.py
index 3fa4d81db..12ffcf36b 100644
--- a/romancal/resample/tests/test_resample.py
+++ b/romancal/resample/tests/test_resample.py
@@ -702,6 +702,10 @@ def test_resampledata_do_drizzle_default_single_exposure_weight_array(
 ):
     """Test that resample methods return non-empty weight arrays."""
 
+    # adding a few zero flux pixels
+    for i, model in enumerate(exposure_1):
+        model.data[10 + i, 40 - i] = 0
+
     input_models = ModelLibrary(exposure_1)
     resample_data = ResampleData(input_models, weight_type=weight_type)
 
@@ -713,6 +717,8 @@ def test_resampledata_do_drizzle_default_single_exposure_weight_array(
         many_to_one_model = output_models_many_to_one.borrow(0)
         assert np.any(many_to_one_model.weight > 0)
         assert np.any(many_to_many_model.weight > 0)
+        assert many_to_one_model.data[10, 40] == 0
+        assert many_to_many_model.data[10, 40] == 0
         output_models_many_to_many.shelve(many_to_many_model, 0, modify=False)
         output_models_many_to_one.shelve(many_to_one_model, 0, modify=False)
 

From 7c3a31c93b46dea6d6b98deb3e06b0428d3e64cf Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Wed, 29 Jan 2025 16:44:44 -0500
Subject: [PATCH 18/21] Fix check style issue.

---
 romancal/resample/resample_utils.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index c47de1bdd..816e1cab5 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -445,6 +445,6 @@ def add_var_sky_array(input_models: ModelLibrary):
                     ok_data
                 ] / img.data[ok_data] * np.median(img.data)
                 img["var_sky"][~ok_data] = img.var_rnoise[~ok_data]
-            except:
-                raise ValueError("Input model contains invalid data array.")
+            except (AttributeError, KeyError, TypeError, ValueError) as e:
+                raise ValueError("Input model contains invalid data array.") from e
             input_models.shelve(img, i, modify=True)

From fcb2c17ab8a9dc3733d019aa5e4850fd30c131d5 Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Wed, 29 Jan 2025 17:29:27 -0500
Subject: [PATCH 19/21] Improve docs and logs.

---
 docs/roman/resample/arguments.rst   | 3 ++-
 romancal/resample/resample_utils.py | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/docs/roman/resample/arguments.rst b/docs/roman/resample/arguments.rst
index 3d00a0c01..8ffb4082f 100644
--- a/docs/roman/resample/arguments.rst
+++ b/docs/roman/resample/arguments.rst
@@ -90,7 +90,8 @@ image.
         \text{VAR\_SKY} = \text{VAR\_RNOISE} + \text{VAR\_POISSON} \, \frac{ med(\text{DATA}) }{ \text{DATA} },
 
     where :math:`\text{DATA}` and :math:`med(\text{DATA})` correspond to the data array and its median, respectively.
-    If the ``VAR_SKY`` array does not exist, the weight is set to 1 for all pixels (equal weighting).
+    If the ``VAR_SKY`` array does not exist (which implies missing ``VAR_RNOISE`` and/or ``VAR_POISSON``), 
+    the weight is set to 1 for all pixels (equal weighting).
 
 ``--in_memory`` (bool, default=True)
     If set to `False`, write output datamodel to disk.
diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py
index 816e1cab5..461e3aa0a 100644
--- a/romancal/resample/resample_utils.py
+++ b/romancal/resample/resample_utils.py
@@ -175,7 +175,7 @@ def build_driz_weight(
             inv_sky_variance[~np.isfinite(inv_sky_variance)] = 0
         else:
             warnings.warn(
-                "var_rnoise and var_poisson arrays are not available. Setting drizzle weight map to 1",
+                "var_sky array is not available. Setting drizzle weight map to 1",
                 stacklevel=2,
             )
             inv_sky_variance = 1.0

From 6507bfd89de9e4f85a1abc5bb3c34ad59f8f53b4 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 29 Jan 2025 22:30:13 +0000
Subject: [PATCH 20/21] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 docs/roman/resample/arguments.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/roman/resample/arguments.rst b/docs/roman/resample/arguments.rst
index 8ffb4082f..ce735251a 100644
--- a/docs/roman/resample/arguments.rst
+++ b/docs/roman/resample/arguments.rst
@@ -90,7 +90,7 @@ image.
         \text{VAR\_SKY} = \text{VAR\_RNOISE} + \text{VAR\_POISSON} \, \frac{ med(\text{DATA}) }{ \text{DATA} },
 
     where :math:`\text{DATA}` and :math:`med(\text{DATA})` correspond to the data array and its median, respectively.
-    If the ``VAR_SKY`` array does not exist (which implies missing ``VAR_RNOISE`` and/or ``VAR_POISSON``), 
+    If the ``VAR_SKY`` array does not exist (which implies missing ``VAR_RNOISE`` and/or ``VAR_POISSON``),
     the weight is set to 1 for all pixels (equal weighting).
 
 ``--in_memory`` (bool, default=True)

From 332943cbb08e8fea7b50dd16fb698bb8dbe72dff Mon Sep 17 00:00:00 2001
From: "M. Teodoro" <mteodoro@stsci.edu>
Date: Thu, 30 Jan 2025 10:09:47 -0500
Subject: [PATCH 21/21] Remove escape characters.

---
 docs/roman/resample/arguments.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/roman/resample/arguments.rst b/docs/roman/resample/arguments.rst
index ce735251a..ff4beb681 100644
--- a/docs/roman/resample/arguments.rst
+++ b/docs/roman/resample/arguments.rst
@@ -87,7 +87,7 @@ image.
 
     .. math::
 
-        \text{VAR\_SKY} = \text{VAR\_RNOISE} + \text{VAR\_POISSON} \, \frac{ med(\text{DATA}) }{ \text{DATA} },
+        \text{VAR_SKY} = \text{VAR_RNOISE} + \text{VAR_POISSON} \, \frac{ med(\text{DATA}) }{ \text{DATA} },
 
     where :math:`\text{DATA}` and :math:`med(\text{DATA})` correspond to the data array and its median, respectively.
     If the ``VAR_SKY`` array does not exist (which implies missing ``VAR_RNOISE`` and/or ``VAR_POISSON``),