From 72e7bb44e56c0edec1c910bd9fabb2b1e04c689a Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Mon, 9 Dec 2024 10:05:50 -0500 Subject: [PATCH 1/5] JP-3757: Use squared drizzle coeffs for variance and error arrays --- changes/8997.resample.rst | 2 + jwst/resample/resample.py | 287 +++------------------- jwst/resample/tests/test_resample_step.py | 33 ++- pyproject.toml | 1 + 4 files changed, 63 insertions(+), 260 deletions(-) create mode 100644 changes/8997.resample.rst diff --git a/changes/8997.resample.rst b/changes/8997.resample.rst new file mode 100644 index 0000000000..9315d83352 --- /dev/null +++ b/changes/8997.resample.rst @@ -0,0 +1,2 @@ +Resample variance and error arrays using squared drizzle coefficients +(error propagation). diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py index 673c4c7f07..e4dafc6315 100644 --- a/jwst/resample/resample.py +++ b/jwst/resample/resample.py @@ -1,6 +1,5 @@ import logging import os -import warnings import json import numpy as np @@ -148,7 +147,6 @@ def __init__(self, input_models, output=None, single=False, blendheaders=True, log.debug(f"Output mosaic size: {tuple(self.output_wcs.pixel_shape)}") - def _create_output_model(self, ref_input_model=None): """ Create a new blank model and update it's meta with info from ``ref_input_model``. """ output_model = datamodels.ImageModel(None) # tuple(self.output_wcs.array_shape)) @@ -245,14 +243,7 @@ def resample_group(self, input_models, indices, compute_error=False): fillval=self.fillval, disable_ctx=True, ) - # Also make a temporary model to hold error data - if compute_error: - driz_error = Drizzle( - out_shape=self.output_array_shape, - kernel=self.kernel, - fillval=self.fillval, - disable_ctx=True, - ) + log.info(f"{len(indices)} exposures to drizzle together") for index in indices: img = input_models.borrow(index) @@ -309,6 +300,7 @@ def resample_group(self, input_models, indices, compute_error=False): 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, scale=iscale, weight_map=inwht, wht_scale=1.0, # hard-coded for JWST count-rate data @@ -321,33 +313,15 @@ def resample_group(self, input_models, indices, compute_error=False): ) del data - # make an approximate error image by drizzling it - # in the same way the image is handled - if compute_error: - driz_error.add_image( - data=img.err, - exptime=img.meta.exposure.exposure_time, # GWCSDrizzle.add_image param default - pixmap=pixmap, - scale=iscale, - weight_map=inwht, - wht_scale=1.0, # hard-coded for JWST count-rate data - pixfrac=self.pixfrac, - in_units="cps", # GWCSDrizzle.add_image param default - xmin=xmin, - xmax=xmax, - ymin=ymin, - ymax=ymax, - ) input_models.shelve(img, index, modify=False) del img output_model.data = driz.out_img output_model.wht = driz.out_wht - del driz # copy the drizzled error into the output model if compute_error: - output_model.err = driz_error.out_img - del driz_error + output_model.err = np.sqrt(driz.out_img2[0]) + del driz return output_model @@ -409,13 +383,15 @@ def resample_many_to_one(self, input_models): max_ctx_id=len(input_models), disable_ctx=False, ) - self._init_variance_arrays() + self._init_exptime_counters() log.info("Resampling science and variance data") - leading_group_idx = [v[0] for v in input_models.group_indices.values()] + invalid_var = 4 * [False] + var_list = ["var_rnoise", "var_flat", "var_poisson", "err"] + leading_group_idx = [v[0] for v in input_models.group_indices.values()] with input_models: for idx, img in enumerate(input_models): if output_model is None: @@ -465,10 +441,28 @@ def resample_many_to_one(self, input_models): data.shape, ) + data2 = [] + for k, name in enumerate(var_list): + var = getattr(img, name) + if var is None or var.shape != data.shape or var.size == 0: + data2.append(None) + invalid_var[k] = True + log.warning( + f"Data shape mismatch for '{name}' for model " + f"{repr(img.meta.filename)}. Skipping ..." + ) + else: + if name == "err": + data2.append(np.square(var)) + else: + data2.append(var) + del var + driz.add_image( data=data, exptime=img.meta.exposure.exposure_time, # GWCSDrizzle.add_image param default pixmap=pixmap, + data2=data2, scale=iscale, weight_map=inwht, wht_scale=1.0, # hard-coded for JWST count-rate data @@ -479,15 +473,6 @@ def resample_many_to_one(self, input_models): ymin=ymin, ymax=ymax, ) - # Resample variance arrays in input_models to output_model - self._resample_variance_arrays( - model=img, - iscale=iscale, - inwht=inwht, - pixmap=pixmap, - in_image_limits=in_image_limits, - output_shape=self.output_array_shape, - ) del data, inwht @@ -499,24 +484,16 @@ 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 - + 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] + setattr(output_model, name, var) del driz - # compute final variances: - self._compute_resample_variance_totals(output_model) - - var_components = [ - output_model.var_rnoise, - output_model.var_poisson, - output_model.var_flat - ] - output_model.err = np.sqrt(np.nansum(var_components, axis=0)) - - # nansum returns zero for input that is all NaN - - # set those values to NaN instead - all_nan = np.all(np.isnan(var_components), axis=0) - output_model.err[all_nan] = np.nan - if self.blendheaders: blender.finalize_model(output_model) @@ -524,202 +501,6 @@ def resample_many_to_one(self, input_models): return ModelLibrary([output_model,], on_disk=False) - def _init_variance_arrays(self): - shape = self.output_array_shape - self._weighted_rn_var = np.full(shape, np.nan, dtype=np.float32) - self._weighted_pn_var = np.full(shape, np.nan, dtype=np.float32) - self._weighted_flat_var = np.full(shape, np.nan, dtype=np.float32) - self._total_weight_rn_var = np.zeros(shape, dtype=np.float32) - self._total_weight_pn_var = np.zeros(shape, dtype=np.float32) - self._total_weight_flat_var = np.zeros(shape, dtype=np.float32) - - def _resample_variance_arrays(self, model, iscale, inwht, pixmap, - in_image_limits, output_shape): - xmin, xmax, ymin, ymax = in_image_limits - - # Do the read noise variance first, so it can be - # used for weights if needed - rn_var = self._resample_one_variance_array( - "var_rnoise", - input_model=model, - iscale=iscale, - inwht=inwht, - pixmap=pixmap, - xmin=xmin, - xmax=xmax, - ymin=ymin, - ymax=ymax, - ) - - # Find valid weighting values in the variance - if rn_var is not None: - mask = (rn_var > 0) & np.isfinite(rn_var) - else: - mask = np.full_like(rn_var, False) - - # Set the weight for the image from the weight type - weight = np.ones(output_shape) - if self.weight_type == "ivm" and rn_var is not None: - weight[mask] = rn_var[mask] ** -1 - elif self.weight_type == "exptime": - if resample_utils.check_for_tmeasure(model): - weight[:] = model.meta.exposure.measurement_time - else: - weight[:] = model.meta.exposure.exposure_time - - # Weight and add the readnoise variance - # Note: floating point overflow is an issue if variance weights - # are used - it can't be squared before multiplication - if rn_var is not None: - mask = (rn_var >= 0) & np.isfinite(rn_var) & (weight > 0) - self._weighted_rn_var[mask] = np.nansum( - [ - self._weighted_rn_var[mask], - rn_var[mask] * weight[mask] * weight[mask] - ], - axis=0 - ) - self._total_weight_rn_var[mask] += weight[mask] - - # Now do poisson and flat variance, updating only valid new values - # (zero is a valid value; negative, inf, or NaN are not) - pn_var = self._resample_one_variance_array( - "var_poisson", - input_model=model, - iscale=iscale, - inwht=inwht, - pixmap=pixmap, - xmin=xmin, - xmax=xmax, - ymin=ymin, - ymax=ymax, - ) - if pn_var is not None: - mask = (pn_var >= 0) & np.isfinite(pn_var) & (weight > 0) - self._weighted_pn_var[mask] = np.nansum( - [ - self._weighted_pn_var[mask], - pn_var[mask] * weight[mask] * weight[mask] - ], - axis=0 - ) - self._total_weight_pn_var[mask] += weight[mask] - - flat_var = self._resample_one_variance_array( - "var_flat", - input_model=model, - iscale=iscale, - inwht=inwht, - pixmap=pixmap, - xmin=xmin, - xmax=xmax, - ymin=ymin, - ymax=ymax, - ) - if flat_var is not None: - mask = (flat_var >= 0) & np.isfinite(flat_var) & (weight > 0) - self._weighted_flat_var[mask] = np.nansum( - [ - self._weighted_flat_var[mask], - flat_var[mask] * weight[mask] * weight[mask] - ], - axis=0 - ) - self._total_weight_flat_var[mask] += weight[mask] - - def _compute_resample_variance_totals(self, output_model): - # Divide by the total weights, squared, and set in the output model. - # Zero weight and missing values are NaN in the output. - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", "invalid value*", RuntimeWarning) - warnings.filterwarnings("ignore", "divide by zero*", RuntimeWarning) - - output_variance = ( - self._weighted_rn_var / self._total_weight_rn_var / - self._total_weight_rn_var - ) - setattr(output_model, "var_rnoise", output_variance) - - output_variance = ( - self._weighted_pn_var / self._total_weight_pn_var / - self._total_weight_pn_var - ) - setattr(output_model, "var_poisson", output_variance) - - output_variance = ( - self._weighted_flat_var / self._total_weight_flat_var / - self._total_weight_flat_var - ) - setattr(output_model, "var_flat", output_variance) - - del ( - self._weighted_rn_var, - self._weighted_pn_var, - self._weighted_flat_var, - self._total_weight_rn_var, - self._total_weight_pn_var, - self._total_weight_flat_var, - ) - - def _resample_one_variance_array(self, name, input_model, iscale, - inwht, pixmap, - xmin=None, xmax=None, ymin=None, ymax=None): - """Resample one variance image from an input model. - - The error image is passed to drizzle instead of the variance, to - better match kernel overlap and user weights to the data, in the - pixel averaging process. The drizzled error image is squared before - returning. - """ - variance = getattr(input_model, name) - if variance is None or variance.size == 0: - log.debug( - f"No data for '{name}' for model " - f"{repr(input_model.meta.filename)}. Skipping ..." - ) - return - - elif variance.shape != input_model.data.shape: - log.warning( - f"Data shape mismatch for '{name}' for model " - f"{repr(input_model.meta.filename)}. Skipping ..." - ) - return - - output_shape = self.output_array_shape - - # Resample the error array. Fill "unpopulated" pixels with NaNs. - driz = Drizzle( - out_shape=output_shape, - kernel=self.kernel, - fillval=np.nan, - disable_ctx=True - ) - - log.debug(f"Pixmap shape: {pixmap[:,:,0].shape}") - log.debug(f"Input Sci shape: {variance.shape}") - log.debug(f"Output Sci shape: {output_shape}") - - # Call 'drizzle' to perform image combination - log.info(f"Drizzling {variance.shape} --> {output_shape}") - - driz.add_image( - data=np.sqrt(variance), - exptime=input_model.meta.exposure.exposure_time, - pixmap=pixmap, - scale=iscale, - weight_map=inwht, - wht_scale=1.0, # hard-coded for JWST count-rate data - pixfrac=self.pixfrac, - in_units="cps", - xmin=xmin, - xmax=xmax, - ymin=ymin, - ymax=ymax, - ) - - return driz.out_img ** 2 - def _init_exptime_counters(self): self._total_exposure_time = 0. self._exposure_times = {'start': [], 'end': []} diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py index 75eddb8530..53e3a6e1db 100644 --- a/jwst/resample/tests/test_resample_step.py +++ b/jwst/resample/tests/test_resample_step.py @@ -788,19 +788,38 @@ def test_resample_variance(nircam_rate, n_images, weight_type): var_poisson = 0.00025 im = AssignWcsStep.call(nircam_rate) _set_photom_kwd(im) + im.data[:, :] = 1.0 im.var_rnoise += var_rnoise im.var_poisson += var_poisson im.err += err im.meta.filename = "foo.fits" + c1 = ModelLibrary([im.copy()]) c = ModelLibrary([im.copy() for _ in range(n_images)]) + result1 = ResampleStep.call(c1, blendheaders=False, weight_type=weight_type) result = ResampleStep.call(c, blendheaders=False, weight_type=weight_type) # Verify that the combined uncertainty goes as 1 / sqrt(N) - assert_allclose(result.err[5:-5, 5:-5].mean(), err / np.sqrt(n_images), atol=1e-5) - assert_allclose(result.var_rnoise[5:-5, 5:-5].mean(), var_rnoise / n_images, atol=1e-7) - assert_allclose(result.var_poisson[5:-5, 5:-5].mean(), var_poisson / n_images, atol=1e-7) + mask = np.isfinite(result.err) + + assert np.all((result1.err[mask] / err) <= 1.0) + assert_allclose( + result.err[mask].mean(), + result1.err[mask].mean() / np.sqrt(n_images), atol=1e-5 + ) + + assert np.all((result1.var_rnoise[mask] / var_rnoise) <= 1.0) + assert_allclose( + result.var_rnoise[mask].mean(), + result1.var_rnoise[mask].mean() / n_images, atol=1e-5 + ) + + assert np.all((result1.var_poisson[mask] / var_poisson) <= 1.0) + assert_allclose( + result.var_poisson[mask].mean(), + result1.var_poisson[mask].mean() / n_images, atol=1e-5 + ) im.close() result.close() @@ -819,8 +838,8 @@ def test_resample_undefined_variance(nircam_rate, shape): with pytest.warns(RuntimeWarning, match="var_rnoise array not available"): result = ResampleStep.call(c, blendheaders=False) - # no valid variance - output error and variance are all NaN - assert_allclose(result.err, np.nan) + # no valid variance - variance are all NaN + assert np.any(np.isfinite(result.err)) assert_allclose(result.var_rnoise, np.nan) assert_allclose(result.var_poisson, np.nan) assert_allclose(result.var_flat, np.nan) @@ -928,8 +947,8 @@ def test_custom_refwcs_resample_imaging(nircam_rate, output_shape2, match, input_mean = np.nanmean(im.data) output_mean_1 = np.nanmean(data1) output_mean_2 = np.nanmean(data2) - assert np.isclose(input_mean * iscale**2, output_mean_1, atol=1e-4) - assert np.isclose(input_mean * iscale**2, output_mean_2, atol=1e-4) + assert np.isclose(input_mean * iscale**2, output_mean_1, atol=2e-4) + assert np.isclose(input_mean * iscale**2, output_mean_2, atol=2e-4) im.close() result.close() diff --git a/pyproject.toml b/pyproject.toml index ef138c0999..72411cbd28 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,7 @@ dependencies = [ "BayesicFitting>=3.0.1", "crds>=12.0.3", "drizzle>=2.0.0", + # "drizzle @ git+https://github.com/mcara/drizzle.git@add-resampling-with-square-weights", "gwcs>=0.21.0,<0.23.0", "numpy>=1.22,<2.0", "opencv-python-headless>=4.6.0.66", From bc372bd0fbdbe5ae55ed108690532fb1771126c3 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Mon, 9 Dec 2024 10:48:13 -0500 Subject: [PATCH 2/5] temporarily change drizzle dependency --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 72411cbd28..5c3be5f1ce 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,8 +22,8 @@ dependencies = [ "astropy>=5.3", "BayesicFitting>=3.0.1", "crds>=12.0.3", - "drizzle>=2.0.0", - # "drizzle @ git+https://github.com/mcara/drizzle.git@add-resampling-with-square-weights", + # "drizzle>=2.1.0", + "drizzle @ git+https://github.com/mcara/drizzle.git@add-resampling-with-square-weights", "gwcs>=0.21.0,<0.23.0", "numpy>=1.22,<2.0", "opencv-python-headless>=4.6.0.66", From d23e4a1daa0bd55ecf18814ea47cdeb6c1b20c48 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Thu, 12 Dec 2024 04:11:25 -0500 Subject: [PATCH 3/5] add unit test for varying pixel scale ratio --- docs/jwst/resample/arguments.rst | 2 +- docs/jwst/resample/main.rst | 48 +-------- jwst/resample/resample.py | 28 ++++-- jwst/resample/tests/test_resample_step.py | 117 ++++++++++++++++++++-- 4 files changed, 134 insertions(+), 61 deletions(-) diff --git a/docs/jwst/resample/arguments.rst b/docs/jwst/resample/arguments.rst index 2b2a47a6a0..5f1b5b5d6d 100644 --- a/docs/jwst/resample/arguments.rst +++ b/docs/jwst/resample/arguments.rst @@ -95,7 +95,7 @@ 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 measurement time (TMEASURE) found in the image header if available; if unavailable, the scaling will be set equal to the exposure time (EFFEXPTM). diff --git a/docs/jwst/resample/main.rst b/docs/jwst/resample/main.rst index 613d3e167a..10c0349870 100644 --- a/docs/jwst/resample/main.rst +++ b/docs/jwst/resample/main.rst @@ -36,56 +36,18 @@ drizzling to create the output product. Error Propagation ----------------- -The error associated with each resampled pixel can in principle be derived +The error associated with each resampled pixel is derived from the variance components associated with each input pixel, weighted by -the square of the input user weights and the square of the overlap between -the input and output pixels. In practice, the cdriz routine does not currently -support propagating variance data alongside science images, so the output -error cannot be precisely calculated. - -To approximate the error on a resampled pixel, the variance arrays associated -with each input image are resampled individually, then combined with a weighted -sum. The process is: - -#. For each input image, take the square root of each of the read noise variance - array to make an error image. - -#. Drizzle the read noise error image onto the output WCS, with drizzle - parameters matching those used for the science data. - -#. Square the resampled read noise to make a variance array. - - a. If the resampling `weight_type` is an inverse variance map (`ivm`), weight - the resampled variance by the square of its own inverse. - - #. If the `weight_type` is the exposure time (`exptime`), weight the - resampled variance by the square of the exposure time for the image. - -#. Add the weighted, resampled read noise variance to a running sum across all - images. Add the weights (unsquared) to a separate running sum across - all images. - -#. Perform the same steps for the Poisson noise variance and the flat variance. - For these components, the weight for the sum is either the resampled read - noise variance or else the exposure time. - -#. For each variance component (read noise, Poisson, and flat), divide the - weighted variance sum by the total weight, squared. +the square of the input user weights, pixel scale ratio, and the square of +the overlap between the input and output pixels. Essentially the code performs +[propagation of uncertainty](https://en.wikipedia.org/wiki/Propagation_of_uncertainty) +for each of the variance arrays in the input models. After each variance component is resampled and summed, the final error array is computed as the square root of the sum of the three independent variance components. This error image is stored in the ``err`` attribute in the output data model or the ``'ERR'`` extension of the FITS file. -It is expected that the output errors computed in this way will -generally overestimate the true error on the resampled data. The magnitude -of the overestimation depends on the details of the pixel weights -and error images. Note, however, that drizzling error images produces -a much better estimate of the output error than directly drizzling -the variance images, since the kernel overlap weights do not need to be -squared for combining error values. - - Context Image ------------- 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..7c4ca1d798 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 @@ -803,22 +804,31 @@ def test_resample_variance(nircam_rate, n_images, weight_type): # Verify that the combined uncertainty goes as 1 / sqrt(N) mask = np.isfinite(result.err) + twht = np.sum(result.wht[mask]) + twht1 = np.sum(result1.wht[mask]) + assert np.all((result1.err[mask] / err) <= 1.0) assert_allclose( - result.err[mask].mean(), - result1.err[mask].mean() / np.sqrt(n_images), atol=1e-5 + np.sum(result.err[mask]**2) * twht**2, + np.sum(result1.err[mask]**2) * twht1**2, + rtol=1e-6, + atol=0.0, ) assert np.all((result1.var_rnoise[mask] / var_rnoise) <= 1.0) assert_allclose( - result.var_rnoise[mask].mean(), - result1.var_rnoise[mask].mean() / n_images, atol=1e-5 + np.sum(result.var_rnoise[mask]) * twht, + np.sum(result1.var_rnoise[mask]) * twht1, + rtol=1e-6, + atol=0.0, ) assert np.all((result1.var_poisson[mask] / var_poisson) <= 1.0) assert_allclose( - result.var_poisson[mask].mean(), - result1.var_poisson[mask].mean() / n_images, atol=1e-5 + np.sum(result.var_poisson[mask]) * twht, + np.sum(result1.var_poisson[mask]) * twht1, + rtol=1e-6, + atol=0.0, ) im.close() @@ -1454,3 +1464,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() From 7531f55118be700549ba661c623b88e0e8fa7624 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Thu, 12 Dec 2024 09:49:35 -0500 Subject: [PATCH 4/5] Fix err computation for undefined variances --- jwst/resample/resample.py | 7 ++++++- jwst/resample/tests/test_resample_step.py | 6 +++--- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py index 356748c4e2..9a5c883b55 100644 --- a/jwst/resample/resample.py +++ b/jwst/resample/resample.py @@ -496,7 +496,12 @@ def resample_many_to_one(self, input_models): valid_var.append(var) setattr(output_model, name, var) - err = np.sum(valid_var, axis=0).astype(np.float32) + if len(valid_var) > 0: + all_nan = np.all(np.isnan(valid_var), axis=0) + err = np.sqrt(np.nansum(valid_var, axis=0)).astype(np.float32) + err[all_nan] = np.nan + else: + err = np.full_like(output_model.data, np.nan) output_model.err = err del driz, err, valid_var, var diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py index 7c4ca1d798..a9c51f9667 100644 --- a/jwst/resample/tests/test_resample_step.py +++ b/jwst/resample/tests/test_resample_step.py @@ -809,8 +809,8 @@ def test_resample_variance(nircam_rate, n_images, weight_type): assert np.all((result1.err[mask] / err) <= 1.0) assert_allclose( - np.sum(result.err[mask]**2) * twht**2, - np.sum(result1.err[mask]**2) * twht1**2, + np.sum(result.err[mask]**2) * twht, + np.sum(result1.err[mask]**2) * twht1, rtol=1e-6, atol=0.0, ) @@ -849,7 +849,7 @@ def test_resample_undefined_variance(nircam_rate, shape): result = ResampleStep.call(c, blendheaders=False) # no valid variance - variance are all NaN - assert np.any(np.isfinite(result.err)) + assert_allclose(result.err, np.nan) assert_allclose(result.var_rnoise, np.nan) assert_allclose(result.var_poisson, np.nan) assert_allclose(result.var_flat, np.nan) From a29c7b20244db7cb6011fbe1504f749a541738e9 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Fri, 13 Dec 2024 13:27:06 -0500 Subject: [PATCH 5/5] Remove outdated code comments --- jwst/resample/tests/test_resample_step.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py index a9c51f9667..b1f4366ffc 100644 --- a/jwst/resample/tests/test_resample_step.py +++ b/jwst/resample/tests/test_resample_step.py @@ -1482,9 +1482,6 @@ def test_nirspec_lamp_pixscale(nirspec_lamp, tmp_path): ) ) 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) @@ -1519,7 +1516,6 @@ def test_variance_arrays(kernel_fc, ps_ratio, weights, nircam_rate): 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,