diff --git a/CHANGES.rst b/CHANGES.rst index 9ebe78b840..0b6da3affb 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -13,6 +13,12 @@ associations - Constraint added to Asn_Lv3Coron to remove background exposures [#5781] +resample +-------- + +- Make inverse variance ``weight_type="ivm"`` the default weighting scheme for + multiple exposures resampled into a single output. [#5738] + 1.0.0 (2021-02-22) ================== diff --git a/JenkinsfileRT b/JenkinsfileRT index 00e92c7029..57e31b37e0 100644 --- a/JenkinsfileRT +++ b/JenkinsfileRT @@ -82,7 +82,6 @@ bc0.build_cmds = [ "pip install -e .[test,ephem]", "pip install pytest-xdist pytest-sugar", "pip freeze", - 'echo "CRDS_CONTEXT = $(crds list --contexts $CRDS_CONTEXT --mappings | grep pmap)"', ] bc0.build_cmds = PipInject(env.OVERRIDE_REQUIREMENTS) + bc0.build_cmds bc0.test_cmds = [ diff --git a/JenkinsfileRT_dev b/JenkinsfileRT_dev index 998f437958..3a7a70890a 100644 --- a/JenkinsfileRT_dev +++ b/JenkinsfileRT_dev @@ -83,7 +83,6 @@ bc0.build_cmds = [ "pip install -e .[test,ephem]", "pip install pytest-xdist pytest-sugar", "pip freeze", - 'echo "CRDS_CONTEXT = $(crds list --contexts $CRDS_CONTEXT --mappings | grep pmap)"', ] bc0.build_cmds = PipInject(env.OVERRIDE_REQUIREMENTS) + bc0.build_cmds bc0.test_cmds = [ diff --git a/jwst/datamodels/schemas/core.schema.yaml b/jwst/datamodels/schemas/core.schema.yaml index a1b9215805..910e7bc683 100644 --- a/jwst/datamodels/schemas/core.schema.yaml +++ b/jwst/datamodels/schemas/core.schema.yaml @@ -2324,7 +2324,7 @@ properties: weight_type: title: Type of drizzle weighting to use in resampling input type: string - enum: [exptime, error] + enum: [exptime, ivm] fits_keyword: RESWHT background: title: Background information diff --git a/jwst/outlier_detection/outlier_detection.py b/jwst/outlier_detection/outlier_detection.py index d59edc2aab..6ff796cc33 100644 --- a/jwst/outlier_detection/outlier_detection.py +++ b/jwst/outlier_detection/outlier_detection.py @@ -110,7 +110,7 @@ def _convert_inputs(self): dq=self.inputs.dq[i]) image.meta = self.inputs.meta image.wht = build_driz_weight(image, - weight_type='exptime', + weight_type='ivm', good_bits=bits) self.input_models.append(image) self.converted = True @@ -202,7 +202,7 @@ def do_detection(self): for i in range(len(self.input_models)): drizzled_models[i].wht = build_driz_weight( self.input_models[i], - weight_type='exptime', + weight_type='ivm', good_bits=pars['good_bits']) # Initialize intermediate products used in the outlier detection diff --git a/jwst/outlier_detection/outlier_detection_scaled.py b/jwst/outlier_detection/outlier_detection_scaled.py index 48c605e20c..496d19eb31 100644 --- a/jwst/outlier_detection/outlier_detection_scaled.py +++ b/jwst/outlier_detection/outlier_detection_scaled.py @@ -106,7 +106,7 @@ def do_detection(self): for image in self.input_models: image.wht = resample_utils.build_driz_weight( image, - weight_type='exptime', + weight_type='ivm', good_bits=pars['good_bits'] ) diff --git a/jwst/outlier_detection/outlier_detection_scaled_step.py b/jwst/outlier_detection/outlier_detection_scaled_step.py index 4d5a8c100f..9e900ac94f 100755 --- a/jwst/outlier_detection/outlier_detection_scaled_step.py +++ b/jwst/outlier_detection/outlier_detection_scaled_step.py @@ -19,7 +19,7 @@ class OutlierDetectionScaledStep(Step): """ spec = """ - weight_type = option('exptime','error',None,default='exptime') + weight_type = option('ivm','exptime',default='ivm') pixfrac = float(default=1.0) kernel = string(default='square') # drizzle kernel fillval = string(default='INDEF') diff --git a/jwst/outlier_detection/outlier_detection_spec.py b/jwst/outlier_detection/outlier_detection_spec.py index 4a6633de17..b2a7d30c16 100644 --- a/jwst/outlier_detection/outlier_detection_spec.py +++ b/jwst/outlier_detection/outlier_detection_spec.py @@ -85,7 +85,7 @@ def do_detection(self): for i in range(len(self.input_models)): drizzled_models[i].wht = resample_utils.build_driz_weight( self.input_models[i], - weight_type='exptime', + weight_type='ivm', good_bits=pars['good_bits']) # Initialize intermediate products used in the outlier detection diff --git a/jwst/outlier_detection/outlier_detection_stack_step.py b/jwst/outlier_detection/outlier_detection_stack_step.py index 07d22b1bf0..a379b61538 100644 --- a/jwst/outlier_detection/outlier_detection_stack_step.py +++ b/jwst/outlier_detection/outlier_detection_stack_step.py @@ -28,7 +28,7 @@ class OutlierDetectionStackStep(Step): """ spec = """ - weight_type = option('exptime','error',None,default='exptime') + weight_type = option('ivm','exptime',default='ivm') pixfrac = float(default=1.0) kernel = string(default='square') # drizzle kernel fillval = string(default='INDEF') diff --git a/jwst/outlier_detection/outlier_detection_step.py b/jwst/outlier_detection/outlier_detection_step.py index 61de9bbf23..7e1a83f4fc 100644 --- a/jwst/outlier_detection/outlier_detection_step.py +++ b/jwst/outlier_detection/outlier_detection_step.py @@ -46,7 +46,7 @@ class OutlierDetectionStep(Step): # by the various versions of the outlier_detection algorithms, and each # version will pick and choose what they need while ignoring the rest. spec = """ - weight_type = option('exptime','error',None,default='exptime') + weight_type = option('ivm','exptime',default='ivm') pixfrac = float(default=1.0) kernel = string(default='square') # drizzle kernel fillval = string(default='INDEF') diff --git a/jwst/resample/gwcs_drizzle.py b/jwst/resample/gwcs_drizzle.py index 34e2461983..a8b8494339 100644 --- a/jwst/resample/gwcs_drizzle.py +++ b/jwst/resample/gwcs_drizzle.py @@ -1,7 +1,6 @@ import numpy as np from drizzle import util -from drizzle import doblot from drizzle import cdrizzle from . import resample_utils @@ -13,9 +12,8 @@ class GWCSDrizzle: """ Combine images using the drizzle algorithm """ - def __init__(self, product, outwcs=None, single=False, - wt_scl="exptime", pixfrac=1.0, kernel="square", - fillval="INDEF"): + def __init__(self, product, outwcs=None, single=False, wt_scl=None, + pixfrac=1.0, kernel="square", fillval="INDEF"): """ Create a new Drizzle output object and set the drizzle parameters. @@ -33,10 +31,10 @@ def __init__(self, product, outwcs=None, single=False, provided, the WCS is taken from product. wt_scl : str, optional - How each input image should be scaled. The choices are `exptime` - which scales each image by its exposure time, `expsq` which scales - each image by the exposure time squared, or an empty string, which - allows each input image to be scaled individually. + How each input image should be scaled. The choices are `exptime`, + which scales each image by its exposure time, or `expsq`, which scales + each image by the exposure time squared. If not set, then each + input image is scaled by its own weight map. pixfrac : float, optional The fraction of a pixel that the pixel flux is confined to. The @@ -64,7 +62,10 @@ def __init__(self, product, outwcs=None, single=False, self.outexptime = 0.0 self.uniqid = 0 - self.wt_scl = wt_scl + if wt_scl is None: + self.wt_scl = "" + else: + self.wt_scl = wt_scl self.kernel = kernel self.fillval = fillval self.pixfrac = pixfrac @@ -90,11 +91,7 @@ def __init__(self, product, outwcs=None, single=False, self.outcon = np.reshape(self.outcon, (1, self.outcon.shape[0], self.outcon.shape[1])) - - elif self.outcon.ndim == 3: - pass - - else: + elif self.outcon.ndim != 3: raise ValueError("Drizzle context image has wrong dimensions: \ {0}".format(product)) @@ -102,11 +99,6 @@ def __init__(self, product, outwcs=None, single=False, if not self.outwcs: raise ValueError("Either an existing file or wcs must be supplied") - if util.is_blank(self.wt_scl): - self.wt_scl = '' - elif self.wt_scl != "exptime" and self.wt_scl != "expsq": - raise ValueError("Illegal value for wt_scl: %s" % self.wt_scl) - if out_units == "counts": np.divide(self.outsci, self.outexptime, self.outsci) elif out_units != "cps": @@ -199,34 +191,6 @@ def add_image(self, insci, inwcs, inwht=None, xmin=0, xmax=0, ymin=0, ymax=0, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, pixfrac=self.pixfrac, kernel=self.kernel, fillval=self.fillval) - def blot_image(self, blotwcs, interp='poly5', sinscl=1.0): - """ - Resample the output image using an input world coordinate system. - - Parameters - ---------- - - blotwcs : wcs - The world coordinate system to resample on. - - interp : str, optional - The type of interpolation used in the resampling. The - possible values are "nearest" (nearest neighbor interpolation), - "linear" (bilinear interpolation), "poly3" (cubic polynomial - interpolation), "poly5" (quintic polynomial interpolation), - "sinc" (sinc interpolation), "lan3" (3rd order Lanczos - interpolation), and "lan5" (5th order Lanczos interpolation). - - sincscl : float, optional - The scaling factor for sinc interpolation. - """ - - util.set_pscale(blotwcs) - self.outsci = doblot.doblot(self.outsci, self.outwcs, blotwcs, - 1.0, interp=interp, sinscl=sinscl) - - self.outwcs = blotwcs - def increment_id(self): """ Increment the id count and add a plane to the context image if needed diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py index 772c032625..74fb5bf082 100644 --- a/jwst/resample/resample.py +++ b/jwst/resample/resample.py @@ -41,7 +41,7 @@ class ResampleData: (eventually) a record of metadata from all input models. """ def __init__(self, input_models, output=None, single=False, blendheaders=True, - pixfrac=1.0, kernel="square", fillval="INDEF", weight_type="exptime", + pixfrac=1.0, kernel="square", fillval="INDEF", weight_type="ivm", good_bits=0, pscale_ratio=1.0, **kwargs): """ Parameters @@ -147,8 +147,7 @@ def do_drizzle(self): single=self.single, pixfrac=self.pixfrac, kernel=self.kernel, - fillval=self.fillval, - wt_scl=self.weight_type) + fillval=self.fillval) for n, img_exp in enumerate(exposure): # Make a copy as this is operated on in-place below @@ -165,7 +164,7 @@ def do_drizzle(self): weight_type=self.weight_type, good_bits=self.good_bits) driz.add_image(img.data, img.meta.wcs, inwht=inwht, - expin=img.meta.exposure.exposure_time) + expin=img.meta.exposure.exposure_time) # Update some basic exposure time values based on all the inputs output_model.meta.exposure.exposure_time = texptime diff --git a/jwst/resample/resample_spec.py b/jwst/resample/resample_spec.py index 5d38afd4d8..6a84fc2a05 100644 --- a/jwst/resample/resample_spec.py +++ b/jwst/resample/resample_spec.py @@ -46,7 +46,7 @@ class ResampleSpecData: def __init__(self, input_models, output=None, single=False, blendheaders=False, pixfrac=1.0, kernel="square", fillval=0, - weight_type="exptime", good_bits=0, pscale_ratio=1.0, **kwargs): + weight_type="ivm", good_bits=0, pscale_ratio=1.0, **kwargs): """ Parameters ---------- diff --git a/jwst/resample/resample_step.py b/jwst/resample/resample_step.py index a5ac0b655b..1692dc4418 100755 --- a/jwst/resample/resample_step.py +++ b/jwst/resample/resample_step.py @@ -34,7 +34,7 @@ class ResampleStep(Step): pixfrac = float(default=1.0) kernel = string(default='square') fillval = string(default='INDEF') - weight_type = option('exptime', default='exptime') + weight_type = option('ivm', 'exptime', default='ivm') pixel_scale_ratio = float(default=1.0) # Ratio of input to output pixel scale single = boolean(default=False) blendheaders = boolean(default=True) @@ -68,8 +68,8 @@ def process(self, input): self.log.info('Drizpars reference file: {}'.format(ref_filename)) kwargs = self.get_drizpars(ref_filename, input_models) else: - # Deal with NIRSpec which currently has no default drizpars reffile - self.log.info("No NIRSpec DIRZPARS reffile") + # If there is no drizpars reffile + self.log.info("No DIRZPARS reffile") kwargs = self._set_spec_defaults() kwargs['allowed_memory'] = self.allowed_memory @@ -113,7 +113,6 @@ def get_drizpars(self, ref_filename, input_models): pixfrac = float(default=None) kernel = string(default=None) fillval = string(default=None) - weight_type = option('exptime', default=None) Once the defaults are set from the reference file, if the user has used a resample.cfg file or run ResampleStep using command line args, @@ -148,15 +147,13 @@ def get_drizpars(self, ref_filename, input_models): self.log.error("No row found in %s matching input data.", ref_filename) raise ValueError - # Define the keys to pull from drizpars reffile table. Note the - # step param 'weight_type' is 'wht_type' in the FITS binary table. + # Define the keys to pull from drizpars reffile table. # All values should be None unless the user set them on the command # line or in the call to the step drizpars = dict( pixfrac=self.pixfrac, kernel=self.kernel, fillval=self.fillval, - wht_type=self.weight_type, pscale_ratio=self.pixel_scale_ratio, ) @@ -177,9 +174,6 @@ def get_drizpars(self, ref_filename, input_models): all_drizpars = {**reffile_drizpars, **user_drizpars} - # Convert the 'wht_type' key to a 'weight_type' key - all_drizpars['weight_type'] = all_drizpars.pop('wht_type') - kwargs = dict( good_bits=GOOD_BITS, single=self.single, @@ -188,12 +182,6 @@ def get_drizpars(self, ref_filename, input_models): kwargs.update(all_drizpars) - if 'wht_type' in kwargs: - raise DeprecationWarning('`wht_type` config keyword has changed ' + - 'to `weight_type`; ' + - 'please update calls to ResampleStep and resample.cfg files') - kwargs.pop('wht_type') - for k,v in kwargs.items(): self.log.debug(' {}={}'.format(k, v)) @@ -220,7 +208,7 @@ def _set_spec_defaults(self): if kwargs['fillval'] is None: kwargs['fillval'] = 'INDEF' if kwargs['weight_type'] is None: - kwargs['weight_type'] = 'exptime' + kwargs['weight_type'] = 'ivm' kwargs['pscale_ratio'] = self.pixel_scale_ratio kwargs.pop('pixel_scale_ratio') diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py index 4795d5bfd5..6559c680db 100644 --- a/jwst/resample/resample_utils.py +++ b/jwst/resample/resample_utils.py @@ -6,7 +6,7 @@ from astropy.modeling import Model from astropy.modeling.models import Mapping from astropy import units as u -from gwcs import WCS, wcstools +import gwcs from jwst.assign_wcs.util import wcs_from_footprints, wcs_bbox_from_shape from jwst.datamodels.dqflags import interpret_bit_flags @@ -72,7 +72,7 @@ def calc_gwcs_pixmap(in_wcs, out_wcs, shape=None): bb = in_wcs.bounding_box log.debug("Bounding box from WCS: {}".format(in_wcs.bounding_box)) - grid = wcstools.grid_from_bounding_box(bb) + grid = gwcs.wcstools.grid_from_bounding_box(bb) pixmap = np.dstack(reproject(in_wcs, out_wcs)(grid[0], grid[1])) return pixmap @@ -99,7 +99,7 @@ def reproject(wcs1, wcs2): if isinstance(wcs1, fitswcs.WCS): forward_transform = wcs1.all_pix2world - elif isinstance(wcs1, WCS): + elif isinstance(wcs1, gwcs.WCS): forward_transform = wcs1.forward_transform elif issubclass(wcs1, Model): forward_transform = wcs1 @@ -109,7 +109,7 @@ def reproject(wcs1, wcs2): if isinstance(wcs2, fitswcs.WCS): backward_transform = wcs2.all_world2pix - elif isinstance(wcs2, WCS): + elif isinstance(wcs2, gwcs.WCS): if not is_sky_like(wcs1.output_frame): # nirspec lamps: simplify backward transformation by omitting the msa_x (it's constant) # and just using the wavelength lookup table [1] and linear msa_y transformation [2] @@ -140,26 +140,35 @@ def _reproject(x, y): def build_driz_weight(model, weight_type=None, good_bits=None): - """ Create input weighting image + """Create a weight map for use by drizzle """ dqmask = build_mask(model.dq, good_bits) - exptime = model.meta.exposure.exposure_time - - if weight_type == 'error': - err_model = np.nan_to_num(model.err) - inwht = (exptime / err_model)**2 * dqmask - log.debug("DEBUG weight mask: {} {}".format(type(inwht), np.sum(inwht))) - # elif weight_type == 'ivm': - # _inwht = img.buildIVMmask(chip._chip,dqarr,pix_ratio) + + if weight_type == 'ivm': + if model.hasattr("var_rnoise"): + with np.errstate(divide="ignore", invalid="ignore"): + inv_variance = model.var_rnoise**-1 + inv_variance[~np.isfinite(inv_variance)] = 1 + else: + warnings.warn("var_rnoise array not available. Setting drizzle weight map to 1", + RuntimeWarning) + inv_variance = 1.0 + result = inv_variance * dqmask elif weight_type == 'exptime': - inwht = exptime * dqmask + exptime = model.meta.exposure.exposure_time + result = exptime * dqmask else: - inwht = np.ones(model.data.shape, dtype=model.data.dtype) - return inwht + warnings.warn("weight_type set to None. Setting drizzle weight map to 1", + RuntimeWarning) + result = np.ones(model.data.shape, dtype=model.data.dtype) + + return result.astype(np.float32) def build_mask(dqarr, bitvalue): - """ Builds a bit-mask from an input DQ array and a bitvalue flag + """Build a bit mask from an input DQ array and a bitvalue flag + + In the returned bit mask, 1 is good, 0 is bad """ bitvalue = interpret_bit_flags(bitvalue) @@ -167,6 +176,7 @@ def build_mask(dqarr, bitvalue): return (np.ones(dqarr.shape, dtype=np.uint8)) return np.logical_not(np.bitwise_and(dqarr, ~bitvalue)).astype(np.uint8) + def is_sky_like(frame): # Differentiate between sky-like and cartesian frames return u.Unit("deg") in frame.unit or u.Unit("arcsec") in frame.unit diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py index 8a4cd2633a..1a514ed25e 100644 --- a/jwst/resample/tests/test_resample_step.py +++ b/jwst/resample/tests/test_resample_step.py @@ -3,14 +3,16 @@ import numpy as np import pytest -from jwst.datamodels import ImageModel +from jwst.datamodels import ImageModel, ModelContainer from jwst.assign_wcs import AssignWcsStep from jwst.extract_2d import Extract2dStep from jwst.resample import ResampleSpecStep, ResampleStep @pytest.fixture def nirspec_rate(): - shape = (2048, 2048) + ysize = 2048 + xsize = 2048 + shape = (ysize, xsize) im = ImageModel(shape) im.meta.wcsinfo = { 'dec_ref': 40, @@ -65,8 +67,10 @@ def nirspec_rate(): def miri_rate(): xsize = 72 ysize = 416 - im = ImageModel((ysize, xsize)) - im.data += 5. + shape = (ysize, xsize) + im = ImageModel(shape) + im.data += 5 + im.var_rnoise = np.random.random(shape) im.meta.wcsinfo = { 'dec_ref': 40, 'ra_ref': 100, @@ -116,8 +120,8 @@ def miri_rate(): def nircam_rate(): xsize = 204 ysize = 204 - im = ImageModel((ysize, xsize)) - im.data += 5. + shape = (ysize, xsize) + im = ImageModel(shape) im.meta.wcsinfo = { 'ctype1': 'RA---TAN', 'ctype2': 'DEC--TAN', @@ -222,6 +226,7 @@ def test_pixel_scale_ratio_spec(miri_rate, ratio): @pytest.mark.parametrize("ratio", [0.5, 0.7, 1.0]) def test_pixel_scale_ratio_imaging(nircam_rate, ratio): im = AssignWcsStep.call(nircam_rate, sip_approx=False) + im.data += 5 result1 = ResampleStep.call(im) result2 = ResampleStep.call(im, pixel_scale_ratio=ratio) @@ -237,6 +242,37 @@ def test_pixel_scale_ratio_imaging(nircam_rate, ratio): assert_allclose(area1 * ratio**2, area2, rtol=1e-6) +def test_weight_type(nircam_rate, _jail): + """Check that weight_type of exptime and ivm work""" + im1 = AssignWcsStep.call(nircam_rate, sip_approx=False) + im2 = im1.copy() + im3 = im1.copy() + im1.data += 10 + im2.data += 5 + im3.data += 5 + im1.var_rnoise += (1 / 10) + im2.var_rnoise += (1 / 5) + im3.var_rnoise += (1 / 5) + im2.meta.observation.sequence_id = "2" + im3.meta.observation.sequence_id = "3" + + c = ModelContainer([im1, im2, im3]) + assert len(c.group_names) == 3 + + result1 = ResampleStep.call(c, weight_type="ivm", blendheaders=False, save_results=True) + + # assert_allclose(result1.data, result2.data) + # assert_allclose(result1.wht, result2.wht) + assert_allclose(result1.data[100:105, 100:105], 7.5, rtol=1e-2) + assert_allclose(result1.wht[100:105, 100:105], 20, rtol=1e-2) + + result2 = ResampleStep.call(c, weight_type="exptime", blendheaders=False) + + assert_allclose(result2.data[100:105, 100:105], 7.5, rtol=1e-2) + assert_allclose(result2.wht[100:105, 100:105], 20, rtol=1e-2) + + + def test_sip_coeffs_do_not_propagate(nircam_rate): im = AssignWcsStep.call(nircam_rate, sip_degree=2) diff --git a/jwst/resample/tests/test_utils.py b/jwst/resample/tests/test_utils.py index c6c4955962..1951c88909 100644 --- a/jwst/resample/tests/test_utils.py +++ b/jwst/resample/tests/test_utils.py @@ -1,11 +1,15 @@ """Test various utility functions""" -import pytest - +from numpy.testing import assert_array_equal import numpy as np +import pytest -from jwst.datamodels import SlitModel +from jwst.datamodels import SlitModel, ImageModel, dqflags from jwst.resample.resample_spec import find_dispersion_axis -from jwst.resample.resample_utils import build_mask +from jwst.resample.resample_utils import build_mask, build_driz_weight + + +DO_NOT_USE = dqflags.pixel["DO_NOT_USE"] +GOOD = dqflags.pixel["GOOD"] DQ = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8]) @@ -40,7 +44,32 @@ def test_build_mask(dq, bitvalues, expected): Expected mask array """ result = build_mask(dq, bitvalues) - assert np.array_equal(result, expected) + assert_array_equal(result, expected) + + +@pytest.mark.parametrize("weight_type", ["ivm", "exptime"]) +def test_build_driz_weight(weight_type): + """Check that correct weight map is returned of different weight types""" + model = ImageModel((10, 10)) + model.dq[0] = DO_NOT_USE + model.meta.exposure.exposure_time = 10.0 + model.var_rnoise += 0.1 + + weight_map = build_driz_weight(model, weight_type=weight_type, good_bits="GOOD") + assert_array_equal(weight_map[0], 0) + assert_array_equal(weight_map[1:], 10.0) + assert weight_map.dtype == np.float32 + + +@pytest.mark.parametrize("weight_type", ["ivm", None]) +def test_build_driz_weight_zeros(weight_type): + """Check that zero or not finite weight maps get set to 1""" + model = ImageModel((10, 10)) + + with pytest.warns(RuntimeWarning): + weight_map = build_driz_weight(model, weight_type=weight_type) + + assert_array_equal(weight_map, 1) def test_find_dispersion_axis():