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/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 673c4c7f07..9a5c883b55 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) @@ -305,10 +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=data2, scale=iscale, weight_map=inwht, wht_scale=1.0, # hard-coded for JWST count-rate data @@ -321,33 +318,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 + if compute_error and driz.out_img2 is not None: + output_model.err = np.sqrt(driz.out_img2[0]) + del driz return output_model @@ -409,13 +388,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()] + 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: for idx, img in enumerate(input_models): if output_model is None: @@ -465,10 +446,25 @@ 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: + 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 +475,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 @@ -500,22 +487,23 @@ def resample_many_to_one(self, input_models): if driz.out_ctx is not None: output_model.con = driz.out_ctx - 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 + valid_var = [] + for k, name in enumerate(var_list): + if invalid_var[k]: + var = np.full_like(output_model.data, np.nan) + else: + var = driz.out_img2[k] + valid_var.append(var) + setattr(output_model, name, var) + + 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 if self.blendheaders: blender.finalize_model(output_model) @@ -524,202 +512,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..b1f4366ffc 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 @@ -788,19 +789,47 @@ 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) + + twht = np.sum(result.wht[mask]) + twht1 = np.sum(result1.wht[mask]) + + assert np.all((result1.err[mask] / err) <= 1.0) + assert_allclose( + np.sum(result.err[mask]**2) * twht, + np.sum(result1.err[mask]**2) * twht1, + rtol=1e-6, + atol=0.0, + ) + + assert np.all((result1.var_rnoise[mask] / var_rnoise) <= 1.0) + assert_allclose( + 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( + np.sum(result.var_poisson[mask]) * twht, + np.sum(result1.var_poisson[mask]) * twht1, + rtol=1e-6, + atol=0.0, + ) im.close() result.close() @@ -819,7 +848,7 @@ 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 + # no valid variance - variance are all NaN assert_allclose(result.err, np.nan) assert_allclose(result.var_rnoise, np.nan) assert_allclose(result.var_poisson, np.nan) @@ -928,8 +957,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() @@ -1435,3 +1464,94 @@ 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): + 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]) + + 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() diff --git a/pyproject.toml b/pyproject.toml index ef138c0999..5c3be5f1ce 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,7 +22,8 @@ dependencies = [ "astropy>=5.3", "BayesicFitting>=3.0.1", "crds>=12.0.3", - "drizzle>=2.0.0", + # "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",