From cc680a3778624aef4edcf7a0d1e7b4a0c102e8a2 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Thu, 12 Dec 2024 04:11:25 -0500 Subject: [PATCH] add unit test for varying pixel scale ratio --- jwst/resample/resample.py | 28 ++++--- jwst/resample/tests/test_resample_step.py | 96 +++++++++++++++++++++++ 2 files changed, 113 insertions(+), 11 deletions(-) diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py index e4dafc6315..356748c4e2 100644 --- a/jwst/resample/resample.py +++ b/jwst/resample/resample.py @@ -296,11 +296,16 @@ def resample_group(self, input_models, indices, compute_error=False): img.data.shape, ) + if compute_error and img.err is not None: + data2 = np.square(img.err) + else: + data2 = None + driz.add_image( data=data, exptime=img.meta.exposure.exposure_time, # GWCSDrizzle.add_image param default was 1.0 pixmap=pixmap, - data2=np.square(img.err) if compute_error else None, + data2=data2, scale=iscale, weight_map=inwht, wht_scale=1.0, # hard-coded for JWST count-rate data @@ -319,7 +324,7 @@ def resample_group(self, input_models, indices, compute_error=False): output_model.data = driz.out_img output_model.wht = driz.out_wht # copy the drizzled error into the output model - if compute_error: + if compute_error and driz.out_img2 is not None: output_model.err = np.sqrt(driz.out_img2[0]) del driz @@ -388,8 +393,8 @@ def resample_many_to_one(self, input_models): log.info("Resampling science and variance data") - invalid_var = 4 * [False] - var_list = ["var_rnoise", "var_flat", "var_poisson", "err"] + var_list = ["var_rnoise", "var_flat", "var_poisson"] + invalid_var = len(var_list) * [False] leading_group_idx = [v[0] for v in input_models.group_indices.values()] with input_models: @@ -452,10 +457,7 @@ def resample_many_to_one(self, input_models): f"{repr(img.meta.filename)}. Skipping ..." ) else: - if name == "err": - data2.append(np.square(var)) - else: - data2.append(var) + data2.append(var) del var driz.add_image( @@ -484,15 +486,19 @@ def resample_many_to_one(self, input_models): output_model.wht = driz.out_wht if driz.out_ctx is not None: output_model.con = driz.out_ctx + + valid_var = [] for k, name in enumerate(var_list): if invalid_var[k]: var = np.full_like(output_model.data, np.nan) - elif name == "err": - var = np.sqrt(driz.out_img2[k]) else: var = driz.out_img2[k] + valid_var.append(var) setattr(output_model, name, var) - del driz + + err = np.sum(valid_var, axis=0).astype(np.float32) + output_model.err = err + del driz, err, valid_var, var if self.blendheaders: blender.finalize_model(output_model) diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py index 53e3a6e1db..ad3f8ff39d 100644 --- a/jwst/resample/tests/test_resample_step.py +++ b/jwst/resample/tests/test_resample_step.py @@ -1,4 +1,5 @@ import pytest +from itertools import product from gwcs.wcstools import grid_from_bounding_box from numpy.testing import assert_allclose @@ -1454,3 +1455,98 @@ def test_nirspec_lamp_pixscale(nirspec_lamp, tmp_path): result2.close() result3.close() result4.close() + + +@pytest.mark.filterwarnings("ignore:Kernel '") +@pytest.mark.parametrize( + 'kernel_fc, ps_ratio, weights', + ( + x for x in product( + [ + ('square', True), + ('point', True), + ('gaussian', False), + ], + [0.25, 0.5, 1, 1.2], + [(0.99, 0.01), (0.9, 1.5), (467, 733)], + ) + ) +) +def test_variance_arrays(kernel_fc, ps_ratio, weights, nircam_rate): + + # check that if both 'pixel_scale_ratio' and 'pixel_scale' are passed in, + # that 'pixel_scale' overrides correctly + im1 = AssignWcsStep.call(nircam_rate, sip_approx=False) + _set_photom_kwd(im1) + im2 = AssignWcsStep.call(nircam_rate, sip_approx=False) + _set_photom_kwd(im2) + + shape = im1.data.shape + xc = shape[1] // 2 + yc = shape[0] // 2 + + # unpack parameters: + kernel, fc = kernel_fc + + # pixel values in input data: + dataval = [1.0, 7.0] + + # pixel values in input variance: + varval = [0.5, 50.0] + + sl = np.s_[yc - 4: yc + 5, xc - 4: xc + 5] + im1.data[sl] = dataval[0] + im2.data[sl] = dataval[1] + + im1.var_poisson[:, :] = 0.0 + im1.var_poisson[sl] = varval[0] + im2.var_poisson[:, :] = 0.0 + im2.var_poisson[sl] = varval[1] + + im1.meta.exposure.exposure_time = weights[0] + im1.meta.exposure.measurement_time = weights[0] + im2.meta.exposure.exposure_time = weights[1] + im2.meta.exposure.measurement_time = weights[1] + + library = ModelLibrary([im1, im2]) + + # check when both pixel_scale and pixel_scale_ratio are passed in + res = ResampleStep.call( + library, + pixel_scale_ratio=ps_ratio, + kernel=kernel, + weight_type="exptime", + blendheaders=False, + save_results=True, + ) + + assert np.any(np.isfinite(res.var_poisson)) + + mask = res.con[0] > 0 + n_nonzero = np.sum(im1.data > 0.0) + res.var_poisson[np.logical_not(mask)] = 0.0 + + rtol = 1.0e-6 if fc else 0.15 + + ideal_output = np.dot(dataval, weights) * n_nonzero + ideal_output2 = np.dot(varval, np.square(weights)) / np.sum(weights)**2 + + tflux = np.sum(res.data[mask] * res.wht[mask]) + tflux2 = np.max(res.var_poisson) + + # check output flux: + assert np.allclose( + tflux, + ideal_output, + rtol=rtol, + atol=0.0 + ) + + # check output variance: + # less restrictive (to account for pixel overlap variations): + assert (np.max(tflux2) <= ideal_output2 * (1 + rtol) and + np.max(tflux2) >= 0.25 * ideal_output2 * (1 - rtol)) + + im1.close() + im2.close() + res.close()