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,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
-    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, 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

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

@@ -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,
@@ -87,14 +87,16 @@ 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.blendheaders = blendheaders
         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:
@@ -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,25 +397,25 @@ 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(
                 [
                     output_model.var_rnoise,
                     output_model.var_poisson,
                     output_model.var_flat,
+                    output_model.var_sky,
                 ],
                 axis=0,
             )
         )
 
         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 +494,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)] = 1
+        else:
+            warnings.warn(
+                "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
     elif weight_type is None:
         result = np.ones(model.data.shape, dtype=model.data.dtype) * dqmask
     else:
@@ -402,3 +419,29 @@
         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):
+            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)

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,
@@ -302,7 +303,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,
@@ -679,15 +680,15 @@ def get_footprint(model, index):
     np.testing.assert_allclose(output_max_value, expected_max_value)
 
 
-@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."""
 
     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()