From a4b23a99cde10f1469b59b09e59440dc87089e58 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Thu, 29 Feb 2024 13:02:54 -0500 Subject: [PATCH 01/15] Update wavecorr to add a zeropoint correction to the wcs pipeline Refactoring wavecorr changes for tests minor refactoring Add wavecorr zeropoint correction to slit wcs --- jwst/wavecorr/wavecorr.py | 115 +++++++++++++++++++++++++++----------- 1 file changed, 81 insertions(+), 34 deletions(-) diff --git a/jwst/wavecorr/wavecorr.py b/jwst/wavecorr/wavecorr.py index 0174eb1f59..7eb0c3715a 100644 --- a/jwst/wavecorr/wavecorr.py +++ b/jwst/wavecorr/wavecorr.py @@ -30,6 +30,8 @@ import logging import numpy as np from gwcs import wcstools +from astropy.modeling import tabular +from astropy.modeling.mappings import Identity from stdatamodels.jwst import datamodels from stdatamodels.jwst.transforms import models as trmodels @@ -127,43 +129,50 @@ def apply_zero_point_correction(slit, reffile): # For the MSA the aperture name is "MOS" aperture_name = "MOS" - lam = slit.wavelength.copy() + lam = slit.wavelength.copy() * 1e-6 dispersion = compute_dispersion(slit.meta.wcs) - corr, dq_lam = compute_zero_point_correction(lam, reffile, source_xpos, - aperture_name, dispersion) - # TODO: set a DQ flag to a TBD value for pixels where dq_lam == 0. - # The only purpose of dq_lam is to set that flag. - # Wavelength is in um, the correction is computed in meters. - slit.wavelength = slit.wavelength - corr * 10 ** 6 - - -def compute_zero_point_correction(lam, freference, source_xpos, - aperture_name, dispersion): - """ Compute the NIRSpec wavelength zero-point correction. + wave2wavecorr = calculate_wavelength_correction_transform(lam, dispersion, + reffile, source_xpos, + aperture_name) + + # Insert the new transform into the slit wcs object + wave2wavecorr = Identity(2) & wave2wavecorr + slit_wcs.insert_transform('slit_frame', transform = wave2wavecorr, after=False) + + # Update the stored wavelengths for the slit + slit.wavelength = compute_wavelength(slit_wcs) + + +def calculate_wavelength_correction_transform(lam, dispersion, freference, + source_xpos, aperture_name): + """ Generate a WCS transform for the NIRSpec wavelength zero-point correction + and add it to the WCS for each slit. Parameters ---------- + slit_wcs : ~gwcs.wcs.WCS` + The WCS for this slit. lam : ndarray - Wavelength array. + Wavelength array [in m]. + dispersion : ndarray + The pixel dispersion [in m]. freference : str ``wavecorr`` reference file name. source_xpos : float X position of the source as a fraction of the slit size. aperture_name : str Aperture name. - dispersion : ndarray - The pixel dispersion [in m]. - + Returns ------- - lambda_corr : ndarray - Wavelength correction. - lam : ndarray - Interpolated wavelengths. Extrapolated values are reset to 0. - This is returned so that the DQ array can be updated with a flag - which indicates that no zero-point correction was done. + model : `~astropy.modeling.tabular.Tabular1D + A model which takes wavelength inputs and returns zeropoint + corrected wavelengths. """ + + + # Open the zerpoint reference model with datamodels.WaveCorrModel(freference) as wavecorr: for ap in wavecorr.apertures: if ap.aperture_name == aperture_name: @@ -175,19 +184,31 @@ def compute_zero_point_correction(lam, freference, source_xpos, break else: log.info(f'No wavelength zero-point correction found for slit {aperture_name}') - - deltax = source_xpos - lam = lam.copy() - lam_no_nans = lam[~np.isnan(lam)] + + # Set lookup table to extrapolate at bounds to recover wavelengths beyond model bounds + # particulary for the red and blue ends of prism observations offset_model.bounds_error = False - correction = offset_model(lam_no_nans * 10 ** -6, [deltax] * lam_no_nans.size) - lam[~np.isnan(lam)] = correction - - # The correction for pixels outside the slit and wavelengths - # outside the wave_range is 0. - lam[np.isnan(lam)] = 0. - lambda_cor = dispersion * lam - return lambda_cor, lam + offset_model.fill_value = None + + # Average the wavelength and dispersion across 2D extracted slit and remove nans + # So that we have a 1D wavelength array for building a 1D lookup table wcs transform + lam_mean = np.nanmean(lam, axis=0) + disp_mean = np.nanmean(dispersion, axis=0) + nan_lams = np.isnan(lam_mean) | np.isnan(disp_mean) + lam_mean = lam_mean[~nan_lams] + disp_mean = disp_mean[~nan_lams] + + # Calculate the corrected wavelengths + pixel_corrections = offset_model(lam_mean, source_xpos) + lam_corrected = lam_mean + (pixel_corrections * disp_mean) + + # Build a look up table to transform between corrected and uncorrected wavelengths + wave2wavecorr = tabular.Tabular1D(points=lam_mean, lookup_table=lam_corrected, + bounds_error=False, fill_value=None) + wave2wavecorr.inverse = tabular.Tabular1D(points=lam_corrected, lookup_table=lam_mean, + bounds_error=False, fill_value=None) + + return wave2wavecorr def compute_dispersion(wcs, xpix=None, ypix=None): @@ -218,6 +239,32 @@ def compute_dispersion(wcs, xpix=None, ypix=None): return (lamright - lamleft) * 10 ** -6 +def compute_wavelength(wcs, xpix=None, ypix=None): + """ Compute the pixel wavelength. + + Parameters + ---------- + wcs : `~gwcs.wcs.WCS` + The WCS object for this slit. + xpix : ndarray, float, optional + ypix : ndarray, float, optional + Compute the wavelength at the x, y pixels. + If not provided the dispersion is computed on a + grid based on ``wcs.bounding_box``. + + Returns + ------- + wavelength : ndarray + The wavelength [in microns]. + + """ + if xpix is None or ypix is None: + xpix, ypix = wcstools.grid_from_bounding_box(wcs.bounding_box, step=(1, 1)) + + _, _, lam = wcs(xpix, ypix) + return lam + + def _is_point_source(slit, exp_type): """ Determine if a source is a point source. From 788ea61eb3dbed7e5284ff265de8e07717e71595 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Fri, 8 Mar 2024 11:33:46 -0500 Subject: [PATCH 02/15] updating wavecorr tests --- jwst/wavecorr/tests/test_wavecorr.py | 83 +++++++++++++++++++++------- 1 file changed, 64 insertions(+), 19 deletions(-) diff --git a/jwst/wavecorr/tests/test_wavecorr.py b/jwst/wavecorr/tests/test_wavecorr.py index e4a06cb625..2a358fae52 100644 --- a/jwst/wavecorr/tests/test_wavecorr.py +++ b/jwst/wavecorr/tests/test_wavecorr.py @@ -27,32 +27,58 @@ def test_wavecorr(): im = datamodels.ImageModel(hdul) im_wcs = AssignWcsStep.call(im) im_ex2d = Extract2dStep.call(im_wcs) - im_ex2d.slits[0].meta.wcs.bounding_box = ((-.5, 1432.5), (-.5, 37.5)) + bbox = ((-.5, 1432.5), (-.5, 37.5)) + im_ex2d.slits[0].meta.wcs.bounding_box = bbox + x, y = wcstools.grid_from_bounding_box(bbox) + ra, dec, lam_before = im_ex2d.slits[0].meta.wcs(x, y) + im_ex2d.slits[0].wavelength = lam_before im_src = SourceTypeStep.call(im_ex2d) + # the mock msa source is an extended source, change to point for testing + im_src.slits[0].source_type = 'POINT' im_wave = WavecorrStep.call(im_src) # test dispersion is of the correct order # there's one slit only slit = im_src.slits[0] - x, y = wcstools.grid_from_bounding_box(slit.meta.wcs.bounding_box) dispersion = wavecorr.compute_dispersion(slit.meta.wcs, x, y) assert_allclose(dispersion[~np.isnan(dispersion)], 1e-9, atol=1e-10) + + # test that the wavelength is on the order of microns + wavelength = wavecorr.compute_wavelength(slit.meta.wcs, x, y) + assert_allclose(np.nanmean(wavelength), 2.5, atol=0.1) + + # Check that the stored wavelengths were corrected + abs_wave_correction = np.abs(im_src.slits[0].wavelength - im_wave.slits[0].wavelength) + assert_allclose(np.nanmean(abs_wave_correction), 0.00046, atol=0.0001) + + # Check that the slit wcs has been updated to provide corrected wavelengths + corrected_wavelength = wavecorr.compute_wavelength(im_wave.slits[0].meta.wcs, x, y) + assert_allclose(im_wave.slits[0].wavelength, corrected_wavelength) + + # test the roundtripping on the wavelength correction transform + ref_name = im_wave.meta.ref_file.wavecorr.name + freference = datamodels.WaveCorrModel(WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) - # the difference in wavelength should be of the order of e-10 in um - assert_allclose(im_src.slits[0].wavelength - im_wave.slits[0].wavelength, 1e-10) - + lam_uncorr = lam_before * 1e-6 + wave2wavecorr = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, + freference, slit.source_xpos, 'MOS') + lam_corr = wave2wavecorr(lam_uncorr) + assert_allclose(lam_uncorr, wave2wavecorr.inverse(lam_corr)) + # test on both sides of the shutter source_xpos1 = -.2 source_xpos2 = .2 - ra, dec, lam = slit.meta.wcs(x, y) - ref_name = im_wave.meta.ref_file.wavecorr.name - freference = datamodels.WaveCorrModel(WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) - zero_point1 = wavecorr.compute_zero_point_correction(lam, freference, source_xpos1, 'MOS', dispersion) - zero_point2 = wavecorr.compute_zero_point_correction(lam, freference, source_xpos2, 'MOS', dispersion) - diff_correction = np.abs(zero_point1[1] - zero_point2[1]) - non_zero = diff_correction[diff_correction != 0] - assert_allclose(np.nanmean(non_zero), 0.75, atol=0.01) + wave_transform1 = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, + freference, source_xpos1, 'MOS') + wave_transform2 = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, + freference, source_xpos2, 'MOS') + + zero_point1 = wave_transform1(lam_uncorr) + zero_point2 = wave_transform2(lam_uncorr) + + diff_correction = np.abs(zero_point1 - zero_point2) + assert_allclose(np.nanmean(diff_correction), 8.0e-10, atol=0.1e-10) def test_ideal_to_v23_fs(): @@ -171,13 +197,32 @@ def test_wavecorr_fs(): dispersion = wavecorr.compute_dispersion(slit.meta.wcs, x, y) assert_allclose(dispersion[~np.isnan(dispersion)], 1e-8, atol=1.04e-8) + # Check that the slit wcs has been updated to provide corrected wavelengths + corrected_wavelength = wavecorr.compute_wavelength(slit.meta.wcs, x, y) + assert_allclose(slit.wavelength, corrected_wavelength) + + # test the roundtripping on the wavelength correction transform + ref_name = result.meta.ref_file.wavecorr.name + freference = datamodels.WaveCorrModel(WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) + + lam_uncorr = lam_before * 1e-6 + wave2wavecorr = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, + freference, slit.source_xpos, 'S200A1') + lam_corr = wave2wavecorr(lam_uncorr) + assert_allclose(lam_uncorr, wave2wavecorr.inverse(lam_corr)) + + # test on both sides of the slit center source_xpos1 = -.2 source_xpos2 = .2 - ref_name = result.meta.ref_file.wavecorr.name - freference = datamodels.WaveCorrModel(WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) - zero_point1 = wavecorr.compute_zero_point_correction(lam_before, freference, source_xpos1, 'S200A1', dispersion) - zero_point2 = wavecorr.compute_zero_point_correction(lam_before, freference, source_xpos2, 'S200A1', dispersion) - diff_correction = np.abs(zero_point1[1] - zero_point2[1]) - assert_allclose(np.nanmean(diff_correction), 0.45, atol=0.01) + wave_transform1 = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, + freference, source_xpos1, 'S200A1') + wave_transform2 = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, + freference, source_xpos2, 'S200A1') + + zero_point1 = wave_transform1(lam_uncorr) + zero_point2 = wave_transform2(lam_uncorr) + + diff_correction = np.abs(zero_point1 - zero_point2) + assert_allclose(np.nanmean(diff_correction), 6.3e-9, atol=0.1e-9) From 882eafc8ef620101cc9e5a014ab392b721c7c083 Mon Sep 17 00:00:00 2001 From: Melanie Clarke Date: Tue, 12 Mar 2024 16:02:45 -0400 Subject: [PATCH 03/15] Code style fixes --- jwst/wavecorr/tests/test_wavecorr.py | 27 ++++++++++++++------------- jwst/wavecorr/wavecorr.py | 20 ++++++++------------ 2 files changed, 22 insertions(+), 25 deletions(-) diff --git a/jwst/wavecorr/tests/test_wavecorr.py b/jwst/wavecorr/tests/test_wavecorr.py index 2a358fae52..047c856637 100644 --- a/jwst/wavecorr/tests/test_wavecorr.py +++ b/jwst/wavecorr/tests/test_wavecorr.py @@ -33,6 +33,7 @@ def test_wavecorr(): ra, dec, lam_before = im_ex2d.slits[0].meta.wcs(x, y) im_ex2d.slits[0].wavelength = lam_before im_src = SourceTypeStep.call(im_ex2d) + # the mock msa source is an extended source, change to point for testing im_src.slits[0].source_type = 'POINT' im_wave = WavecorrStep.call(im_src) @@ -55,13 +56,14 @@ def test_wavecorr(): corrected_wavelength = wavecorr.compute_wavelength(im_wave.slits[0].meta.wcs, x, y) assert_allclose(im_wave.slits[0].wavelength, corrected_wavelength) - # test the roundtripping on the wavelength correction transform + # test the round-tripping on the wavelength correction transform ref_name = im_wave.meta.ref_file.wavecorr.name - freference = datamodels.WaveCorrModel(WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) + freference = datamodels.WaveCorrModel( + WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) lam_uncorr = lam_before * 1e-6 - wave2wavecorr = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, - freference, slit.source_xpos, 'MOS') + wave2wavecorr = wavecorr.calculate_wavelength_correction_transform( + lam_uncorr, dispersion, freference, slit.source_xpos, 'MOS') lam_corr = wave2wavecorr(lam_uncorr) assert_allclose(lam_uncorr, wave2wavecorr.inverse(lam_corr)) @@ -69,10 +71,10 @@ def test_wavecorr(): source_xpos1 = -.2 source_xpos2 = .2 - wave_transform1 = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, - freference, source_xpos1, 'MOS') - wave_transform2 = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, - freference, source_xpos2, 'MOS') + wave_transform1 = wavecorr.calculate_wavelength_correction_transform( + lam_uncorr, dispersion, freference, source_xpos1, 'MOS') + wave_transform2 = wavecorr.calculate_wavelength_correction_transform( + lam_uncorr, dispersion, freference, source_xpos2, 'MOS') zero_point1 = wave_transform1(lam_uncorr) zero_point2 = wave_transform2(lam_uncorr) @@ -211,15 +213,14 @@ def test_wavecorr_fs(): lam_corr = wave2wavecorr(lam_uncorr) assert_allclose(lam_uncorr, wave2wavecorr.inverse(lam_corr)) - # test on both sides of the slit center source_xpos1 = -.2 source_xpos2 = .2 - wave_transform1 = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, - freference, source_xpos1, 'S200A1') - wave_transform2 = wavecorr.calculate_wavelength_correction_transform(lam_uncorr, dispersion, - freference, source_xpos2, 'S200A1') + wave_transform1 = wavecorr.calculate_wavelength_correction_transform( + lam_uncorr, dispersion, freference, source_xpos1, 'S200A1') + wave_transform2 = wavecorr.calculate_wavelength_correction_transform( + lam_uncorr, dispersion, freference, source_xpos2, 'S200A1') zero_point1 = wave_transform1(lam_uncorr) zero_point2 = wave_transform2(lam_uncorr) diff --git a/jwst/wavecorr/wavecorr.py b/jwst/wavecorr/wavecorr.py index 7eb0c3715a..67f8e56c35 100644 --- a/jwst/wavecorr/wavecorr.py +++ b/jwst/wavecorr/wavecorr.py @@ -138,7 +138,7 @@ def apply_zero_point_correction(slit, reffile): # Insert the new transform into the slit wcs object wave2wavecorr = Identity(2) & wave2wavecorr - slit_wcs.insert_transform('slit_frame', transform = wave2wavecorr, after=False) + slit_wcs.insert_transform('slit_frame', transform=wave2wavecorr, after=False) # Update the stored wavelengths for the slit slit.wavelength = compute_wavelength(slit_wcs) @@ -151,7 +151,7 @@ def calculate_wavelength_correction_transform(lam, dispersion, freference, Parameters ---------- - slit_wcs : ~gwcs.wcs.WCS` + slit_wcs : `~gwcs.wcs.WCS` The WCS for this slit. lam : ndarray Wavelength array [in m]. @@ -166,27 +166,23 @@ def calculate_wavelength_correction_transform(lam, dispersion, freference, Returns ------- - model : `~astropy.modeling.tabular.Tabular1D - A model which takes wavelength inputs and returns zeropoint + model : `~astropy.modeling.tabular.Tabular1D` + A model which takes wavelength inputs and returns zero-point corrected wavelengths. """ - - - # Open the zerpoint reference model + # Open the zero point reference model with datamodels.WaveCorrModel(freference) as wavecorr: for ap in wavecorr.apertures: if ap.aperture_name == aperture_name: log.info(f'Using wavelength zero-point correction for aperture {ap.aperture_name}') offset_model = ap.zero_point_offset.copy() - # TODO: implement variance - # variance = ap.variance.copy() - # width = ap.width break else: log.info(f'No wavelength zero-point correction found for slit {aperture_name}') - # Set lookup table to extrapolate at bounds to recover wavelengths beyond model bounds - # particulary for the red and blue ends of prism observations + # Set lookup table to extrapolate at bounds to recover wavelengths + # beyond model bounds, particularly for the red and blue ends of + # prism observations offset_model.bounds_error = False offset_model.fill_value = None From 86075ba1c3dbf53249d284d741b2c0cdbcb29826 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Fri, 15 Mar 2024 10:56:37 -0400 Subject: [PATCH 04/15] Updating comments --- jwst/wavecorr/wavecorr.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/jwst/wavecorr/wavecorr.py b/jwst/wavecorr/wavecorr.py index 67f8e56c35..7f33b928f0 100644 --- a/jwst/wavecorr/wavecorr.py +++ b/jwst/wavecorr/wavecorr.py @@ -151,8 +151,6 @@ def calculate_wavelength_correction_transform(lam, dispersion, freference, Parameters ---------- - slit_wcs : `~gwcs.wcs.WCS` - The WCS for this slit. lam : ndarray Wavelength array [in m]. dispersion : ndarray @@ -182,7 +180,9 @@ def calculate_wavelength_correction_transform(lam, dispersion, freference, # Set lookup table to extrapolate at bounds to recover wavelengths # beyond model bounds, particularly for the red and blue ends of - # prism observations + # prism observations. fill_value = None sets the lookup tables + # to use the default extrapolation which is a linear extrapolation + # from scipy.interpolate.interpn offset_model.bounds_error = False offset_model.fill_value = None @@ -200,7 +200,7 @@ def calculate_wavelength_correction_transform(lam, dispersion, freference, # Build a look up table to transform between corrected and uncorrected wavelengths wave2wavecorr = tabular.Tabular1D(points=lam_mean, lookup_table=lam_corrected, - bounds_error=False, fill_value=None) + bounds_error=False, fill_value=None, name='wave2wavecorr') wave2wavecorr.inverse = tabular.Tabular1D(points=lam_corrected, lookup_table=lam_mean, bounds_error=False, fill_value=None) From 257d341223d180bb958ff7e154746dd8b11a8492 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Tue, 16 Apr 2024 11:07:04 -0400 Subject: [PATCH 05/15] Added wcorr wcs frame and check on transform inverse --- jwst/wavecorr/tests/test_wavecorr.py | 22 +++++++- jwst/wavecorr/wavecorr.py | 84 +++++++++++++++++++++------- 2 files changed, 84 insertions(+), 22 deletions(-) diff --git a/jwst/wavecorr/tests/test_wavecorr.py b/jwst/wavecorr/tests/test_wavecorr.py index 047c856637..bd459bdf03 100644 --- a/jwst/wavecorr/tests/test_wavecorr.py +++ b/jwst/wavecorr/tests/test_wavecorr.py @@ -147,8 +147,26 @@ def test_skipped(): outw = WavecorrStep.call(outs) source_pos = (0.004938526981283373, -0.02795306204991911) - assert_allclose((outw.slits[ind].source_xpos, outw.slits[ind].source_ypos), source_pos) + assert_allclose((outw.slits[ind].source_xpos, outw.slits[ind].source_ypos),source_pos) + # Test that no correction transform is returned (a skip criterion) + # if the corrected wavelengths are not monotonically increasing + + # This case is not expected with real data, test with simple case + # of flipped wavelength solutions, which produces a monotonically + # decreasing wavelengths + lam = np.tile(np.flip(np.arange(0.6, 5.5, 0.01)*1e-6), (22, 1)) + disp = np.tile(np.full(lam.shape[-1], -0.01)*1e-6, (22, 1)) + + ref_name = outw.meta.ref_file.wavecorr.name + reffile = datamodels.WaveCorrModel( + WavecorrStep.reference_uri_to_cache_path(ref_name, im.crds_observatory)) + source_xpos = 0.1 + aperture_name = 'S200A1' + + transform = wavecorr.calculate_wavelength_correction_transform( + lam, disp, reffile, source_xpos, aperture_name) + assert transform is None def test_wavecorr_fs(): hdul = create_nirspec_fs_file(grating="PRISM", filter="CLEAR") @@ -199,7 +217,7 @@ def test_wavecorr_fs(): dispersion = wavecorr.compute_dispersion(slit.meta.wcs, x, y) assert_allclose(dispersion[~np.isnan(dispersion)], 1e-8, atol=1.04e-8) - # Check that the slit wcs has been updated to provide corrected wavelengths + # Check that the slit wavelengths are consistent with the slit wcs corrected_wavelength = wavecorr.compute_wavelength(slit.meta.wcs, x, y) assert_allclose(slit.wavelength, corrected_wavelength) diff --git a/jwst/wavecorr/wavecorr.py b/jwst/wavecorr/wavecorr.py index 7f33b928f0..59533ffb74 100644 --- a/jwst/wavecorr/wavecorr.py +++ b/jwst/wavecorr/wavecorr.py @@ -30,6 +30,8 @@ import logging import numpy as np from gwcs import wcstools +from gwcs import coordinate_frames as cf +from astropy import units as u from astropy.modeling import tabular from astropy.modeling.mappings import Identity @@ -91,16 +93,27 @@ def do_correction(input_model, wavecorr_file): input_model.meta.cal_step.wavecorr = 'SKIPPED' break if _is_point_source(slit, exp_type): - apply_zero_point_correction(slit, wavecorr_file) - output_model.meta.cal_step.wavecorr = 'COMPLETE' + completed = apply_zero_point_correction(slit, wavecorr_file) + if completed: + output_model.meta.cal_step.wavecorr = 'COMPLETE' + else: + log.warning(f'Corrections are not invertible for slit {slit.name}') + log.warning('Skipping wavecorr correction') + output_model.meta.cal_step.wavecorr = 'SKIPPED' + break # For MOS work on all slits containing a point source else: for slit in output_model.slits: if _is_point_source(slit, exp_type): - apply_zero_point_correction(slit, wavecorr_file) - output_model.meta.cal_step.wavecorr = 'COMPLETE' + completed = apply_zero_point_correction(slit, wavecorr_file) + if completed: + output_model.meta.cal_step.wavecorr = 'COMPLETE' + else: + log.warning(f'Corrections are not invertible for slit {slit.name}') + log.warning('Skipping wavecorr correction') + output_model.meta.cal_step.wavecorr = 'SKIPPED' return output_model @@ -114,6 +127,11 @@ def apply_zero_point_correction(slit, reffile): Slit data to be corrected. reffile : str The ``wavecorr`` reference file. + + Returns + ------- + completed : bool + A flag to report whether the zero-point correction was added or skipped. """ log.info(f'slit name {slit.name}') slit_wcs = slit.meta.wcs @@ -132,16 +150,33 @@ def apply_zero_point_correction(slit, reffile): lam = slit.wavelength.copy() * 1e-6 dispersion = compute_dispersion(slit.meta.wcs) - wave2wavecorr = calculate_wavelength_correction_transform(lam, dispersion, - reffile, source_xpos, + wave2wavecorr = calculate_wavelength_correction_transform(lam, + dispersion, + reffile, + source_xpos, aperture_name) - # Insert the new transform into the slit wcs object - wave2wavecorr = Identity(2) & wave2wavecorr - slit_wcs.insert_transform('slit_frame', transform=wave2wavecorr, after=False) - - # Update the stored wavelengths for the slit - slit.wavelength = compute_wavelength(slit_wcs) + if wave2wavecorr is None: # Should not occur for real data + completed = False + return completed + else: + # Make a new frame to insert into the slit wcs object + slit_spatial = cf.Frame2D(name='slit_spatial', axes_order=(0, 1), + unit=("", ""), axes_names=('x_slit', 'y_slit')) + spec = cf.SpectralFrame(name='spectral', axes_order=(2,), unit=(u.micron,), + axes_names=('wavelength',)) + wcorr_frame = cf.CompositeFrame( + [slit_spatial, spec], name='wcorr_frame') + + # Insert the new transform into the slit wcs object + wave2wavecorr = Identity(2) & wave2wavecorr + slit_wcs.insert_frame('slit_frame', wave2wavecorr, wcorr_frame) + + # Update the stored wavelengths for the slit + slit.wavelength = compute_wavelength(slit_wcs) + + completed = True + return completed def calculate_wavelength_correction_transform(lam, dispersion, freference, @@ -164,9 +199,10 @@ def calculate_wavelength_correction_transform(lam, dispersion, freference, Returns ------- - model : `~astropy.modeling.tabular.Tabular1D` + model : `~astropy.modeling.tabular.Tabular1D`or None A model which takes wavelength inputs and returns zero-point - corrected wavelengths. + corrected wavelengths. Returns None if an invertible model + cannot be generated. """ # Open the zero point reference model with datamodels.WaveCorrModel(freference) as wavecorr: @@ -198,13 +234,21 @@ def calculate_wavelength_correction_transform(lam, dispersion, freference, pixel_corrections = offset_model(lam_mean, source_xpos) lam_corrected = lam_mean + (pixel_corrections * disp_mean) - # Build a look up table to transform between corrected and uncorrected wavelengths - wave2wavecorr = tabular.Tabular1D(points=lam_mean, lookup_table=lam_corrected, - bounds_error=False, fill_value=None, name='wave2wavecorr') - wave2wavecorr.inverse = tabular.Tabular1D(points=lam_corrected, lookup_table=lam_mean, - bounds_error=False, fill_value=None) + # Check to make sure that the corrected wavelengths are monotonically increasing + if np.all(np.diff(lam_corrected) > 0): + # monotonically increasing + # Build a look up table to transform between corrected and uncorrected wavelengths + wave2wavecorr = tabular.Tabular1D(points=lam_mean, + lookup_table=lam_corrected, + bounds_error=False, + fill_value=None, + name='wave2wavecorr') + + return wave2wavecorr - return wave2wavecorr + else: + # output wavelengths are not monotonically increasing + return None def compute_dispersion(wcs, xpix=None, ypix=None): From e57f82f0219031230231e9b0d5044f5085241434 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Tue, 16 Apr 2024 13:04:43 -0400 Subject: [PATCH 06/15] documentation updates --- docs/jwst/wavecorr/description.rst | 22 ++++++++++++++-------- jwst/wavecorr/tests/test_wavecorr.py | 8 ++++---- jwst/wavecorr/wavecorr.py | 2 +- 3 files changed, 19 insertions(+), 13 deletions(-) diff --git a/docs/jwst/wavecorr/description.rst b/docs/jwst/wavecorr/description.rst index 999ee6a6db..fa7e311d5a 100644 --- a/docs/jwst/wavecorr/description.rst +++ b/docs/jwst/wavecorr/description.rst @@ -22,16 +22,22 @@ These are recorded in the "SRCXPOS" and "SRCYPOS" keywords in the SCI extension header of each slitlet in a FITS product. The ``wavecorr`` step loops over all slit instances in the input -science product and applies a wavelength correction to slits that -contain a point source. The point source determination is based on the -value of the "SRCTYPE" keyword populated for each slit by the +science product and updates the WCS models of slits that contain a point +source to include a wavelength correction. The point source determination is +based on the value of the "SRCTYPE" keyword populated for each slit by the :ref:`srctype ` step. The computation of the correction is based on the "SRCXPOS" value. A value of 0.0 indicates a perfectly centered source, and ranges from -0.5 to +0.5 for sources at the extreme edges of a slit. The computation uses calibration data from the ``WAVECORR`` -reference file. The correction is computed as a 2-D grid of -wavelength offsets, which is applied to the original 2-D grid of -wavelengths associated with each slit. +reference file, which contains pixel shifts as a function of source position +and wavelength, and can be converted to wavelength shifts with the dispersion. +For each slit, the ``wavecorr`` step uses the average wavelengths and +dispersions in a slit (averaged across the cross-dispersion direction) to +calculate corresponding corrected wavelengths. It then uses the average +wavelengths and their corrections to generate a transform that interpolates +between "center of slit" wavelengths and corrected wavelengths. This +transformation is added to the slit WCS after the ``slit_frame`` and +produces a new wavelength corrected slit frame, ``wavecorr_frame``. NIRSpec Fixed Slit (FS) ----------------------- @@ -47,8 +53,8 @@ the wavelength correction is only applied to the primary slit. The estimated position of the source within the primary slit (in the dispersion direction) is then used in the same manner as described above -for MOS slitlets to compute offsets to be added to the nominal wavelength -grid for the primary slit. +for MOS slitlets to update the slit WCS and compute corrected wavelengths +for the primary slit. Upon successful completion of the step, the status keyword "S_WAVCOR" is set to "COMPLETE". diff --git a/jwst/wavecorr/tests/test_wavecorr.py b/jwst/wavecorr/tests/test_wavecorr.py index bd459bdf03..4755f26c57 100644 --- a/jwst/wavecorr/tests/test_wavecorr.py +++ b/jwst/wavecorr/tests/test_wavecorr.py @@ -149,10 +149,10 @@ def test_skipped(): source_pos = (0.004938526981283373, -0.02795306204991911) assert_allclose((outw.slits[ind].source_xpos, outw.slits[ind].source_ypos),source_pos) - # Test that no correction transform is returned (a skip criterion) - # if the corrected wavelengths are not monotonically increasing - - # This case is not expected with real data, test with simple case + # Test if the corrected wavelengths are not monotonically increasing + + # This case is not expected with real data, test that no correction + # transform is returned (a skip criterion) with simple case # of flipped wavelength solutions, which produces a monotonically # decreasing wavelengths lam = np.tile(np.flip(np.arange(0.6, 5.5, 0.01)*1e-6), (22, 1)) diff --git a/jwst/wavecorr/wavecorr.py b/jwst/wavecorr/wavecorr.py index 59533ffb74..d852936de0 100644 --- a/jwst/wavecorr/wavecorr.py +++ b/jwst/wavecorr/wavecorr.py @@ -166,7 +166,7 @@ def apply_zero_point_correction(slit, reffile): spec = cf.SpectralFrame(name='spectral', axes_order=(2,), unit=(u.micron,), axes_names=('wavelength',)) wcorr_frame = cf.CompositeFrame( - [slit_spatial, spec], name='wcorr_frame') + [slit_spatial, spec], name='wavecorr_frame') # Insert the new transform into the slit wcs object wave2wavecorr = Identity(2) & wave2wavecorr From 4779ca8a92873f9907fb6a645223508764feea63 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Tue, 30 Apr 2024 17:22:18 -0400 Subject: [PATCH 07/15] updated changelog --- CHANGES.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index b9bdbf8fce..70b4896a26 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -215,6 +215,12 @@ tweakreg - Change code default to use IRAF StarFinder instead of DAO StarFinder [#8487] +wavecorr +-------- + +- Changed the NIRSpec wavelength correction algorithm to include it in slit WCS + models and resampling. Fixed the sign of the wavelength corrections. [#8376] + wfss_contam ----------- From e4205064d1b95e3ce19857b6f4999b6afd3428a7 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Fri, 31 May 2024 11:54:23 -0400 Subject: [PATCH 08/15] update get_wavelengths to use new wavecorr --- jwst/lib/tests/test_wcs_utils.py | 157 +++++++++++++++++++++++++++++++ jwst/lib/wcs_utils.py | 20 ++-- 2 files changed, 169 insertions(+), 8 deletions(-) create mode 100644 jwst/lib/tests/test_wcs_utils.py diff --git a/jwst/lib/tests/test_wcs_utils.py b/jwst/lib/tests/test_wcs_utils.py new file mode 100644 index 0000000000..0cd731e771 --- /dev/null +++ b/jwst/lib/tests/test_wcs_utils.py @@ -0,0 +1,157 @@ +import numpy as np +from numpy.testing import assert_allclose + +from astropy import units as u +from astropy import coordinates as coord +from astropy.modeling.models import Mapping, Identity, Shift, Scale +from gwcs import wcstools, wcs +from gwcs import coordinate_frames as cf + +from stdatamodels.jwst import datamodels +from stdatamodels.jwst.transforms.models import NirissSOSSModel +from jwst.lib.wcs_utils import get_wavelengths +from jwst.assign_wcs import util + + +def create_model(): + + det = cf.Frame2D(name="detector", axes_order=(0, 1)) + + sky = cf.CelestialFrame(name="sky", axes_order=(0, 1), reference_frame=coord.ICRS()) + slit_spatial = cf.Frame2D( + name="slit_spatial", + axes_order=(0, 1), + unit=("", ""), + axes_names=("x_slit", "y_slit"), + ) + + spec = cf.SpectralFrame( + name="spectral", axes_order=(2,), unit=(u.micron,), axes_names=("wavelength",) + ) + slit_frame = cf.CompositeFrame([slit_spatial, spec], name="slit_frame") + world = cf.CompositeFrame([sky, spec], name="world") + + det2slit = Mapping((0, 1, 1)) | (Identity(2) & (Scale(0.5) | Shift(0.5))) + slit2sky = Identity(3) + + slit_wcs = wcs.WCS([(det, det2slit), (slit_frame, slit2sky), (world, None)]) + + # compute wavelengths + + data = np.full((10, 10), fill_value=5.0) + + bounding_box = util.wcs_bbox_from_shape(data.shape) + + x, y = wcstools.grid_from_bounding_box(bounding_box, step=(1, 1)) + _, _, lam = slit_wcs(x, y) + lam = lam.astype(np.float32) + model = datamodels.SlitModel(data=data, wavelength=lam) + model.meta.wcs = slit_wcs + + return model + + +def create_mock_wl(): + wl = np.arange(10.0) + wl = wl[:, np.newaxis] + wl = np.repeat(wl, 10, axis=1) + wl = (wl * 0.5) + 0.5 + return wl + + +def test_get_wavelengths(): + + # create a mock SlitModel + model = create_model() + + # calculate what the wavelength array should be + wl_og = create_mock_wl() + + # Test that the get wavelengths returns the wavelength grid + wl = get_wavelengths(model) + assert_allclose(wl, wl_og) + + del model.wavelength + + # Check that wavelengths can be generated from wcs when the + # wavelength attribute is unavailable + wl = get_wavelengths(model) + assert_allclose(wl, wl_og) + + # Check that wavelengths are generated correctly when given a WFSS exp_type + wl = get_wavelengths(model, exp_type="NRC_TSGRISM") + assert_allclose(wl, wl_og) + + +def test_get_wavelengths_soss(): + + # create a mock SlitModel + model = create_model() + + del model.wavelength + model.meta.exposure.type = "NIS_SOSS" + + wcs = model.meta.wcs + new_wcs = NirissSOSSModel( + [ + 1, + ], + [ + wcs, + ], + ) + model.meta.wcs = new_wcs + + # calculate what the wavelength array should be + wl_og = create_mock_wl() + + wl = get_wavelengths(model, order=1) + assert_allclose(wl, wl_og) + + +def test_get_wavelength_wavecorr(): + + # create a mock SlitModel + model = create_model() + + wl_og = create_mock_wl() + + # Test use_wavecorr with no wavelength correction modificiation + # get_wavelengths should return the same wavelengths for use_wavecorr + # True and False + + wl_corr = get_wavelengths(model, use_wavecorr=True) + assert_allclose(wl_corr, wl_og) + + wl_uncorr = get_wavelengths(model, use_wavecorr=False) + assert_allclose(wl_corr, wl_uncorr) + + # Update the model wcs to add a wavelength corrected slit frame + slit_spatial = cf.Frame2D( + name="slit_spatial", + axes_order=(0, 1), + unit=("", ""), + axes_names=("x_slit", "y_slit"), + ) + spec = cf.SpectralFrame( + name="spectral", axes_order=(2,), unit=(u.micron,), axes_names=("wavelength",) + ) + wcorr_frame = cf.CompositeFrame([slit_spatial, spec], name="wavecorr_frame") + + # Insert the new transform into the slit wcs object + wave2wavecorr = Identity(2) & Shift(0.1) + model.meta.wcs.insert_frame("slit_frame", wave2wavecorr, wcorr_frame) + + bounding_box = util.wcs_bbox_from_shape(model.data.shape) + x, y = wcstools.grid_from_bounding_box(bounding_box, step=(1, 1)) + _, _, lam = model.meta.wcs(x, y) + model.wavelength = lam + + # calculate what the corrected wavelength array should be + wl_corr_og = wl_og + 0.1 + + wl_corr = get_wavelengths(model, use_wavecorr=True) + assert_allclose(wl_corr, wl_corr_og) + + wl_uncorr = get_wavelengths(model, use_wavecorr=False) + assert_allclose(wl_uncorr, wl_og) diff --git a/jwst/lib/wcs_utils.py b/jwst/lib/wcs_utils.py index 5dbfceeffb..a89e729f4d 100644 --- a/jwst/lib/wcs_utils.py +++ b/jwst/lib/wcs_utils.py @@ -41,19 +41,23 @@ def get_wavelengths(model, exp_type="", order=None, use_wavecorr=None): got_wavelength = False wl_array = None - # If we've been asked to use the corrected wavelengths stored in - # the wavelength array, return those wavelengths. Otherwise, the - # results computed from the WCS object (below) will be returned. - if use_wavecorr is not None: - if use_wavecorr: - return wl_array - else: - got_wavelength = False # force wl computation below # Evaluate the WCS on the grid of pixel indexes, capturing only the # resulting wavelength values shape = model.data.shape grid = np.indices(shape[-2:], dtype=np.float64) + + # If we've been asked to use the uncorrected wavelengths we need to + # recalculate them from the wcs by skipping the transformation between + # the slit frame and the wavelength corrected slit frame. If the wavecorr_frame + # is not in the wcs assume that the wavelength correction has not been applied. + if use_wavecorr is not None: + if not use_wavecorr and 'wavecorr_frame' in model.meta.wcs.available_frames: + wcs = model.meta.wcs + detector2slit = wcs.get_transform('detector', 'slit_frame') + wavecorr2world = wcs.get_transform("wavecorr_frame", "world") + wl_array = (detector2slit | wavecorr2world)(grid[1], grid[0])[2] + return wl_array # If no existing wavelength array, compute one if hasattr(model.meta, "wcs") and not got_wavelength: From 280f6dfdd99f4f835c1634cb19e00d0d9a42b64f Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Fri, 31 May 2024 12:02:54 -0400 Subject: [PATCH 09/15] refactor FS flatfield and pathloss to use get_wavelengths --- jwst/flatfield/flat_field.py | 58 ++++-------------------------------- jwst/pathloss/pathloss.py | 45 +++++++++++++++++++--------- jwst/wavecorr/wavecorr.py | 8 ++--- 3 files changed, 40 insertions(+), 71 deletions(-) diff --git a/jwst/flatfield/flat_field.py b/jwst/flatfield/flat_field.py index 1c4c1db6ea..576bf1dead 100644 --- a/jwst/flatfield/flat_field.py +++ b/jwst/flatfield/flat_field.py @@ -11,7 +11,7 @@ from stdatamodels.jwst import datamodels from stdatamodels.jwst.datamodels import dqflags -from ..lib import reffile_utils +from ..lib import reffile_utils, wcs_utils from ..assign_wcs import nirspec log = logging.getLogger(__name__) @@ -1924,59 +1924,11 @@ def flat_for_nirspec_slit(slit, f_flat_model, s_flat_model, d_flat_model, # Get the wavelength at each pixel in the extracted slit data. # If the wavelength attribute exists and is populated, use it # in preference to the wavelengths returned by the wcs function. - got_wl_attribute = True - try: - wl = slit.wavelength.copy() # a 2-D array - except AttributeError: - got_wl_attribute = False - if not got_wl_attribute or len(wl) == 0: - got_wl_attribute = False - return_dummy = False - # Has the use_wavecorr param been set? - if use_wavecorr is not None: - if use_wavecorr: - # Need to use the 2D wavelength array, because that's where - # the corrected wavelengths are stored - if got_wl_attribute: - # We've got the "wl" wavelength array we need - pass - else: - # Can't do the computation without the 2D wavelength array - log.error(f"The wavelength array for slit {slit.name} is not populated") - log.error("Skipping flat-field correction") - return_dummy = True - elif not use_wavecorr: - # Need to use the WCS object to create an uncorrected 2D wavelength array - if got_wcs: - log.info(f"Creating wavelength array from WCS for slit {slit.name}") - bb = slit.meta.wcs.bounding_box - grid = grid_from_bounding_box(bb) - wl = slit.meta.wcs(*grid)[2] - del grid - else: - # Can't create the uncorrected wavelengths without the WCS - log.error(f"Slit {slit.name} has no WCS object") - log.error("Skipping flat-field correction") - return_dummy = True - else: - # use_wavecorr was not specified, so use default processing - if not got_wl_attribute or np.nanmin(wl) == 0. and np.nanmax(wl) == 0.: - got_wl_attribute = False - log.warning(f"The wavelength array for slit {slit.name} has not been populated") - # Try to create it from the WCS - if got_wcs: - bb = slit.meta.wcs.bounding_box - grid = grid_from_bounding_box(bb) - wl = slit.meta.wcs(*grid)[2] - del grid - else: - log.warning("and this slit does not have a 'wcs' attribute") - log.warning("likely because assign_wcs has not been run.") - log.error("skipping ...") - return_dummy = True - else: - log.debug("Wavelengths are from the wavelength array.") + return_dummy = False + wl = wcs_utils.get_wavelengths(slit, use_wavecorr=use_wavecorr) + if wl is None: + return_dummy = True # Create and return a dummy flat as a placeholder, if necessary if return_dummy: diff --git a/jwst/pathloss/pathloss.py b/jwst/pathloss/pathloss.py index 4ad2f551b8..9586bc320e 100644 --- a/jwst/pathloss/pathloss.py +++ b/jwst/pathloss/pathloss.py @@ -976,24 +976,41 @@ def _corrections_for_fixedslit(slit, pathloss, exp_type, source_type): wavelength_pointsource *= 1.0e6 wavelength_uniformsource *= 1.0e6 - wavelength_array = slit.wavelength - - # Compute the point source pathloss 2D correction - pathloss_2d_ps = interpolate_onto_grid( - wavelength_array, - wavelength_pointsource, - pathloss_pointsource_vector) - - # Compute the uniform source pathloss 2D correction - pathloss_2d_un = interpolate_onto_grid( - wavelength_array, - wavelength_uniformsource, - pathloss_uniform_vector) - # Use the appropriate correction for this slit if is_pointsource(source_type or slit.source_type): + # calculate the point source corrected wavelengths and uncorrected wavelengths for the slit + wavelength_array_corr = get_wavelengths(slit, use_wavecorr=True) + wavelength_array_uncorr = get_wavelengths(slit, use_wavecorr=False) + + # Compute the point source pathloss 2D correction + pathloss_2d_ps = interpolate_onto_grid( + wavelength_array_corr, + wavelength_pointsource, + pathloss_pointsource_vector) + + # Compute the uniform source pathloss 2D correction + pathloss_2d_un = interpolate_onto_grid( + wavelength_array_uncorr, + wavelength_uniformsource, + pathloss_uniform_vector) + pathloss_2d = pathloss_2d_ps + else: + wavelength_array = slit.wavelength + + # Compute the point source pathloss 2D correction + pathloss_2d_ps = interpolate_onto_grid( + wavelength_array, + wavelength_pointsource, + pathloss_pointsource_vector) + + # Compute the uniform source pathloss 2D correction + pathloss_2d_un = interpolate_onto_grid( + wavelength_array, + wavelength_uniformsource, + pathloss_uniform_vector) + pathloss_2d = pathloss_2d_un # Save the corrections. The `data` portion is the correction used. diff --git a/jwst/wavecorr/wavecorr.py b/jwst/wavecorr/wavecorr.py index d852936de0..c4efb74299 100644 --- a/jwst/wavecorr/wavecorr.py +++ b/jwst/wavecorr/wavecorr.py @@ -96,7 +96,7 @@ def do_correction(input_model, wavecorr_file): completed = apply_zero_point_correction(slit, wavecorr_file) if completed: output_model.meta.cal_step.wavecorr = 'COMPLETE' - else: + else: # pragma: no cover log.warning(f'Corrections are not invertible for slit {slit.name}') log.warning('Skipping wavecorr correction') output_model.meta.cal_step.wavecorr = 'SKIPPED' @@ -110,7 +110,7 @@ def do_correction(input_model, wavecorr_file): completed = apply_zero_point_correction(slit, wavecorr_file) if completed: output_model.meta.cal_step.wavecorr = 'COMPLETE' - else: + else: # pragma: no cover log.warning(f'Corrections are not invertible for slit {slit.name}') log.warning('Skipping wavecorr correction') output_model.meta.cal_step.wavecorr = 'SKIPPED' @@ -155,8 +155,8 @@ def apply_zero_point_correction(slit, reffile): reffile, source_xpos, aperture_name) - - if wave2wavecorr is None: # Should not occur for real data + # wave2wavecorr should not be None for real data + if wave2wavecorr is None: # pragma: no cover completed = False return completed else: From 70e72a054734f1233a49f9224a113ed8329ed99e Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Mon, 3 Jun 2024 11:38:28 -0400 Subject: [PATCH 10/15] updated documentation and the changelog --- CHANGES.rst | 17 +++++++++++++++++ docs/jwst/pathloss/description.rst | 4 ++++ 2 files changed, 21 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 70b4896a26..cb39203416 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -93,6 +93,9 @@ flat_field - Update NIRSpec flatfield code for all modes to ensure SCI=ERR=NaN wherever the DO_NOT_USE flag is set in the DQ array. [#8463] +- Updated the flatfield code to use the new format of the ``wavecorr`` +wavelength zero-point corrections for point sources. [#8376] + general ------- @@ -101,6 +104,13 @@ general - Increase minimum required scipy. [#8441] +lib +--- + +- Updated the ``wcs_utils.get_wavelength`` to use the new format +ofthe ``wavecorr`` wavelength zero-point corrections for point +sources. [#8376] + outlier_detection ----------------- @@ -127,6 +137,13 @@ outlier_detection to detect outliers in TSO data, with user-defined rolling window width via the ``rolling_window_width`` parameter. [#8473] +pathloss +-------- + +- Updated pathloss calculations for fixed slit mode to use the appropriate +wavelengths for point and uniform sources if the ``wavecorr`` wavelength +zero-point corrections for point sources have been applied. [#8376] + photom ------ diff --git a/docs/jwst/pathloss/description.rst b/docs/jwst/pathloss/description.rst index 2e10ad4cbe..65529b736a 100644 --- a/docs/jwst/pathloss/description.rst +++ b/docs/jwst/pathloss/description.rst @@ -45,6 +45,10 @@ Once the 1-D correction arrays have been computed, both forms of the correction (point and uniform) are interpolated, as a function of wavelength, into the 2-D space of the slit or IFU data and attached to the output data model (extensions "PATHLOSS_PS" and "PATHLOSS_UN") as a record of what was computed. +For fixed slit data, if the ``wavecorr`` step has been run to provide wavelength +corrections to point sources, the corrected wavelengths will be used to +calculate the point source pathloss, whereas the uncorrected wavelengths (appropriate +for uniform sources) will be used to calculate the uniform source pathlosses. The form of the 2-D correction (point or uniform) that's appropriate for the data is divided into the SCI and ERR arrays and propagated into the variance arrays of the science data. From a6f1175a8588e90fe55644d60cff3ebe66485a1a Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Mon, 3 Jun 2024 11:53:42 -0400 Subject: [PATCH 11/15] remove unused import --- jwst/flatfield/flat_field.py | 1 - 1 file changed, 1 deletion(-) diff --git a/jwst/flatfield/flat_field.py b/jwst/flatfield/flat_field.py index 576bf1dead..01da31668d 100644 --- a/jwst/flatfield/flat_field.py +++ b/jwst/flatfield/flat_field.py @@ -6,7 +6,6 @@ import math import numpy as np -from gwcs.wcstools import grid_from_bounding_box from stdatamodels.jwst import datamodels from stdatamodels.jwst.datamodels import dqflags From 910210f9e06312479b761d598b6097869e146228 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Mon, 3 Jun 2024 14:07:43 -0400 Subject: [PATCH 12/15] record wavecorr completion for each mos slit --- jwst/wavecorr/tests/test_wavecorr.py | 39 ++++++++++++++++++++++++++++ jwst/wavecorr/wavecorr.py | 8 ++++-- 2 files changed, 45 insertions(+), 2 deletions(-) diff --git a/jwst/wavecorr/tests/test_wavecorr.py b/jwst/wavecorr/tests/test_wavecorr.py index 4755f26c57..5ff61fc241 100644 --- a/jwst/wavecorr/tests/test_wavecorr.py +++ b/jwst/wavecorr/tests/test_wavecorr.py @@ -167,6 +167,45 @@ def test_skipped(): transform = wavecorr.calculate_wavelength_correction_transform( lam, disp, reffile, source_xpos, aperture_name) assert transform is None + + +def test_mos_slit_status(): + """ Test conditions that are skipped for mos slitlets.""" + + hdul = create_nirspec_mos_file() + msa_meta = os.path.join(jwst.__path__[0], *['assign_wcs', 'tests', 'data', 'msa_configuration.fits']) + hdul[0].header['MSAMETFL'] = msa_meta + hdul[0].header['MSAMETID'] = 12 + im = datamodels.ImageModel(hdul) + im_wcs = AssignWcsStep.call(im) + im_ex2d = Extract2dStep.call(im_wcs) + bbox = ((-.5, 1432.5), (-.5, 37.5)) + im_ex2d.slits[0].meta.wcs.bounding_box = bbox + x, y = wcstools.grid_from_bounding_box(bbox) + ra, dec, lam_before = im_ex2d.slits[0].meta.wcs(x, y) + im_ex2d.slits[0].wavelength = lam_before + im_src = SourceTypeStep.call(im_ex2d) + + # test the mock msa source as an extended source + im_src.slits[0].source_type = 'EXTENDED' + im_wave = WavecorrStep.call(im_src) + + # check that the step is recorded as completed + assert im_wave.meta.cal_step.wavecorr == 'COMPLETE' + + # check that the step is listed as skipped for extended mos sources + assert im_wave.slits[0].meta.cal_step.wavecorr == 'SKIPPED' + + # test the mock msa source as a point source + im_src.slits[0].source_type = 'POINT' + im_wave = WavecorrStep.call(im_src) + + # check that the step is recorded as completed + assert im_wave.meta.cal_step.wavecorr == 'COMPLETE' + + # check that the step is listed as complete for mos point sources + assert im_wave.slits[0].meta.cal_step.wavecorr == 'COMPLETE' + def test_wavecorr_fs(): hdul = create_nirspec_fs_file(grating="PRISM", filter="CLEAR") diff --git a/jwst/wavecorr/wavecorr.py b/jwst/wavecorr/wavecorr.py index c4efb74299..46a92513b2 100644 --- a/jwst/wavecorr/wavecorr.py +++ b/jwst/wavecorr/wavecorr.py @@ -109,11 +109,15 @@ def do_correction(input_model, wavecorr_file): if _is_point_source(slit, exp_type): completed = apply_zero_point_correction(slit, wavecorr_file) if completed: - output_model.meta.cal_step.wavecorr = 'COMPLETE' + slit.meta.cal_step.wavecorr = 'COMPLETE' else: # pragma: no cover log.warning(f'Corrections are not invertible for slit {slit.name}') log.warning('Skipping wavecorr correction') - output_model.meta.cal_step.wavecorr = 'SKIPPED' + slit.meta.cal_step.wavecorr = 'SKIPPED' + else: + slit.meta.cal_step.wavecorr = 'SKIPPED' + + output_model.meta.cal_step.wavecorr = 'COMPLETE' return output_model From d03ef6f42a01bea107bd8a2318da4cb76fdfd329 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Tue, 4 Jun 2024 17:44:53 -0400 Subject: [PATCH 13/15] Minor change log updates Co-authored-by: Melanie Clarke --- CHANGES.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index cb39203416..d6050da468 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -94,7 +94,7 @@ flat_field DO_NOT_USE flag is set in the DQ array. [#8463] - Updated the flatfield code to use the new format of the ``wavecorr`` -wavelength zero-point corrections for point sources. [#8376] + wavelength zero-point corrections for point sources. [#8376] general ------- @@ -108,8 +108,8 @@ lib --- - Updated the ``wcs_utils.get_wavelength`` to use the new format -ofthe ``wavecorr`` wavelength zero-point corrections for point -sources. [#8376] + of the ``wavecorr`` wavelength zero-point corrections for point + sources. [#8376] outlier_detection ----------------- @@ -141,8 +141,8 @@ pathloss -------- - Updated pathloss calculations for fixed slit mode to use the appropriate -wavelengths for point and uniform sources if the ``wavecorr`` wavelength -zero-point corrections for point sources have been applied. [#8376] + wavelengths for point and uniform sources if the ``wavecorr`` wavelength + zero-point corrections for point sources have been applied. [#8376] photom ------ From e66388f1678b28dce40e32cf5935088431d57c50 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Wed, 5 Jun 2024 09:26:52 -0400 Subject: [PATCH 14/15] catch missing wcs case in get_wavelengths Co-authored-by: Melanie Clarke --- jwst/lib/wcs_utils.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/jwst/lib/wcs_utils.py b/jwst/lib/wcs_utils.py index a89e729f4d..3070756f13 100644 --- a/jwst/lib/wcs_utils.py +++ b/jwst/lib/wcs_utils.py @@ -52,7 +52,8 @@ def get_wavelengths(model, exp_type="", order=None, use_wavecorr=None): # the slit frame and the wavelength corrected slit frame. If the wavecorr_frame # is not in the wcs assume that the wavelength correction has not been applied. if use_wavecorr is not None: - if not use_wavecorr and 'wavecorr_frame' in model.meta.wcs.available_frames: + if (not use_wavecorr and hasattr(model.meta, "wcs") + and 'wavecorr_frame' in model.meta.wcs.available_frames): wcs = model.meta.wcs detector2slit = wcs.get_transform('detector', 'slit_frame') wavecorr2world = wcs.get_transform("wavecorr_frame", "world") From 772ed05b8dacd7ef2c3fe1ee4c70723abad76578 Mon Sep 17 00:00:00 2001 From: Christian Hayes Date: Fri, 7 Jun 2024 11:15:22 -0400 Subject: [PATCH 15/15] specify NIRSpec in change log Co-authored-by: Howard Bushouse --- CHANGES.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 61be417cec..08c5e3fe30 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -108,7 +108,7 @@ flat_field - Update NIRSpec flatfield code for all modes to ensure SCI=ERR=NaN wherever the DO_NOT_USE flag is set in the DQ array. [#8463] -- Updated the flatfield code to use the new format of the ``wavecorr`` +- Updated the NIRSpec flatfield code to use the new format of the ``wavecorr`` wavelength zero-point corrections for point sources. [#8376] general @@ -124,7 +124,7 @@ lib - Updated the ``wcs_utils.get_wavelength`` to use the new format of the ``wavecorr`` wavelength zero-point corrections for point - sources. [#8376] + sources in NIRSpec slit data. [#8376] master_background_mos --------------------- @@ -162,7 +162,7 @@ outlier_detection pathloss -------- -- Updated pathloss calculations for fixed slit mode to use the appropriate +- Updated pathloss calculations for NIRSpec fixed slit mode to use the appropriate wavelengths for point and uniform sources if the ``wavecorr`` wavelength zero-point corrections for point sources have been applied. [#8376]