Hackday mteodoro sky variance readnoise

changes/1556.resample.rst

@@ -0,0 +1 @@
+Allow resample to use inverse sky variance when building the drizzle weight map.

docs/roman/resample/arguments.rst

@@ -77,10 +77,21 @@ 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
-    not exist, the variance is set to 1 for all pixels (equal weighting).
+    (``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. ``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 (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.

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)
@@ -200,14 +204,16 @@ 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[img_1_input_coords["x"], img_1_input_coords["y"]] = 100000
+    img_1.data[:] = background_level
+    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"

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(

romancal/resample/resample.py

@@ -51,7 +51,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,
@@ -86,13 +86,15 @@ 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
         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:
@@ -178,6 +180,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)
 
@@ -374,11 +378,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(
                 [
@@ -392,7 +396,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])
@@ -471,7 +474,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)

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)

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)
@@ -111,7 +112,7 @@
     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.
@@ -163,6 +164,22 @@
     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)] = 0
+        else:
+            warnings.warn(
+                "var_sky array is not available. Setting drizzle weight map to 1",
+                stacklevel=2,
+            )
+            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
     else:
@@ -402,3 +419,32 @@
         ymax = min(data_shape[0] - 1, int(y2 + 0.5))
 
     return xmin, xmax, ymin, ymax
+
+
+def add_var_sky_array(input_models: ModelLibrary):
+    """
+    Add sky variance array to each model of a ModelLibrary.
+
+    Parameters
+    ----------
+    input_models : ModelLibrary
+        A library of models to which the sky variance array will be added.
+
+    Returns
+    -------
+    None
+    """
+    with input_models:
+        ref_img = input_models.borrow(index=0)
+        input_models.shelve(model=ref_img, index=0)
+        for i, img in enumerate(input_models):
+            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 (AttributeError, KeyError, TypeError, ValueError) as e:
+                raise ValueError("Input model contains invalid data array.") from e
+            input_models.shelve(img, i, modify=True)

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,
@@ -305,7 +306,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,
@@ -694,15 +695,19 @@ 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,
 ):
     """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, 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()
@@ -712,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)