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. diff --git a/docs/roman/resample/arguments.rst b/docs/roman/resample/arguments.rst index b72bb630a..668ef208d 100644 --- a/docs/roman/resample/arguments.rst +++ b/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 diff --git a/romancal/outlier_detection/tests/test_outlier_detection.py b/romancal/outlier_detection/tests/test_outlier_detection.py index 7779dadc9..4233b8c4d 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) @@ -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" 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 0e9ad543c..f814c57f6 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, @@ -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 @@ -94,7 +96,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: @@ -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,17 +397,18 @@ 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, ) @@ -414,7 +416,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 +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) 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) diff --git a/romancal/resample/resample_utils.py b/romancal/resample/resample_utils.py index 4b8c78aa0..199aec7d4 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) @@ -111,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. @@ -163,6 +164,22 @@ 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", + 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 @@ 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): + """ + 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) diff --git a/romancal/resample/tests/test_resample.py b/romancal/resample/tests/test_resample.py index ec8a27b69..e956a5394 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, @@ -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,7 +680,7 @@ 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, @@ -687,7 +688,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()