From 95385d0ab45989b80f89cbd78d15ef410dab0cf1 Mon Sep 17 00:00:00 2001 From: James Davies Date: Thu, 11 Feb 2021 08:02:18 -0500 Subject: [PATCH 01/10] Make ivm default weight_type for resample --- jwst/datamodels/schemas/core.schema.yaml | 2 +- jwst/resample/gwcs_drizzle.py | 29 ----------------------- jwst/resample/resample.py | 5 ++-- jwst/resample/resample_step.py | 17 +++---------- jwst/resample/resample_utils.py | 9 +++---- jwst/resample/tests/test_resample_step.py | 27 +++++++++++++++++---- 6 files changed, 31 insertions(+), 58 deletions(-) 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/resample/gwcs_drizzle.py b/jwst/resample/gwcs_drizzle.py index 34e2461983..70ee37cf7b 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 @@ -199,34 +198,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..6ed1c3611f 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 diff --git a/jwst/resample/resample_step.py b/jwst/resample/resample_step.py index a5ac0b655b..2b93b90116 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('exptime', 'ivm', default='ivm') pixel_scale_ratio = float(default=1.0) # Ratio of input to output pixel scale single = boolean(default=False) blendheaders = boolean(default=True) @@ -148,15 +148,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 +175,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 +183,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 +209,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..a71c0a7ceb 100644 --- a/jwst/resample/resample_utils.py +++ b/jwst/resample/resample_utils.py @@ -145,16 +145,13 @@ def build_driz_weight(model, weight_type=None, good_bits=None): 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': + inwht = np.nan_to_num(1 / model.var_rnoise * dqmask) elif weight_type == 'exptime': inwht = exptime * dqmask else: inwht = np.ones(model.data.shape, dtype=model.data.dtype) + return inwht diff --git a/jwst/resample/tests/test_resample_step.py b/jwst/resample/tests/test_resample_step.py index 8a4cd2633a..c90775447e 100644 --- a/jwst/resample/tests/test_resample_step.py +++ b/jwst/resample/tests/test_resample_step.py @@ -10,8 +10,11 @@ @pytest.fixture def nirspec_rate(): - shape = (2048, 2048) + ysize = 2048 + xsize = 2048 + shape = (ysize, xsize) im = ImageModel(shape) + im.var_rnoise = np.random.random(shape) im.meta.wcsinfo = { 'dec_ref': 40, 'ra_ref': 100, @@ -65,8 +68,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 +121,10 @@ def miri_rate(): def nircam_rate(): xsize = 204 ysize = 204 - 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 = { 'ctype1': 'RA---TAN', 'ctype2': 'DEC--TAN', @@ -237,6 +244,16 @@ def test_pixel_scale_ratio_imaging(nircam_rate, ratio): assert_allclose(area1 * ratio**2, area2, rtol=1e-6) +def test_weight_type(nircam_rate): + """Check that weight_type of exptime and ivm work""" + im = AssignWcsStep.call(nircam_rate, sip_approx=False) + result1 = ResampleStep.call(im, weight_type="ivm") + result2 = ResampleStep.call(im, weight_type="exptime") + + assert_allclose(result1.data, result2.data) + assert_allclose(result1.wht, result2.wht) + + def test_sip_coeffs_do_not_propagate(nircam_rate): im = AssignWcsStep.call(nircam_rate, sip_degree=2) From 7f8bca2d87f88a26fc00de1aa4448270192d83aa Mon Sep 17 00:00:00 2001 From: James Davies Date: Sat, 13 Feb 2021 10:37:12 -0500 Subject: [PATCH 02/10] Remove CRDS printing from Jenkinsfiles --- JenkinsfileRT | 1 - JenkinsfileRT_dev | 1 - 2 files changed, 2 deletions(-) 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 = [ From 5ce0ffae9bce37162b4f714a3bd6a061172b3e30 Mon Sep 17 00:00:00 2001 From: James Davies Date: Wed, 17 Feb 2021 16:49:04 -0500 Subject: [PATCH 03/10] =?UTF-8?q?Clean=20up=20clean=20up=20=F0=9F=8E=B5?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- jwst/resample/gwcs_drizzle.py | 27 ++++++---------- jwst/resample/resample.py | 3 +- jwst/resample/resample_spec.py | 3 +- jwst/resample/resample_step.py | 5 ++- jwst/resample/resample_utils.py | 25 +++++++++------ jwst/resample/tests/test_resample_step.py | 39 +++++++++++++++++------ 6 files changed, 58 insertions(+), 44 deletions(-) diff --git a/jwst/resample/gwcs_drizzle.py b/jwst/resample/gwcs_drizzle.py index 70ee37cf7b..25dbc2f62f 100644 --- a/jwst/resample/gwcs_drizzle.py +++ b/jwst/resample/gwcs_drizzle.py @@ -12,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,9 +32,9 @@ def __init__(self, product, outwcs=None, single=False, 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. + 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 @@ -63,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 @@ -89,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)) @@ -101,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": diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py index 6ed1c3611f..44c4a12116 100644 --- a/jwst/resample/resample.py +++ b/jwst/resample/resample.py @@ -163,8 +163,7 @@ def do_drizzle(self): inwht = resample_utils.build_driz_weight(img, 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) + driz.add_image(img.data, img.meta.wcs, inwht=inwht) # 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..a03e01a9f6 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 ---------- @@ -414,7 +414,6 @@ def do_drizzle(self, xmin=0, xmax=0, ymin=0, ymax=0, **pars): in_wcs = img.meta.wcs driz.add_image(img.data, in_wcs, inwht=inwht, - expin=img.meta.exposure.exposure_time, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) # Update some basic exposure time values based on all the inputs diff --git a/jwst/resample/resample_step.py b/jwst/resample/resample_step.py index 2b93b90116..0a3b37a91a 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', 'ivm', default='ivm') + 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,7 +68,7 @@ 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 + # Deal with FGS which currently has no default drizpars reffile self.log.info("No NIRSpec DIRZPARS reffile") kwargs = self._set_spec_defaults() @@ -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, diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py index a71c0a7ceb..02c1821764 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,23 +140,27 @@ 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 == 'ivm': - inwht = np.nan_to_num(1 / model.var_rnoise * dqmask) + with np.errstate(divide="ignore", invalid="ignore"): + result = 1 / model.var_rnoise * dqmask + result[~np.isfinite(result)] = 0 elif weight_type == 'exptime': - inwht = exptime * dqmask + result = exptime * dqmask else: - inwht = np.ones(model.data.shape, dtype=model.data.dtype) + result = np.ones(model.data.shape, dtype=model.data.dtype) - return inwht + return result 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) @@ -164,6 +168,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 c90775447e..1a514ed25e 100644 --- a/jwst/resample/tests/test_resample_step.py +++ b/jwst/resample/tests/test_resample_step.py @@ -3,7 +3,7 @@ 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 @@ -14,7 +14,6 @@ def nirspec_rate(): xsize = 2048 shape = (ysize, xsize) im = ImageModel(shape) - im.var_rnoise = np.random.random(shape) im.meta.wcsinfo = { 'dec_ref': 40, 'ra_ref': 100, @@ -123,8 +122,6 @@ def nircam_rate(): ysize = 204 shape = (ysize, xsize) im = ImageModel(shape) - im.data += 5 - im.var_rnoise = np.random.random(shape) im.meta.wcsinfo = { 'ctype1': 'RA---TAN', 'ctype2': 'DEC--TAN', @@ -229,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) @@ -244,14 +242,35 @@ def test_pixel_scale_ratio_imaging(nircam_rate, ratio): assert_allclose(area1 * ratio**2, area2, rtol=1e-6) -def test_weight_type(nircam_rate): +def test_weight_type(nircam_rate, _jail): """Check that weight_type of exptime and ivm work""" - im = AssignWcsStep.call(nircam_rate, sip_approx=False) - result1 = ResampleStep.call(im, weight_type="ivm") - result2 = ResampleStep.call(im, weight_type="exptime") + 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) - assert_allclose(result1.data, result2.data) - assert_allclose(result1.wht, result2.wht) def test_sip_coeffs_do_not_propagate(nircam_rate): From 75ff5e21de881482e6741cb61437b730da6f56fd Mon Sep 17 00:00:00 2001 From: James Davies Date: Wed, 17 Feb 2021 18:44:14 -0500 Subject: [PATCH 04/10] Update CHANGES.rst --- CHANGES.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 9ebe78b840..53612c307b 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -146,6 +146,12 @@ ramp_fitting - Update documentation to define optimal weighting algorithm [#5682] +resample +-------- + +- Make inverse variance ``weight_type="ivm"`` the default weighting scheme for + multiple exposures resampled into a single output. [#5738] + source_catalog -------------- From 5e41405ed44bae090a9f49f2468e910619893425 Mon Sep 17 00:00:00 2001 From: James Davies Date: Thu, 18 Feb 2021 23:46:33 -0500 Subject: [PATCH 05/10] Address review comments --- jwst/resample/gwcs_drizzle.py | 4 ++-- jwst/resample/resample.py | 3 ++- jwst/resample/resample_spec.py | 1 + jwst/resample/resample_step.py | 4 ++-- 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/jwst/resample/gwcs_drizzle.py b/jwst/resample/gwcs_drizzle.py index 25dbc2f62f..a8b8494339 100644 --- a/jwst/resample/gwcs_drizzle.py +++ b/jwst/resample/gwcs_drizzle.py @@ -31,8 +31,8 @@ def __init__(self, product, outwcs=None, single=False, wt_scl=None, 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 or `expsq` which scales + 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. diff --git a/jwst/resample/resample.py b/jwst/resample/resample.py index 44c4a12116..74fb5bf082 100644 --- a/jwst/resample/resample.py +++ b/jwst/resample/resample.py @@ -163,7 +163,8 @@ def do_drizzle(self): inwht = resample_utils.build_driz_weight(img, weight_type=self.weight_type, good_bits=self.good_bits) - driz.add_image(img.data, img.meta.wcs, inwht=inwht) + driz.add_image(img.data, img.meta.wcs, inwht=inwht, + 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 a03e01a9f6..6a84fc2a05 100644 --- a/jwst/resample/resample_spec.py +++ b/jwst/resample/resample_spec.py @@ -414,6 +414,7 @@ def do_drizzle(self, xmin=0, xmax=0, ymin=0, ymax=0, **pars): in_wcs = img.meta.wcs driz.add_image(img.data, in_wcs, inwht=inwht, + expin=img.meta.exposure.exposure_time, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax) # Update some basic exposure time values based on all the inputs diff --git a/jwst/resample/resample_step.py b/jwst/resample/resample_step.py index 0a3b37a91a..1692dc4418 100755 --- a/jwst/resample/resample_step.py +++ b/jwst/resample/resample_step.py @@ -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 FGS 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 From 6504363b9875e1aba2b5e82903848cc0ee251f75 Mon Sep 17 00:00:00 2001 From: James Davies Date: Fri, 19 Feb 2021 14:49:57 -0500 Subject: [PATCH 06/10] Make outlier_detection use ivm --- jwst/outlier_detection/outlier_detection.py | 4 ++-- jwst/outlier_detection/outlier_detection_scaled.py | 2 +- jwst/outlier_detection/outlier_detection_scaled_step.py | 2 +- jwst/outlier_detection/outlier_detection_spec.py | 2 +- jwst/outlier_detection/outlier_detection_stack_step.py | 2 +- jwst/outlier_detection/outlier_detection_step.py | 2 +- 6 files changed, 7 insertions(+), 7 deletions(-) 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') From a04b82c5edf0a48fc6adb7a07fd30969dcc22853 Mon Sep 17 00:00:00 2001 From: James Davies Date: Mon, 22 Feb 2021 12:50:42 -0500 Subject: [PATCH 07/10] Handle the case where there's no var_rnoise --- jwst/resample/resample_utils.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py index 02c1821764..2d0204ea39 100644 --- a/jwst/resample/resample_utils.py +++ b/jwst/resample/resample_utils.py @@ -143,13 +143,17 @@ def build_driz_weight(model, weight_type=None, good_bits=None): """Create a weight map for use by drizzle """ dqmask = build_mask(model.dq, good_bits) - exptime = model.meta.exposure.exposure_time if weight_type == 'ivm': - with np.errstate(divide="ignore", invalid="ignore"): - result = 1 / model.var_rnoise * dqmask - result[~np.isfinite(result)] = 0 + 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: + inv_variance = 1.0 + result = inv_variance * dqmask elif weight_type == 'exptime': + exptime = model.meta.exposure.exposure_time result = exptime * dqmask else: result = np.ones(model.data.shape, dtype=model.data.dtype) From 1a91c8d598f15706f03ea76de3b9684558740b3c Mon Sep 17 00:00:00 2001 From: James Davies Date: Wed, 24 Feb 2021 23:33:49 -0500 Subject: [PATCH 08/10] Move entry to unreleased part of CHANGES.rst --- CHANGES.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 53612c307b..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) ================== @@ -146,12 +152,6 @@ ramp_fitting - Update documentation to define optimal weighting algorithm [#5682] -resample --------- - -- Make inverse variance ``weight_type="ivm"`` the default weighting scheme for - multiple exposures resampled into a single output. [#5738] - source_catalog -------------- From ad9d56a407fc52ec0c5c6148f2203d3f7b147edd Mon Sep 17 00:00:00 2001 From: James Davies Date: Wed, 24 Feb 2021 23:38:46 -0500 Subject: [PATCH 09/10] Add tests for build_driz_weight --- jwst/resample/tests/test_utils.py | 40 ++++++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/jwst/resample/tests/test_utils.py b/jwst/resample/tests/test_utils.py index c6c4955962..c5ae38c197 100644 --- a/jwst/resample/tests/test_utils.py +++ b/jwst/resample/tests/test_utils.py @@ -3,9 +3,13 @@ import numpy as np -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,37 @@ def test_build_mask(dq, bitvalues, expected): Expected mask array """ result = build_mask(dq, bitvalues) - assert np.array_equal(result, expected) + np.testing.assert_array_equal(result, expected) + + +@pytest.mark.parametrize("weight_type, expected_value", [ + ("ivm", 0.1), + ("exptime", 10), + (None, 1.0), + ] +) +def test_build_driz_weight(weight_type, expected_value): + """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 += 10.0 + + weight_map = build_driz_weight(model, weight_type=weight_type, good_bits="GOOD") + np.testing.assert_array_equal(weight_map[0], 0) + np.testing.assert_array_equal(weight_map[1:], expected_value) + assert weight_map.dtype == np.float32 + + +@pytest.mark.parametrize("weight_type", ["ivm", "exptime"]) +def test_build_driz_weight_zeros(weight_type): + """Check that zero or not finite weight maps raise an error""" + model = ImageModel((10, 10)) + model.var_rnoise += 0 + model.meta.exposure.exposure_time = 0.0 + + with pytest.raises(ValueError): + build_driz_weight(model, weight_type=weight_type) def test_find_dispersion_axis(): From c283048043ee351b4be8452d88e19e88de432025 Mon Sep 17 00:00:00 2001 From: James Davies Date: Thu, 25 Feb 2021 00:29:57 -0500 Subject: [PATCH 10/10] Issue warning when var_rnoise doesn't exist --- jwst/resample/resample_utils.py | 6 +++++- jwst/resample/tests/test_utils.py | 33 +++++++++++++------------------ 2 files changed, 19 insertions(+), 20 deletions(-) diff --git a/jwst/resample/resample_utils.py b/jwst/resample/resample_utils.py index 2d0204ea39..6559c680db 100644 --- a/jwst/resample/resample_utils.py +++ b/jwst/resample/resample_utils.py @@ -150,15 +150,19 @@ def build_driz_weight(model, weight_type=None, good_bits=None): 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': exptime = model.meta.exposure.exposure_time result = exptime * dqmask else: + 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 + return result.astype(np.float32) def build_mask(dqarr, bitvalue): diff --git a/jwst/resample/tests/test_utils.py b/jwst/resample/tests/test_utils.py index c5ae38c197..1951c88909 100644 --- a/jwst/resample/tests/test_utils.py +++ b/jwst/resample/tests/test_utils.py @@ -1,7 +1,7 @@ """Test various utility functions""" -import pytest - +from numpy.testing import assert_array_equal import numpy as np +import pytest from jwst.datamodels import SlitModel, ImageModel, dqflags from jwst.resample.resample_spec import find_dispersion_axis @@ -44,37 +44,32 @@ def test_build_mask(dq, bitvalues, expected): Expected mask array """ result = build_mask(dq, bitvalues) - np.testing.assert_array_equal(result, expected) + assert_array_equal(result, expected) -@pytest.mark.parametrize("weight_type, expected_value", [ - ("ivm", 0.1), - ("exptime", 10), - (None, 1.0), - ] -) -def test_build_driz_weight(weight_type, expected_value): +@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 += 10.0 + model.var_rnoise += 0.1 weight_map = build_driz_weight(model, weight_type=weight_type, good_bits="GOOD") - np.testing.assert_array_equal(weight_map[0], 0) - np.testing.assert_array_equal(weight_map[1:], expected_value) + 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", "exptime"]) +@pytest.mark.parametrize("weight_type", ["ivm", None]) def test_build_driz_weight_zeros(weight_type): - """Check that zero or not finite weight maps raise an error""" + """Check that zero or not finite weight maps get set to 1""" model = ImageModel((10, 10)) - model.var_rnoise += 0 - model.meta.exposure.exposure_time = 0.0 - with pytest.raises(ValueError): - build_driz_weight(model, weight_type=weight_type) + 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():