From c83b30c2933335d19a83728b340007bf4a179e11 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 10 Jan 2024 14:40:47 -0500 Subject: [PATCH 01/69] Update jump.py --- src/stcal/jump/jump.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 6139d103..a98d70c3 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -2,6 +2,7 @@ import multiprocessing import time +import astropy.stats import cv2 as cv import numpy as np from astropy import stats @@ -838,7 +839,7 @@ def find_faint_extended( data = indata.copy() data[gdq == sat_flag] = np.nan data[gdq == 1] = np.nan - data[gdq == jump_flag] = np.nan +# data[gdq == jump_flag] = np.nan all_ellipses = [] first_diffs = np.diff(data, axis=1) first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) @@ -879,7 +880,9 @@ def find_faint_extended( nrows = ratio.shape[1] ncols = ratio.shape[2] extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) - exty, extx = np.where(masked_smoothed_ratio > snr_threshold) + mean, median, stddev = astropy.stats.sigma_clipped_stats(masked_smoothed_ratio) + cutoff = median + snr_threshold * stddev + exty, extx = np.where(masked_smoothed_ratio > cutoff) extended_emission[exty, extx] = 1 # find the contours of the extended emission contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) From d84d449249e682d4ac56e42ca1ba9c64f851a1a8 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Tue, 16 Jan 2024 15:52:11 -0500 Subject: [PATCH 02/69] fix missing pdq info --- src/stcal/jump/jump.py | 21 ++++++++++++++------- tests/test_jump.py | 12 +++++++++--- 2 files changed, 23 insertions(+), 10 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index a98d70c3..05f123bc 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -299,6 +299,7 @@ def detect_jumps( gdq, num_showers = find_faint_extended( data, gdq, + pdq, readnoise_2d, frames_per_group, minimum_sigclip_groups, @@ -449,6 +450,7 @@ def detect_jumps( gdq, num_showers = find_faint_extended( data, gdq, + pdq, readnoise_2d, frames_per_group, minimum_sigclip_groups, @@ -782,6 +784,7 @@ def near_edge(jump, low_threshold, high_threshold): def find_faint_extended( indata, gdq, + pdq, readnoise_2d, nframes, minimum_sigclip_groups, @@ -868,20 +871,24 @@ def find_faint_extended( # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] masked_ratio = ratio[grp - 1].copy() - jumpy, jumpx = np.where(gdq[intg, grp, :, :] == jump_flag) - # mask pix. that are already flagged as jump - masked_ratio[jumpy, jumpx] = np.nan - saty, satx = np.where(gdq[intg, grp, :, :] == sat_flag) + # mask pixels that are already flagged as jump + combined_pixel_mask = np.bitwise_or(gdq[intg, grp, :, :], pdq[intg, grp, :, :]) + jump_pixels_array = np.bitwise_and(combined_pixel_mask, jump_flag) + jumpy, jumpx = np.where(jump_pixels_array == jump_flag) + masked_ratio[jumpy, jumpx] = np.nan - # mask pix. that are already flagged as sat. + # mask pixels that are already flagged as sat. + sat_pixels_array = np.bitwise_and(combined_pixel_mask, sat_flag) + saty, satx = np.where(sat_pixels_array == sat_flag) masked_ratio[saty, satx] = np.nan + masked_smoothed_ratio = convolve(masked_ratio, ring_2D_kernel) nrows = ratio.shape[1] ncols = ratio.shape[2] extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) - mean, median, stddev = astropy.stats.sigma_clipped_stats(masked_smoothed_ratio) - cutoff = median + snr_threshold * stddev + mean2, median2, stddev2 = astropy.stats.sigma_clipped_stats(masked_smoothed_ratio) + cutoff = median2 + snr_threshold * stddev2 exty, extx = np.where(masked_smoothed_ratio > cutoff) extended_emission[exty, extx] = 1 # find the contours of the extended emission diff --git a/tests/test_jump.py b/tests/test_jump.py index 0ddbefb1..d1e6ae70 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -1,6 +1,6 @@ import numpy as np import pytest - +from astropy.io import fits from stcal.jump.jump import ( calc_num_slices, extend_saturation, @@ -214,14 +214,17 @@ def test_find_faint_extended(): nint, ngrps, ncols, nrows = 1, 6, 30, 30 data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) gdq = np.zeros_like(data, dtype=np.uint8) + pdq = np.zeros_like(data, dtype=np.uint8) gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain rng = np.random.default_rng(12345) - data[0, 1:, 14:20, 15:20] = 6 * gain * 1.7 + data[0, 1:2, 14:20, 15:20] = 6 * gain * 2.0 data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise + fits.writeto("data.fits", data, overwrite=True) gdq, num_showers = find_faint_extended( data, gdq, + pdq, readnoise, 1, 100, @@ -232,8 +235,9 @@ def test_find_faint_extended(): sat_flag=2, jump_flag=4, ellipse_expand=1.1, - num_grps_masked=3, + num_grps_masked=0, ) + fits.writeto("outputgdq.fits", gdq, overwrite=True) # Check that all the expected samples in group 2 are flagged as jump and # that they are not flagged outside assert num_showers == 3 @@ -274,6 +278,7 @@ def test_find_faint_extended_sigclip(): nint, ngrps, ncols, nrows = 101, 6, 30, 30 data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) gdq = np.zeros_like(data, dtype=np.uint8) + pdq = np.zeros_like(data, dtype=np.uint8) gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain rng = np.random.default_rng(12345) @@ -282,6 +287,7 @@ def test_find_faint_extended_sigclip(): gdq, num_showers = find_faint_extended( data, gdq, + pdq, readnoise, 1, 100, From 01f930a6e808c254ab73381cebb6dc61307d33e1 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Tue, 16 Jan 2024 16:20:02 -0500 Subject: [PATCH 03/69] test --- src/stcal/jump/jump.py | 4 +++- tests/test_jump.py | 14 +++++++------- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 05f123bc..359c5c48 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -1,7 +1,7 @@ import logging import multiprocessing import time - +from astropy.io import fits import astropy.stats import cv2 as cv import numpy as np @@ -891,6 +891,8 @@ def find_faint_extended( cutoff = median2 + snr_threshold * stddev2 exty, extx = np.where(masked_smoothed_ratio > cutoff) extended_emission[exty, extx] = 1 + if grp == 2: + fits.writeto("extendedemission.fits", extended_emission, overwrite=True) # find the contours of the extended emission contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) # get the contours that are above the minimum size diff --git a/tests/test_jump.py b/tests/test_jump.py index d1e6ae70..28799314 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -211,27 +211,27 @@ def test_flag_large_events_withsnowball_noextension(): def test_find_faint_extended(): - nint, ngrps, ncols, nrows = 1, 6, 30, 30 + nint, ngrps, ncols, nrows = 1, 66, 30, 30 data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) gdq = np.zeros_like(data, dtype=np.uint8) pdq = np.zeros_like(data, dtype=np.uint8) gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain rng = np.random.default_rng(12345) - data[0, 1:2, 14:20, 15:20] = 6 * gain * 2.0 + data[0, 1:2, 14:20, 15:20] = 6 * gain * 4.0 data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise fits.writeto("data.fits", data, overwrite=True) gdq, num_showers = find_faint_extended( data, gdq, pdq, - readnoise, + readnoise * np.sqrt(2), 1, 100, - snr_threshold=1.3, - min_shower_area=20, + snr_threshold=3, + min_shower_area=10, inner=1, - outer=2, + outer=2.6, sat_flag=2, jump_flag=4, ellipse_expand=1.1, @@ -240,7 +240,7 @@ def test_find_faint_extended(): fits.writeto("outputgdq.fits", gdq, overwrite=True) # Check that all the expected samples in group 2 are flagged as jump and # that they are not flagged outside - assert num_showers == 3 + assert num_showers == 1 assert np.all(gdq[0, 1, 22, 14:23] == 0) assert np.all(gdq[0, 1, 21, 16:20] == DQFLAGS["JUMP_DET"]) assert np.all(gdq[0, 1, 20, 15:22] == DQFLAGS["JUMP_DET"]) From 1a4424c6ce89c5baaf2d3386325d1d55a490b014 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Tue, 16 Jan 2024 17:08:05 -0500 Subject: [PATCH 04/69] updates --- src/stcal/jump/jump.py | 15 +++++++++++---- tests/test_jump.py | 4 +++- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 359c5c48..862c201a 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -840,11 +840,14 @@ def find_faint_extended( """ read_noise_2 = readnoise_2d**2 data = indata.copy() - data[gdq == sat_flag] = np.nan - data[gdq == 1] = np.nan -# data[gdq == jump_flag] = np.nan - all_ellipses = [] + bad_pixels_array = np.bitwise_and(pdq, 1) + dnu_int, dnu_grp, dnuy, dnux = np.where(bad_pixels_array == 1) + data[dnu_int, dnu_grp, dnuy, dnux] = np.nan + first_diffs = np.diff(data, axis=1) + fits.writeto("first_diffs.fits", first_diffs, overwrite=True) + all_ellipses = [] + first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) nints = data.shape[0] if nints > minimum_sigclip_groups: @@ -891,7 +894,11 @@ def find_faint_extended( cutoff = median2 + snr_threshold * stddev2 exty, extx = np.where(masked_smoothed_ratio > cutoff) extended_emission[exty, extx] = 1 + if grp == 1: + fits.writeto("masked_smoothed_ratio1.fits", masked_smoothed_ratio, overwrite=True) + fits.writeto("extendedemission1.fits", extended_emission, overwrite=True) if grp == 2: + fits.writeto("masked_smoothed_ratio.fits", masked_smoothed_ratio, overwrite=True) fits.writeto("extendedemission.fits", extended_emission, overwrite=True) # find the contours of the extended emission contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) diff --git a/tests/test_jump.py b/tests/test_jump.py index 28799314..91e36d9c 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -215,10 +215,12 @@ def test_find_faint_extended(): data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) gdq = np.zeros_like(data, dtype=np.uint8) pdq = np.zeros_like(data, dtype=np.uint8) + pdq[:, :, 0, 0] = 1 +# pdq = np.zeros(shape=(data.shape[2], data.shape[3]), dtype=np.uint8) gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain rng = np.random.default_rng(12345) - data[0, 1:2, 14:20, 15:20] = 6 * gain * 4.0 + data[0, 1:, 14:20, 15:20] = 6 * gain * 4.0 data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise fits.writeto("data.fits", data, overwrite=True) gdq, num_showers = find_faint_extended( From bc223a5b8bababff345a067a1cdca16f1f39f7fd Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 11:45:43 -0500 Subject: [PATCH 05/69] Update test_jump.py --- tests/test_jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_jump.py b/tests/test_jump.py index 91e36d9c..e954efc4 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -220,7 +220,7 @@ def test_find_faint_extended(): gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain rng = np.random.default_rng(12345) - data[0, 1:, 14:20, 15:20] = 6 * gain * 4.0 + data[0, 1:, 14:20, 15:20] = 6 * gain * 6.0 * np.sqrt(2) data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise fits.writeto("data.fits", data, overwrite=True) gdq, num_showers = find_faint_extended( From 70662f89ea4ac245cc46281d243eb3542de0c151 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 11:46:44 -0500 Subject: [PATCH 06/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 862c201a..05d99f31 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -674,7 +674,7 @@ def extend_ellipses( ) jump_ellipse = image[:, :, 2] ngrps = gdq_cube.shape[1] - last_grp = min(grp + num_grps_masked, ngrps) + last_grp = find_last_grp(grp, ngrps, num_grps_masked) # This loop will flag the number of groups for flg_grp in range(grp, last_grp): sat_pix = np.bitwise_and(gdq_cube[intg, flg_grp, :, :], sat_flag) From 6d14d3b656e0c473d2a4c5da3d63fbcc2733b11c Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 11:53:52 -0500 Subject: [PATCH 07/69] Update jump.py --- src/stcal/jump/jump.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 05d99f31..69c4149d 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -683,6 +683,10 @@ def extend_ellipses( gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) return gdq_cube, num_ellipses +def find_last_grp(grp, ngrps, num_grps_masked_after): + num_grps_masked_after += 1 + last_grp = min(grp + num_grps_masked_after, ngrps) + return last_grp def find_circles(dqplane, bitmask, min_area): # Using an input DQ plane this routine will find the groups of pixels with at least the minimum From 5832948340b1d984aa02df357d9a4e38437e7f97 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 12:02:05 -0500 Subject: [PATCH 08/69] fixed pdq --- src/stcal/jump/jump.py | 6 +++--- tests/test_jump.py | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 69c4149d..88659879 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -845,8 +845,8 @@ def find_faint_extended( read_noise_2 = readnoise_2d**2 data = indata.copy() bad_pixels_array = np.bitwise_and(pdq, 1) - dnu_int, dnu_grp, dnuy, dnux = np.where(bad_pixels_array == 1) - data[dnu_int, dnu_grp, dnuy, dnux] = np.nan + dnuy, dnux = np.where(bad_pixels_array == 1) + data[:, :, dnuy, dnux] = np.nan first_diffs = np.diff(data, axis=1) fits.writeto("first_diffs.fits", first_diffs, overwrite=True) @@ -880,7 +880,7 @@ def find_faint_extended( masked_ratio = ratio[grp - 1].copy() # mask pixels that are already flagged as jump - combined_pixel_mask = np.bitwise_or(gdq[intg, grp, :, :], pdq[intg, grp, :, :]) + combined_pixel_mask = np.bitwise_or(gdq[intg, grp, :, :], pdq[:, :]) jump_pixels_array = np.bitwise_and(combined_pixel_mask, jump_flag) jumpy, jumpx = np.where(jump_pixels_array == jump_flag) masked_ratio[jumpy, jumpx] = np.nan diff --git a/tests/test_jump.py b/tests/test_jump.py index e954efc4..9f177da5 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -214,8 +214,8 @@ def test_find_faint_extended(): nint, ngrps, ncols, nrows = 1, 66, 30, 30 data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) gdq = np.zeros_like(data, dtype=np.uint8) - pdq = np.zeros_like(data, dtype=np.uint8) - pdq[:, :, 0, 0] = 1 + pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint8) + pdq[0, 0] = 1 # pdq = np.zeros(shape=(data.shape[2], data.shape[3]), dtype=np.uint8) gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain From 963f9d206bbeb514d196e7903f291942bcd8d416 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 12:10:23 -0500 Subject: [PATCH 09/69] Update jump.py --- src/stcal/jump/jump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 88659879..5c40526d 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -847,6 +847,7 @@ def find_faint_extended( bad_pixels_array = np.bitwise_and(pdq, 1) dnuy, dnux = np.where(bad_pixels_array == 1) data[:, :, dnuy, dnux] = np.nan + gdq[:, :, dnuy, dnux] = np.nan first_diffs = np.diff(data, axis=1) fits.writeto("first_diffs.fits", first_diffs, overwrite=True) From e12f98cb4ac7005b7e2fd874e2c56ba4eab0aec5 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 12:20:26 -0500 Subject: [PATCH 10/69] Update jump.py --- src/stcal/jump/jump.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 5c40526d..afc5b7d4 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -787,7 +787,7 @@ def near_edge(jump, low_threshold, high_threshold): def find_faint_extended( indata, - gdq, + ingdq, pdq, readnoise_2d, nframes, @@ -842,12 +842,13 @@ def find_faint_extended( Total number of showers detected. """ - read_noise_2 = readnoise_2d**2 + gdq = ingdq.copy() data = indata.copy() + read_noise_2 = readnoise_2d**2 bad_pixels_array = np.bitwise_and(pdq, 1) dnuy, dnux = np.where(bad_pixels_array == 1) data[:, :, dnuy, dnux] = np.nan - gdq[:, :, dnuy, dnux] = np.nan + gdq = np.bitwise_or(bad_pixels_array[np.newaxis, np.newaxis, :, :], gdq) first_diffs = np.diff(data, axis=1) fits.writeto("first_diffs.fits", first_diffs, overwrite=True) From 616785f41ca5ffcec2874e77688938bf12bba17d Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 12:31:02 -0500 Subject: [PATCH 11/69] Update jump.py --- src/stcal/jump/jump.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index afc5b7d4..6e3d5d84 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -842,14 +842,17 @@ def find_faint_extended( Total number of showers detected. """ + fits.writeto("ingdq.fits", ingdq, overwrite=True) + fits.writeto("inpdg.fits", pdq, overwrite=True) gdq = ingdq.copy() data = indata.copy() read_noise_2 = readnoise_2d**2 bad_pixels_array = np.bitwise_and(pdq, 1) + fits.writeto("badpixelsarray.fits", bad_pixels_array, overwrite=True) dnuy, dnux = np.where(bad_pixels_array == 1) data[:, :, dnuy, dnux] = np.nan gdq = np.bitwise_or(bad_pixels_array[np.newaxis, np.newaxis, :, :], gdq) - + fits.writeto("updategdq.fits". gdq, overwrite=True) first_diffs = np.diff(data, axis=1) fits.writeto("first_diffs.fits", first_diffs, overwrite=True) all_ellipses = [] From eb196a881fb73bd77a86625f39af056342572e16 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 12:34:28 -0500 Subject: [PATCH 12/69] Update jump.py --- src/stcal/jump/jump.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 6e3d5d84..c9ce0f31 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -895,6 +895,11 @@ def find_faint_extended( saty, satx = np.where(sat_pixels_array == sat_flag) masked_ratio[saty, satx] = np.nan + # mask pixels that are already flagged as do not use + sat_pixels_array = np.bitwise_and(combined_pixel_mask, 1) + dnuy, dnux = np.where(sat_pixels_array == 1) + masked_ratio[saty, satx] = np.nan + masked_smoothed_ratio = convolve(masked_ratio, ring_2D_kernel) nrows = ratio.shape[1] ncols = ratio.shape[2] From 5b95710d0df6daaa388259e31fcaee3ab157c3e6 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 12:38:43 -0500 Subject: [PATCH 13/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index c9ce0f31..8c491ab6 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -852,7 +852,7 @@ def find_faint_extended( dnuy, dnux = np.where(bad_pixels_array == 1) data[:, :, dnuy, dnux] = np.nan gdq = np.bitwise_or(bad_pixels_array[np.newaxis, np.newaxis, :, :], gdq) - fits.writeto("updategdq.fits". gdq, overwrite=True) + fits.writeto("updategdq.fits", gdq, overwrite=True) first_diffs = np.diff(data, axis=1) fits.writeto("first_diffs.fits", first_diffs, overwrite=True) all_ellipses = [] From a3cbd359562dc4245c6318ab7ae159cbf014303e Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 14:49:36 -0500 Subject: [PATCH 14/69] tests --- src/stcal/jump/jump.py | 15 +++++++++------ tests/test_jump.py | 7 ++++--- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 8c491ab6..e7068ade 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -847,11 +847,13 @@ def find_faint_extended( gdq = ingdq.copy() data = indata.copy() read_noise_2 = readnoise_2d**2 - bad_pixels_array = np.bitwise_and(pdq, 1) - fits.writeto("badpixelsarray.fits", bad_pixels_array, overwrite=True) - dnuy, dnux = np.where(bad_pixels_array == 1) - data[:, :, dnuy, dnux] = np.nan - gdq = np.bitwise_or(bad_pixels_array[np.newaxis, np.newaxis, :, :], gdq) +# ref_pixels_array = np.zeros_like(pdq) + ref_pixels_array = np.bitwise_and(pdq, 2147483648) // 2147483648 + test_array = ref_pixels_array // 2147483648 + fits.writeto("refpixelsarray.fits", ref_pixels_array, overwrite=True) + refy, refx = np.where(ref_pixels_array == 1) + data[:, :, refy, refx] = np.nan + gdq = np.bitwise_or(ref_pixels_array[np.newaxis, np.newaxis, :, :], gdq) fits.writeto("updategdq.fits", gdq, overwrite=True) first_diffs = np.diff(data, axis=1) fits.writeto("first_diffs.fits", first_diffs, overwrite=True) @@ -899,7 +901,8 @@ def find_faint_extended( sat_pixels_array = np.bitwise_and(combined_pixel_mask, 1) dnuy, dnux = np.where(sat_pixels_array == 1) masked_ratio[saty, satx] = np.nan - +# if grp == 1: +# fits.writeto("masked_ratio1.fits", masked_ratio., overwrite=True) masked_smoothed_ratio = convolve(masked_ratio, ring_2D_kernel) nrows = ratio.shape[1] ncols = ratio.shape[2] diff --git a/tests/test_jump.py b/tests/test_jump.py index 9f177da5..ec47d3a7 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -211,11 +211,12 @@ def test_flag_large_events_withsnowball_noextension(): def test_find_faint_extended(): - nint, ngrps, ncols, nrows = 1, 66, 30, 30 + nint, ngrps, ncols, nrows = 1, 66, 5, 5 data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) - gdq = np.zeros_like(data, dtype=np.uint8) - pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint8) + gdq = np.zeros_like(data, dtype=np.uint32) + pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) pdq[0, 0] = 1 + pdq[1, 1] = 2147483648 # pdq = np.zeros(shape=(data.shape[2], data.shape[3]), dtype=np.uint8) gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain From 3cf3a51b5bae1ef6d3da52a26f961d6a28684ab8 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 17 Jan 2024 15:21:09 -0500 Subject: [PATCH 15/69] Update jump.py --- src/stcal/jump/jump.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index e7068ade..6b3783aa 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -898,9 +898,9 @@ def find_faint_extended( masked_ratio[saty, satx] = np.nan # mask pixels that are already flagged as do not use - sat_pixels_array = np.bitwise_and(combined_pixel_mask, 1) - dnuy, dnux = np.where(sat_pixels_array == 1) - masked_ratio[saty, satx] = np.nan + dnu_pixels_array = np.bitwise_and(combined_pixel_mask, 1) + dnuy, dnux = np.where(dnu_pixels_array == 1) + masked_ratio[dnuy, dnux] = np.nan # if grp == 1: # fits.writeto("masked_ratio1.fits", masked_ratio., overwrite=True) masked_smoothed_ratio = convolve(masked_ratio, ring_2D_kernel) From 39b9962cfb06233ded3f4b2a52edff0a2db30e91 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 18 Jan 2024 10:17:44 -0500 Subject: [PATCH 16/69] Update jump.py --- src/stcal/jump/jump.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 6b3783aa..ad25f608 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -446,6 +446,9 @@ def detect_jumps( ) log.info("Total snowballs = %i", total_snowballs) number_extended_events = total_snowballs + fits.writeto("incoming_data.fits", data, overwrite=True) + fits.writeto("incoming_pdq.fits", pdq, overwrite=True) + fits.writeto("incoming_gdq.fits", gdq, overwrite=True) if find_showers: gdq, num_showers = find_faint_extended( data, From 3df80013a4c5b0cabf94666220454bc2f7bf0d99 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 18 Jan 2024 11:26:54 -0500 Subject: [PATCH 17/69] fixes --- src/stcal/jump/jump.py | 32 ++++++++++++++++-------- tests/test_jump.py | 56 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 78 insertions(+), 10 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index ad25f608..501562b3 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -218,6 +218,7 @@ def detect_jumps( constants.update_dqflags(dqflags) # populate dq flags sat_flag = dqflags["SATURATED"] jump_flag = dqflags["JUMP_DET"] + refpix_flag = dqflags["REFERENCE_PIXEL"] number_extended_events = 0 # Flag the pixeldq where the gain is <=0 or NaN so they will be ignored wh_g = np.where(gain_2d <= 0.0) @@ -309,6 +310,7 @@ def detect_jumps( outer=extend_outer_radius, sat_flag=sat_flag, jump_flag=jump_flag, + refpix_flag=refpix_flag, ellipse_expand=extend_ellipse_expand_ratio, num_grps_masked=grps_masked_after_shower, max_extended_radius=max_extended_radius, @@ -463,6 +465,7 @@ def detect_jumps( outer=extend_outer_radius, sat_flag=sat_flag, jump_flag=jump_flag, + refpix_flag=refpix_flag, ellipse_expand=extend_ellipse_expand_ratio, num_grps_masked=grps_masked_after_shower, max_extended_radius=max_extended_radius, @@ -799,8 +802,10 @@ def find_faint_extended( min_shower_area=40, inner=1, outer=2, + donotuse_flag = 1, sat_flag=2, jump_flag=4, + refpix_flag=2147483648, ellipse_expand=1.1, num_grps_masked=25, max_extended_radius=200, @@ -845,21 +850,24 @@ def find_faint_extended( Total number of showers detected. """ - fits.writeto("ingdq.fits", ingdq, overwrite=True) - fits.writeto("inpdg.fits", pdq, overwrite=True) +# fits.writeto("ingdq.fits", ingdq, overwrite=True) +# fits.writeto("inpdg.fits", pdq, overwrite=True) gdq = ingdq.copy() data = indata.copy() read_noise_2 = readnoise_2d**2 # ref_pixels_array = np.zeros_like(pdq) - ref_pixels_array = np.bitwise_and(pdq, 2147483648) // 2147483648 - test_array = ref_pixels_array // 2147483648 - fits.writeto("refpixelsarray.fits", ref_pixels_array, overwrite=True) - refy, refx = np.where(ref_pixels_array == 1) - data[:, :, refy, refx] = np.nan - gdq = np.bitwise_or(ref_pixels_array[np.newaxis, np.newaxis, :, :], gdq) +# ref_pixels_array = np.bitwise_and(pdq, 2147483648) // 2147483648 +# fits.writeto("refpixelsarray.fits", ref_pixels_array, overwrite=True) +# refy, refx = np.where(ref_pixels_array == 1) +# data[:, :, refy, refx] = np.nan + gdq = np.bitwise_or(pdq[np.newaxis, np.newaxis, :, :], gdq) + refint, refgrp, refy, refx = np.where(gdq == refpix_flag) + gdq[refint, refgrp, refy, refx] = donotuse_flag fits.writeto("updategdq.fits", gdq, overwrite=True) + x = gdq[0, 2, 1, 1] first_diffs = np.diff(data, axis=1) fits.writeto("first_diffs.fits", first_diffs, overwrite=True) + all_ellipses = [] first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) @@ -904,9 +912,13 @@ def find_faint_extended( dnu_pixels_array = np.bitwise_and(combined_pixel_mask, 1) dnuy, dnux = np.where(dnu_pixels_array == 1) masked_ratio[dnuy, dnux] = np.nan -# if grp == 1: -# fits.writeto("masked_ratio1.fits", masked_ratio., overwrite=True) + if grp == 2: + fits.writeto("dnu_pixel_array.fits", dnu_pixels_array, overwrite=True) + fits.writeto("masked_ratio1.fits", masked_ratio.filled(np.nan), overwrite=True) masked_smoothed_ratio = convolve(masked_ratio, ring_2D_kernel) + # mask out the pixels that got refilled by the convolution + masked_smoothed_ratio[dnuy, dnux] = np.nan + nrows = ratio.shape[1] ncols = ratio.shape[2] extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) diff --git a/tests/test_jump.py b/tests/test_jump.py index ec47d3a7..e4c4f032 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -275,6 +275,62 @@ def test_find_faint_extended(): # Check that the flags are not applied in the 3rd group after the event assert np.all(gdq[0, 4, 12:22, 14:23]) == 0 +def test_shower_enhancements(): + data = fits.getdata("data/incoming_data.fits") + pdq = fits.getdata("data/incoming_pdq.fits") + gdq = fits.getdata("data/incoming_gdq.fits") + readnoise_DN = fits.getdata("/grp/crds/jwst/references/jwst/jwst_miri_readnoise_0085.fits") + gain = fits.getdata("/grp/crds/jwst/references/jwst/jwst_miri_gain_0020.fits") + readnoise = gain * readnoise_DN + gdq, num_showers = find_faint_extended( + data, + gdq, + pdq, + readnoise * np.sqrt(2), + 1, + 100, + snr_threshold=3, + min_shower_area=10, + inner=1, + outer=2.6, + sat_flag=2, + jump_flag=4, + ellipse_expand=1.1, + num_grps_masked=0, + ) + fits.writeto("outputgdq.fits", gdq, overwrite=True) + def test_find_faint_extended(): + nint, ngrps, ncols, nrows = 1, 66, 5, 5 + data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) + gdq = np.zeros_like(data, dtype=np.uint32) + pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) + pdq[0, 0] = 1 + pdq[1, 1] = 2147483648 + # pdq = np.zeros(shape=(data.shape[2], data.shape[3]), dtype=np.uint8) + gain = 4 + readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain + rng = np.random.default_rng(12345) + data[0, 1:, 14:20, 15:20] = 6 * gain * 6.0 * np.sqrt(2) + data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise + fits.writeto("data.fits", data, overwrite=True) + gdq, num_showers = find_faint_extended( + data, + gdq, + pdq, + readnoise * np.sqrt(2), + 1, + 100, + snr_threshold=3, + min_shower_area=10, + inner=1, + outer=2.6, + sat_flag=2, + jump_flag=4, + ellipse_expand=1.1, + num_grps_masked=0, + ) + fits.writeto("outputgdq.fits", gdq, overwrite=True) + # No shower is found because the event is identical in all ints def test_find_faint_extended_sigclip(): From 7d9571601598a19cf00248a62fbc6855952ce47d Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 18 Jan 2024 15:00:17 -0500 Subject: [PATCH 18/69] updates --- src/stcal/jump/jump.py | 8 +++++++- src/stcal/jump/twopoint_difference.py | 3 ++- tests/test_jump.py | 2 +- 3 files changed, 10 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 501562b3..4ac755c8 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -1,6 +1,7 @@ import logging import multiprocessing import time +import warnings from astropy.io import fits import astropy.stats import cv2 as cv @@ -874,6 +875,8 @@ def find_faint_extended( nints = data.shape[0] if nints > minimum_sigclip_groups: mean, median, stddev = stats.sigma_clipped_stats(first_diffs_masked, sigma=5, axis=0) +# warnings.filterwarnings("ignore", ".*Input data contains invalid values.*", RuntimeWarning) + warnings.filterwarnings("ignore", ".*invalid value.*", RuntimeWarning) for intg in range(nints): # calculate sigma for each pixel if nints <= minimum_sigclip_groups: @@ -915,10 +918,11 @@ def find_faint_extended( if grp == 2: fits.writeto("dnu_pixel_array.fits", dnu_pixels_array, overwrite=True) fits.writeto("masked_ratio1.fits", masked_ratio.filled(np.nan), overwrite=True) + masked_smoothed_ratio = convolve(masked_ratio, ring_2D_kernel) # mask out the pixels that got refilled by the convolution masked_smoothed_ratio[dnuy, dnux] = np.nan - + nrows = ratio.shape[1] ncols = ratio.shape[2] extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) @@ -984,6 +988,8 @@ def find_faint_extended( if len(ellipses) > 0: # add all the showers for this integration to the list all_ellipses.append([intg, grp, ellipses]) + # Reset the warnings filter to its original state + warnings.resetwarnings() if all_ellipses: # Now we actually do the flagging of the pixels inside showers. # This is deferred until all showers are detected. because the showers diff --git a/src/stcal/jump/twopoint_difference.py b/src/stcal/jump/twopoint_difference.py index 62d44b1a..11109b5a 100644 --- a/src/stcal/jump/twopoint_difference.py +++ b/src/stcal/jump/twopoint_difference.py @@ -401,7 +401,8 @@ def find_crs( # the transient seen after ramp jumps flag_e_threshold = [after_jump_flag_e1, after_jump_flag_e2] flag_groups = [after_jump_flag_n1, after_jump_flag_n2] - + print("elec thresholds", flag_e_threshold) + print("groups flagged", flag_groups) for cthres, cgroup in zip(flag_e_threshold, flag_groups): if cgroup > 0: cr_intg, cr_group, cr_row, cr_col = np.where(np.bitwise_and(gdq, jump_flag)) diff --git a/tests/test_jump.py b/tests/test_jump.py index e4c4f032..4fe30e4a 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -290,7 +290,7 @@ def test_shower_enhancements(): 1, 100, snr_threshold=3, - min_shower_area=10, + min_shower_area=50, inner=1, outer=2.6, sat_flag=2, From 0195ac080d8ef7f63d246aa8ddf27c7f25d8c6c9 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 18 Jan 2024 15:14:23 -0500 Subject: [PATCH 19/69] Update jump.py --- src/stcal/jump/jump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 4ac755c8..09a96fc8 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -237,6 +237,7 @@ def detect_jumps( data *= gain_2d err *= gain_2d readnoise_2d *= gain_2d + print("median gain", np.median(gain_2d)) # also apply to the after_jump thresholds after_jump_flag_e1 = after_jump_flag_dn1 * gain_2d after_jump_flag_e2 = after_jump_flag_dn2 * gain_2d From f14b1b3611c5fe891f18831ebb345c507c694d4e Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 19 Jan 2024 10:21:05 -0500 Subject: [PATCH 20/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 09a96fc8..61d4d696 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -877,7 +877,7 @@ def find_faint_extended( if nints > minimum_sigclip_groups: mean, median, stddev = stats.sigma_clipped_stats(first_diffs_masked, sigma=5, axis=0) # warnings.filterwarnings("ignore", ".*Input data contains invalid values.*", RuntimeWarning) - warnings.filterwarnings("ignore", ".*invalid value.*", RuntimeWarning) + warnings.filterwarnings("ignore", ".*nan_treatment.*", RuntimeWarning) for intg in range(nints): # calculate sigma for each pixel if nints <= minimum_sigclip_groups: From 75f939951acc59e65243e361e5ef5e4c0586c52d Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 19 Jan 2024 10:45:47 -0500 Subject: [PATCH 21/69] Update twopoint_difference.py --- src/stcal/jump/twopoint_difference.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/stcal/jump/twopoint_difference.py b/src/stcal/jump/twopoint_difference.py index 11109b5a..b6d92b59 100644 --- a/src/stcal/jump/twopoint_difference.py +++ b/src/stcal/jump/twopoint_difference.py @@ -401,8 +401,6 @@ def find_crs( # the transient seen after ramp jumps flag_e_threshold = [after_jump_flag_e1, after_jump_flag_e2] flag_groups = [after_jump_flag_n1, after_jump_flag_n2] - print("elec thresholds", flag_e_threshold) - print("groups flagged", flag_groups) for cthres, cgroup in zip(flag_e_threshold, flag_groups): if cgroup > 0: cr_intg, cr_group, cr_row, cr_col = np.where(np.bitwise_and(gdq, jump_flag)) From f696a8c441200a895794959ded1438371b56c84e Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 19 Jan 2024 13:30:22 -0500 Subject: [PATCH 22/69] Update twopoint_difference.py --- src/stcal/jump/twopoint_difference.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/twopoint_difference.py b/src/stcal/jump/twopoint_difference.py index b6d92b59..a45b68a0 100644 --- a/src/stcal/jump/twopoint_difference.py +++ b/src/stcal/jump/twopoint_difference.py @@ -401,6 +401,7 @@ def find_crs( # the transient seen after ramp jumps flag_e_threshold = [after_jump_flag_e1, after_jump_flag_e2] flag_groups = [after_jump_flag_n1, after_jump_flag_n2] + print("flag groups", flag_groups) for cthres, cgroup in zip(flag_e_threshold, flag_groups): if cgroup > 0: cr_intg, cr_group, cr_row, cr_col = np.where(np.bitwise_and(gdq, jump_flag)) From 0c18abf4320da76da379a306985d948c4d137c9f Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 19 Jan 2024 13:33:50 -0500 Subject: [PATCH 23/69] Update twopoint_difference.py --- src/stcal/jump/twopoint_difference.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/twopoint_difference.py b/src/stcal/jump/twopoint_difference.py index a45b68a0..955c1c43 100644 --- a/src/stcal/jump/twopoint_difference.py +++ b/src/stcal/jump/twopoint_difference.py @@ -123,6 +123,7 @@ def find_crs( pixels above current row also to be flagged as a CR """ + print("n1, n2", after_jump_flag_n1, after_jump_flag_n2) # copy data and group DQ array if copy_arrs: dat = dataa.copy() From e2dfd702cbf6659d21817c9edbced3992ace8111 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 19 Jan 2024 14:31:36 -0500 Subject: [PATCH 24/69] Update twopoint_difference.py --- src/stcal/jump/twopoint_difference.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/twopoint_difference.py b/src/stcal/jump/twopoint_difference.py index 955c1c43..184a9a19 100644 --- a/src/stcal/jump/twopoint_difference.py +++ b/src/stcal/jump/twopoint_difference.py @@ -124,6 +124,7 @@ def find_crs( """ print("n1, n2", after_jump_flag_n1, after_jump_flag_n2) + print("e1, e2", after_jump_flag_e1, after_jump_flag_e2) # copy data and group DQ array if copy_arrs: dat = dataa.copy() From 5a6b120c892657c2e2616ad554ea49df6e36ce4f Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 19 Jan 2024 15:57:53 -0500 Subject: [PATCH 25/69] skip flagged groups --- src/stcal/jump/jump.py | 14 +++++++++++++- src/stcal/jump/twopoint_difference.py | 2 -- tests/test_jump.py | 14 ++++++++++++++ 3 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 61d4d696..4b2bdd59 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -891,7 +891,8 @@ def find_faint_extended( # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] - for grp in range(1, ngrps): + first_good_group = find_first_good_group(gdq[intg, :, :, :]) + for grp in range(first_good_group, ngrps): if nints > minimum_sigclip_groups: median_diffs = median[grp - 1] sigma = stddev[grp - 1] @@ -1014,6 +1015,17 @@ def find_faint_extended( ) return gdq, len(all_ellipses) +def find_first_good_group(int_gdq, do_not_use): + ngrps = int_gdq.shape[0] + skip_grp = True + first_good_group = 0 + for grp in range(ngrps): + mask = np.logical_and(int_gdq[grp], do_not_use) + skip_grp = np.unique(mask) + if not skip_grp: + first_good_group = grp + break + return first_good_group def calc_num_slices(n_rows, max_cores, max_available): n_slices = 1 diff --git a/src/stcal/jump/twopoint_difference.py b/src/stcal/jump/twopoint_difference.py index 184a9a19..a45b68a0 100644 --- a/src/stcal/jump/twopoint_difference.py +++ b/src/stcal/jump/twopoint_difference.py @@ -123,8 +123,6 @@ def find_crs( pixels above current row also to be flagged as a CR """ - print("n1, n2", after_jump_flag_n1, after_jump_flag_n2) - print("e1, e2", after_jump_flag_e1, after_jump_flag_e2) # copy data and group DQ array if copy_arrs: dat = dataa.copy() diff --git a/tests/test_jump.py b/tests/test_jump.py index 4fe30e4a..20a4231f 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -8,6 +8,7 @@ find_faint_extended, flag_large_events, point_inside_ellipse, + find_first_good_group, ) DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8} @@ -417,3 +418,16 @@ def test_calc_num_slices(): assert calc_num_slices(n_rows, "3/4", max_available_cores) == 1 n_rows = 9 assert calc_num_slices(n_rows, "21", max_available_cores) == 9 + +def test_find_first_good_grp(): + ngrps = 5 + ncols = 2 + nrows = 2 + intg_gdq = np.zeros(shape=(ngrps, ncols, nrows), dtype=np.uint32) + assert find_first_good_group(intg_gdq, DQFLAGS['DO_NOT_USE']) == 0 + intg_gdq[0, :, :] = 5 + assert find_first_good_group(intg_gdq, DQFLAGS['DO_NOT_USE']) == 1 + intg_gdq[1, :, :] = 5 + assert find_first_good_group(intg_gdq, DQFLAGS['DO_NOT_USE']) == 2 + intg_gdq[0, 0, 1] = 4 + assert find_first_good_group(intg_gdq, DQFLAGS['DO_NOT_USE']) == 0 From 7d6280a408bd9c3d486275ec73874cb137739471 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 19 Jan 2024 16:06:47 -0500 Subject: [PATCH 26/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 4b2bdd59..86b3d37c 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -891,7 +891,7 @@ def find_faint_extended( # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] - first_good_group = find_first_good_group(gdq[intg, :, :, :]) + first_good_group = find_first_good_group(gdq[intg, :, :, :], donotuse_flag) for grp in range(first_good_group, ngrps): if nints > minimum_sigclip_groups: median_diffs = median[grp - 1] From 39b03bc766363c26b0d31e42b47c9c06f7a4054a Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 19 Jan 2024 16:15:30 -0500 Subject: [PATCH 27/69] Update jump.py --- src/stcal/jump/jump.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 86b3d37c..1fcceb24 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -1020,8 +1020,8 @@ def find_first_good_group(int_gdq, do_not_use): skip_grp = True first_good_group = 0 for grp in range(ngrps): - mask = np.logical_and(int_gdq[grp], do_not_use) - skip_grp = np.unique(mask) + mask = np.bitwise_and(int_gdq[grp], do_not_use) + skip_grp = np.all(mask) if not skip_grp: first_good_group = grp break From dd9c50176ea7a90fda02722274daec53f8c56130 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 19 Jan 2024 16:20:47 -0500 Subject: [PATCH 28/69] fixes --- src/stcal/jump/jump.py | 2 +- src/stcal/jump/twopoint_difference.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 1fcceb24..56356820 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -892,7 +892,7 @@ def find_faint_extended( ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] first_good_group = find_first_good_group(gdq[intg, :, :, :], donotuse_flag) - for grp in range(first_good_group, ngrps): + for grp in range(first_good_group + 1, ngrps): if nints > minimum_sigclip_groups: median_diffs = median[grp - 1] sigma = stddev[grp - 1] diff --git a/src/stcal/jump/twopoint_difference.py b/src/stcal/jump/twopoint_difference.py index a45b68a0..b6d92b59 100644 --- a/src/stcal/jump/twopoint_difference.py +++ b/src/stcal/jump/twopoint_difference.py @@ -401,7 +401,6 @@ def find_crs( # the transient seen after ramp jumps flag_e_threshold = [after_jump_flag_e1, after_jump_flag_e2] flag_groups = [after_jump_flag_n1, after_jump_flag_n2] - print("flag groups", flag_groups) for cthres, cgroup in zip(flag_e_threshold, flag_groups): if cgroup > 0: cr_intg, cr_group, cr_row, cr_col = np.where(np.bitwise_and(gdq, jump_flag)) From ac35fda387dd2b2747e521baf61a48fa79562345 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Mon, 22 Jan 2024 13:12:21 -0500 Subject: [PATCH 29/69] updates --- src/stcal/jump/jump.py | 26 +++++++++++++------------- tests/test_jump.py | 8 +++++--- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 56356820..00676a9c 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -862,13 +862,15 @@ def find_faint_extended( # fits.writeto("refpixelsarray.fits", ref_pixels_array, overwrite=True) # refy, refx = np.where(ref_pixels_array == 1) # data[:, :, refy, refx] = np.nan - gdq = np.bitwise_or(pdq[np.newaxis, np.newaxis, :, :], gdq) +# gdq = np.bitwise_or(pdq[np.newaxis, np.newaxis, :, :], gdq) refint, refgrp, refy, refx = np.where(gdq == refpix_flag) gdq[refint, refgrp, refy, refx] = donotuse_flag fits.writeto("updategdq.fits", gdq, overwrite=True) - x = gdq[0, 2, 1, 1] first_diffs = np.diff(data, axis=1) fits.writeto("first_diffs.fits", first_diffs, overwrite=True) + masked_ratio_cube = np.zeros_like(data) + masked_smoothed_ratio_cube = np.zeros_like(data) + extended_emission_cube = np.zeros_like(data) all_ellipses = [] @@ -917,14 +919,12 @@ def find_faint_extended( dnu_pixels_array = np.bitwise_and(combined_pixel_mask, 1) dnuy, dnux = np.where(dnu_pixels_array == 1) masked_ratio[dnuy, dnux] = np.nan - if grp == 2: - fits.writeto("dnu_pixel_array.fits", dnu_pixels_array, overwrite=True) - fits.writeto("masked_ratio1.fits", masked_ratio.filled(np.nan), overwrite=True) + masked_ratio_cube[intg, grp, :, :] = masked_ratio - masked_smoothed_ratio = convolve(masked_ratio, ring_2D_kernel) + masked_smoothed_ratio = convolve(masked_ratio.filled(np.nan), ring_2D_kernel) # mask out the pixels that got refilled by the convolution masked_smoothed_ratio[dnuy, dnux] = np.nan - + masked_smoothed_ratio_cube[intg, grp, :, :] = masked_smoothed_ratio nrows = ratio.shape[1] ncols = ratio.shape[2] extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) @@ -932,16 +932,13 @@ def find_faint_extended( cutoff = median2 + snr_threshold * stddev2 exty, extx = np.where(masked_smoothed_ratio > cutoff) extended_emission[exty, extx] = 1 - if grp == 1: - fits.writeto("masked_smoothed_ratio1.fits", masked_smoothed_ratio, overwrite=True) - fits.writeto("extendedemission1.fits", extended_emission, overwrite=True) - if grp == 2: - fits.writeto("masked_smoothed_ratio.fits", masked_smoothed_ratio, overwrite=True) - fits.writeto("extendedemission.fits", extended_emission, overwrite=True) + extended_emission_cube[intg, grp, :, :] = extended_emission # find the contours of the extended emission contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) # get the contours that are above the minimum size bigcontours = [con for con in contours if cv.contourArea(con) > min_shower_area] + if grp == 5: + a=1 # get the minimum enclosing rectangle which is the same as the # minimum enclosing ellipse ellipses = [cv.minAreaRect(con) for con in bigcontours] @@ -992,6 +989,9 @@ def find_faint_extended( all_ellipses.append([intg, grp, ellipses]) # Reset the warnings filter to its original state warnings.resetwarnings() + fits.writeto("masked_ratio_cube.fits", masked_ratio_cube, overwrite=True) + fits.writeto("masked_smoothed_ratio_cube.fits", masked_smoothed_ratio_cube, overwrite=True) + fits.writeto("extended_emission_cube.fits", extended_emission_cube, overwrite=True) if all_ellipses: # Now we actually do the flagging of the pixels inside showers. # This is deferred until all showers are detected. because the showers diff --git a/tests/test_jump.py b/tests/test_jump.py index 20a4231f..51fc2c81 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -282,6 +282,7 @@ def test_shower_enhancements(): gdq = fits.getdata("data/incoming_gdq.fits") readnoise_DN = fits.getdata("/grp/crds/jwst/references/jwst/jwst_miri_readnoise_0085.fits") gain = fits.getdata("/grp/crds/jwst/references/jwst/jwst_miri_gain_0020.fits") + print("median gain", np.nanmedian(gain)) readnoise = gain * readnoise_DN gdq, num_showers = find_faint_extended( data, @@ -290,8 +291,8 @@ def test_shower_enhancements(): readnoise * np.sqrt(2), 1, 100, - snr_threshold=3, - min_shower_area=50, + snr_threshold=3.0, + min_shower_area=30, inner=1, outer=2.6, sat_flag=2, @@ -344,6 +345,7 @@ def test_find_faint_extended_sigclip(): rng = np.random.default_rng(12345) data[0, 1:, 14:20, 15:20] = 6 * gain * 1.7 data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise + min_shower_area=20 gdq, num_showers = find_faint_extended( data, gdq, @@ -352,7 +354,7 @@ def test_find_faint_extended_sigclip(): 1, 100, snr_threshold=1.3, - min_shower_area=20, + min_shower_area=min_shower_area, inner=1, outer=2, sat_flag=2, From 36ed8e1eed2465e853248b8404818c6046e6afb2 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Tue, 23 Jan 2024 15:31:50 -0500 Subject: [PATCH 30/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 00676a9c..2f9e49c1 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -879,7 +879,7 @@ def find_faint_extended( if nints > minimum_sigclip_groups: mean, median, stddev = stats.sigma_clipped_stats(first_diffs_masked, sigma=5, axis=0) # warnings.filterwarnings("ignore", ".*Input data contains invalid values.*", RuntimeWarning) - warnings.filterwarnings("ignore", ".*nan_treatment.*", RuntimeWarning) + warnings.filterwarnings("ignore") for intg in range(nints): # calculate sigma for each pixel if nints <= minimum_sigclip_groups: From c39e85fea88ccc7a0136c816b3a343bb1ff545ea Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 24 Jan 2024 10:13:28 -0500 Subject: [PATCH 31/69] Update jump.py --- src/stcal/jump/jump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 2f9e49c1..207b589c 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -1013,6 +1013,7 @@ def find_faint_extended( num_grps_masked=num_grps_masked, max_extended_radius=max_extended_radius, ) + fits.writeto("aftergdq.fits", gdq, overwrite=True) return gdq, len(all_ellipses) def find_first_good_group(int_gdq, do_not_use): From b7c67b594693155784e50af0ff35d78abc0b2905 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 24 Jan 2024 10:39:22 -0500 Subject: [PATCH 32/69] Update jump.py --- src/stcal/jump/jump.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 207b589c..e60855d9 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -542,7 +542,7 @@ def flag_large_events( Nothing, gdq array is modified. """ - log.info("Flagging large Snowballs") + log.info("Flagging Snowballs") n_showers_grp = [] total_snowballs = 0 @@ -854,6 +854,7 @@ def find_faint_extended( """ # fits.writeto("ingdq.fits", ingdq, overwrite=True) # fits.writeto("inpdg.fits", pdq, overwrite=True) + log.info("Flagging Showers") gdq = ingdq.copy() data = indata.copy() read_noise_2 = readnoise_2d**2 @@ -863,8 +864,8 @@ def find_faint_extended( # refy, refx = np.where(ref_pixels_array == 1) # data[:, :, refy, refx] = np.nan # gdq = np.bitwise_or(pdq[np.newaxis, np.newaxis, :, :], gdq) - refint, refgrp, refy, refx = np.where(gdq == refpix_flag) - gdq[refint, refgrp, refy, refx] = donotuse_flag + refy, refx = np.where(pdq == refpix_flag) + gdq[:, :, refy, refx] = donotuse_flag fits.writeto("updategdq.fits", gdq, overwrite=True) first_diffs = np.diff(data, axis=1) fits.writeto("first_diffs.fits", first_diffs, overwrite=True) From f96c8bde82feedd7ed77a69840a4e2a9fea36da5 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 24 Jan 2024 15:30:56 -0500 Subject: [PATCH 33/69] tests --- src/stcal/jump/jump.py | 2 +- tests/test_jump.py | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index e60855d9..7da45c00 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -320,7 +320,7 @@ def detect_jumps( log.info("Total showers= %i", num_showers) number_extended_events = num_showers else: - yinc = int(n_rows / n_slices) + yinc = int(n_rows // n_slices) slices = [] # Slice up data, gdq, readnoise_2d into slices # Each element of slices is a tuple of diff --git a/tests/test_jump.py b/tests/test_jump.py index 51fc2c81..fa1a4434 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -160,6 +160,7 @@ def test_flag_large_events_groupedsnowball(): cube[0, 2, 5, 1:6] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 1] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 5] = DQFLAGS["JUMP_DET"] + fits.writeto("incube.fits", cube, overwrite=True) flag_large_events( cube, DQFLAGS["JUMP_DET"], @@ -174,6 +175,7 @@ def test_flag_large_events_groupedsnowball(): ) # assert cube[0, 1, 2, 2] == 0 # assert cube[0, 1, 3, 5] == 0 + fits.writeto("outcube.fits", cube, overwrite=True) assert cube[0, 2, 0, 0] == 0 assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended From dfb793e33e58702c98e9731d8e6f2fb4a27547e4 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 24 Jan 2024 16:54:37 -0500 Subject: [PATCH 34/69] Update jump.py --- src/stcal/jump/jump.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 7da45c00..88662b40 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -433,6 +433,12 @@ def detect_jumps( # save the neighbors to be flagged that will be in the next slice previous_row_above_gdq = row_above_gdq.copy() k += 1 + # remove redundant flagging of pixels that have jump flagged but were + # already flagged as do_not_use or saturated. + gdq[gdq == np.logical_or(dqflags['DO_NOT_USE'], dqflags['JUMP_DET'])] = \ + dqflags['DO_NOT_USE'] + gdq[gdq == np.logical_or(dqflags['SATURATED'], dqflags['JUMP_DET'])] = \ + dqflags['SATURATED'] # This is the flag that controls the flagging of snowballs. if expand_large_events: total_snowballs = flag_large_events( From d728225bdd35c2c9d17d90e2a5c52d79992c2a18 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 25 Jan 2024 08:15:56 -0500 Subject: [PATCH 35/69] Update jump.py --- src/stcal/jump/jump.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 88662b40..3b2b5070 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -435,9 +435,13 @@ def detect_jumps( k += 1 # remove redundant flagging of pixels that have jump flagged but were # already flagged as do_not_use or saturated. - gdq[gdq == np.logical_or(dqflags['DO_NOT_USE'], dqflags['JUMP_DET'])] = \ + # ints, grps, rows, cols = np.where(np.logical_or(dqflags['DO_NOT_USE'], dqflags['JUMP_DET'])) + # gdq[ints, grps, rows, cols] = dqflags['DO_NOT_USE'] + # ints, grps, rows, cols = np.where(np.logical_or(dqflags['SATURATED'], dqflags['JUMP_DET'])) + # gdq[ints, grps, rows, cols] = dqflags['SATURATED'] + gdq[gdq == np.bitwise_or(dqflags['DO_NOT_USE'], dqflags['JUMP_DET'])] = \ dqflags['DO_NOT_USE'] - gdq[gdq == np.logical_or(dqflags['SATURATED'], dqflags['JUMP_DET'])] = \ + gdq[gdq == np.bitwise_or(dqflags['SATURATED'], dqflags['JUMP_DET'])] = \ dqflags['SATURATED'] # This is the flag that controls the flagging of snowballs. if expand_large_events: From d66d64e110a5452a1bc5bf51a809a844a843c89f Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 25 Jan 2024 11:00:00 -0500 Subject: [PATCH 36/69] Update jump.py --- src/stcal/jump/jump.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 3b2b5070..f360ee20 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -433,16 +433,13 @@ def detect_jumps( # save the neighbors to be flagged that will be in the next slice previous_row_above_gdq = row_above_gdq.copy() k += 1 - # remove redundant flagging of pixels that have jump flagged but were + # remove redundant bits in pixels that have jump flagged but were # already flagged as do_not_use or saturated. - # ints, grps, rows, cols = np.where(np.logical_or(dqflags['DO_NOT_USE'], dqflags['JUMP_DET'])) - # gdq[ints, grps, rows, cols] = dqflags['DO_NOT_USE'] - # ints, grps, rows, cols = np.where(np.logical_or(dqflags['SATURATED'], dqflags['JUMP_DET'])) - # gdq[ints, grps, rows, cols] = dqflags['SATURATED'] gdq[gdq == np.bitwise_or(dqflags['DO_NOT_USE'], dqflags['JUMP_DET'])] = \ dqflags['DO_NOT_USE'] gdq[gdq == np.bitwise_or(dqflags['SATURATED'], dqflags['JUMP_DET'])] = \ dqflags['SATURATED'] + # This is the flag that controls the flagging of snowballs. if expand_large_events: total_snowballs = flag_large_events( From 47028138bd6f9491d2585f8384e4863f6212bbc7 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 21 Feb 2024 14:27:23 -0500 Subject: [PATCH 37/69] Update jump.py --- src/stcal/jump/jump.py | 45 +++++++++++++++++------------------------- 1 file changed, 18 insertions(+), 27 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index f360ee20..31c194be 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -2,7 +2,6 @@ import multiprocessing import time import warnings -from astropy.io import fits import astropy.stats import cv2 as cv import numpy as np @@ -338,9 +337,9 @@ def detect_jumps( slices.insert( i, ( - data[:, :, i * yinc : (i + 1) * yinc, :], - gdq[:, :, i * yinc : (i + 1) * yinc, :], - readnoise_2d[i * yinc : (i + 1) * yinc, :], + data[:, :, i * yinc: (i + 1) * yinc, :], + gdq[:, :, i * yinc: (i + 1) * yinc, :], + readnoise_2d[i * yinc: (i + 1) * yinc, :], rejection_thresh, three_grp_thresh, four_grp_thresh, @@ -364,9 +363,9 @@ def detect_jumps( slices.insert( n_slices - 1, ( - data[:, :, (n_slices - 1) * yinc : n_rows, :], - gdq[:, :, (n_slices - 1) * yinc : n_rows, :], - readnoise_2d[(n_slices - 1) * yinc : n_rows, :], + data[:, :, (n_slices - 1) * yinc: n_rows, :], + gdq[:, :, (n_slices - 1) * yinc: n_rows, :], + readnoise_2d[(n_slices - 1) * yinc: n_rows, :], rejection_thresh, three_grp_thresh, four_grp_thresh, @@ -406,15 +405,15 @@ def detect_jumps( stddev = np.zeros((nrows, ncols), dtype=np.float32) for resultslice in real_result: if len(real_result) == k + 1: # last result - gdq[:, :, k * yinc : n_rows, :] = resultslice[0] + gdq[:, :, k * yinc: n_rows, :] = resultslice[0] if only_use_ints: - stddev[:, k * yinc : n_rows, :] = resultslice[4] + stddev[:, k * yinc: n_rows, :] = resultslice[4] else: - stddev[k * yinc : n_rows, :] = resultslice[4] + stddev[k * yinc: n_rows, :] = resultslice[4] else: - gdq[:, :, k * yinc : (k + 1) * yinc, :] = resultslice[0] + gdq[:, :, k * yinc: (k + 1) * yinc, :] = resultslice[0] if only_use_ints: - stddev[:, k * yinc : (k + 1) * yinc, :] = resultslice[4] + stddev[:, k * yinc: (k + 1) * yinc, :] = resultslice[4] else: stddev[k * yinc : (k + 1) * yinc, :] = resultslice[4] row_below_gdq[:, :, :] = resultslice[1] @@ -457,9 +456,6 @@ def detect_jumps( ) log.info("Total snowballs = %i", total_snowballs) number_extended_events = total_snowballs - fits.writeto("incoming_data.fits", data, overwrite=True) - fits.writeto("incoming_pdq.fits", pdq, overwrite=True) - fits.writeto("incoming_gdq.fits", gdq, overwrite=True) if find_showers: gdq, num_showers = find_faint_extended( data, @@ -865,17 +861,14 @@ def find_faint_extended( gdq = ingdq.copy() data = indata.copy() read_noise_2 = readnoise_2d**2 -# ref_pixels_array = np.zeros_like(pdq) -# ref_pixels_array = np.bitwise_and(pdq, 2147483648) // 2147483648 -# fits.writeto("refpixelsarray.fits", ref_pixels_array, overwrite=True) -# refy, refx = np.where(ref_pixels_array == 1) -# data[:, :, refy, refx] = np.nan -# gdq = np.bitwise_or(pdq[np.newaxis, np.newaxis, :, :], gdq) + gdq[gdq == np.bitwise_or(donotuse_flag, jump_flag)] = donotuse_flag + gdq[gdq == np.bitwise_or(donotuse_flag, sat_flag)] = donotuse_flag + data[gdq == sat_flag] = np.nan + data[gdq == jump_flag] = np.nan + data[gdq == donotuse_flag] = np.nan refy, refx = np.where(pdq == refpix_flag) gdq[:, :, refy, refx] = donotuse_flag - fits.writeto("updategdq.fits", gdq, overwrite=True) first_diffs = np.diff(data, axis=1) - fits.writeto("first_diffs.fits", first_diffs, overwrite=True) masked_ratio_cube = np.zeros_like(data) masked_smoothed_ratio_cube = np.zeros_like(data) extended_emission_cube = np.zeros_like(data) @@ -997,9 +990,7 @@ def find_faint_extended( all_ellipses.append([intg, grp, ellipses]) # Reset the warnings filter to its original state warnings.resetwarnings() - fits.writeto("masked_ratio_cube.fits", masked_ratio_cube, overwrite=True) - fits.writeto("masked_smoothed_ratio_cube.fits", masked_smoothed_ratio_cube, overwrite=True) - fits.writeto("extended_emission_cube.fits", extended_emission_cube, overwrite=True) + if all_ellipses: # Now we actually do the flagging of the pixels inside showers. # This is deferred until all showers are detected. because the showers @@ -1021,7 +1012,7 @@ def find_faint_extended( num_grps_masked=num_grps_masked, max_extended_radius=max_extended_radius, ) - fits.writeto("aftergdq.fits", gdq, overwrite=True) + return gdq, len(all_ellipses) def find_first_good_group(int_gdq, do_not_use): From 26a63960f743b4fcca3a63c7e98dff3e3c2bdac7 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 15:10:55 -0500 Subject: [PATCH 38/69] Update jump.py --- src/stcal/jump/jump.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 31c194be..8b20b517 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -3,6 +3,7 @@ import time import warnings import astropy.stats +from astropy.io import fits import cv2 as cv import numpy as np from astropy import stats @@ -989,6 +990,7 @@ def find_faint_extended( # add all the showers for this integration to the list all_ellipses.append([intg, grp, ellipses]) # Reset the warnings filter to its original state + fits.writeto("masked_smoothed_ratio_cube.fits",masked_smoothed_ratio_cube, overwrite=True) warnings.resetwarnings() if all_ellipses: From e3dd0aeb70bc9374503e19d3849e112868bd7693 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 15:22:14 -0500 Subject: [PATCH 39/69] Update jump.py --- src/stcal/jump/jump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 8b20b517..9b4a5165 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -990,6 +990,7 @@ def find_faint_extended( # add all the showers for this integration to the list all_ellipses.append([intg, grp, ellipses]) # Reset the warnings filter to its original state + print("writing msr cube") fits.writeto("masked_smoothed_ratio_cube.fits",masked_smoothed_ratio_cube, overwrite=True) warnings.resetwarnings() From 31019a2126a2df576ce7f6a1aaf31de2ddde152b Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 15:33:05 -0500 Subject: [PATCH 40/69] Update jump.py --- src/stcal/jump/jump.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 9b4a5165..13767f2a 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -873,7 +873,8 @@ def find_faint_extended( masked_ratio_cube = np.zeros_like(data) masked_smoothed_ratio_cube = np.zeros_like(data) extended_emission_cube = np.zeros_like(data) - + e_jump_cube = np.zeros_like(data) + median_cube = np.zeros_like(shape=(data.shape[0], data.shape[2], data.shape[3])) all_ellipses = [] first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) @@ -889,6 +890,8 @@ def find_faint_extended( sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # The difference from the median difference for each group e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] + e_jump_cube[intg, :, :, :] = e_jump + median_cube[intg, :, :] = median_diffs # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] @@ -991,6 +994,8 @@ def find_faint_extended( all_ellipses.append([intg, grp, ellipses]) # Reset the warnings filter to its original state print("writing msr cube") + fits.writeto("median_cube.fits", median_cube, overwrite=True) + fits.writeto("e_jump_cube.fits", e_jump_cube, overwrite=True) fits.writeto("masked_smoothed_ratio_cube.fits",masked_smoothed_ratio_cube, overwrite=True) warnings.resetwarnings() From 7f622dd293f4191ea000699fa5c8c968c32f1197 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 15:35:43 -0500 Subject: [PATCH 41/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 13767f2a..035c17fa 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -874,7 +874,7 @@ def find_faint_extended( masked_smoothed_ratio_cube = np.zeros_like(data) extended_emission_cube = np.zeros_like(data) e_jump_cube = np.zeros_like(data) - median_cube = np.zeros_like(shape=(data.shape[0], data.shape[2], data.shape[3])) + median_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) all_ellipses = [] first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) From 55ea0b0246257878335b87f3c518da398924c5d3 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 15:38:28 -0500 Subject: [PATCH 42/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 035c17fa..509e8d7c 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -873,7 +873,7 @@ def find_faint_extended( masked_ratio_cube = np.zeros_like(data) masked_smoothed_ratio_cube = np.zeros_like(data) extended_emission_cube = np.zeros_like(data) - e_jump_cube = np.zeros_like(data) + e_jump_cube = np.zeros_like(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) median_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) all_ellipses = [] From f4816c257fa01b805208819e5012560e29542375 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 15:39:56 -0500 Subject: [PATCH 43/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 509e8d7c..d044d435 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -873,7 +873,7 @@ def find_faint_extended( masked_ratio_cube = np.zeros_like(data) masked_smoothed_ratio_cube = np.zeros_like(data) extended_emission_cube = np.zeros_like(data) - e_jump_cube = np.zeros_like(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) + e_jump_cube = np.zeros(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) median_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) all_ellipses = [] From a52a302d3a7b1e6d77b8680ae297672770d2f877 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 15:49:24 -0500 Subject: [PATCH 44/69] Update jump.py --- src/stcal/jump/jump.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index d044d435..fa39bcfa 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -869,12 +869,14 @@ def find_faint_extended( data[gdq == donotuse_flag] = np.nan refy, refx = np.where(pdq == refpix_flag) gdq[:, :, refy, refx] = donotuse_flag + fits.writeto("masked_gdq.fits", gdq, overwrite=True) first_diffs = np.diff(data, axis=1) masked_ratio_cube = np.zeros_like(data) masked_smoothed_ratio_cube = np.zeros_like(data) extended_emission_cube = np.zeros_like(data) e_jump_cube = np.zeros(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) median_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) + sigma_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) all_ellipses = [] first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) @@ -892,6 +894,7 @@ def find_faint_extended( e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] e_jump_cube[intg, :, :, :] = e_jump median_cube[intg, :, :] = median_diffs + sigma_cube[intg, :, :] = sigma # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] @@ -994,6 +997,7 @@ def find_faint_extended( all_ellipses.append([intg, grp, ellipses]) # Reset the warnings filter to its original state print("writing msr cube") + fits.writeto("sigma_cube.fits", sigma_cube, overwrite=True) fits.writeto("median_cube.fits", median_cube, overwrite=True) fits.writeto("e_jump_cube.fits", e_jump_cube, overwrite=True) fits.writeto("masked_smoothed_ratio_cube.fits",masked_smoothed_ratio_cube, overwrite=True) From e59ec6676c09d46bb979e597b14d110d27307c98 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 16:05:21 -0500 Subject: [PATCH 45/69] Update jump.py --- src/stcal/jump/jump.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index fa39bcfa..743d5da1 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -877,6 +877,7 @@ def find_faint_extended( e_jump_cube = np.zeros(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) median_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) sigma_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) + ratio_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) all_ellipses = [] first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) @@ -897,7 +898,7 @@ def find_faint_extended( sigma_cube[intg, :, :] = sigma # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] - + ratio_cube[intg, :, :] = ratio # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] @@ -998,6 +999,7 @@ def find_faint_extended( # Reset the warnings filter to its original state print("writing msr cube") fits.writeto("sigma_cube.fits", sigma_cube, overwrite=True) + fits.writeto("ratio_cube.fits", ratio_cube, overwrite=True) fits.writeto("median_cube.fits", median_cube, overwrite=True) fits.writeto("e_jump_cube.fits", e_jump_cube, overwrite=True) fits.writeto("masked_smoothed_ratio_cube.fits",masked_smoothed_ratio_cube, overwrite=True) From f83a2a759025dcd41fcac20e240ede862273d45c Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 16:08:55 -0500 Subject: [PATCH 46/69] Update jump.py --- src/stcal/jump/jump.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 743d5da1..48844627 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -888,7 +888,7 @@ def find_faint_extended( warnings.filterwarnings("ignore") for intg in range(nints): # calculate sigma for each pixel - if nints <= minimum_sigclip_groups: + if nints < minimum_sigclip_groups: median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # The difference from the median difference for each group @@ -898,13 +898,13 @@ def find_faint_extended( sigma_cube[intg, :, :] = sigma # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] - ratio_cube[intg, :, :] = ratio + ratio_cube[intg, :, :, :] = ratio # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] first_good_group = find_first_good_group(gdq[intg, :, :, :], donotuse_flag) for grp in range(first_good_group + 1, ngrps): - if nints > minimum_sigclip_groups: + if nints >= minimum_sigclip_groups: median_diffs = median[grp - 1] sigma = stddev[grp - 1] # The difference from the median difference for each group @@ -946,8 +946,6 @@ def find_faint_extended( contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) # get the contours that are above the minimum size bigcontours = [con for con in contours if cv.contourArea(con) > min_shower_area] - if grp == 5: - a=1 # get the minimum enclosing rectangle which is the same as the # minimum enclosing ellipse ellipses = [cv.minAreaRect(con) for con in bigcontours] From b616564aa4b4843af4d89975be2b18d1adbe3526 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 16:10:51 -0500 Subject: [PATCH 47/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 48844627..307928e0 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -877,7 +877,7 @@ def find_faint_extended( e_jump_cube = np.zeros(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) median_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) sigma_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) - ratio_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) + ratio_cube = np.zeros(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) all_ellipses = [] first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) From 9a6e1fba74491dc6322df38fa33950e0c471d99c Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 16:18:34 -0500 Subject: [PATCH 48/69] Update jump.py --- src/stcal/jump/jump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 307928e0..ad610d20 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -996,6 +996,7 @@ def find_faint_extended( all_ellipses.append([intg, grp, ellipses]) # Reset the warnings filter to its original state print("writing msr cube") + fits.writeto("extended_emission_cube.fits", extended_emission_cube, overwrite=True) fits.writeto("sigma_cube.fits", sigma_cube, overwrite=True) fits.writeto("ratio_cube.fits", ratio_cube, overwrite=True) fits.writeto("median_cube.fits", median_cube, overwrite=True) From bcce2b7a2eca9c67d4d1c54988222044f0f28ed8 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 22 Feb 2024 16:28:05 -0500 Subject: [PATCH 49/69] Update jump.py --- src/stcal/jump/jump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index ad610d20..dc88e64c 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -939,6 +939,7 @@ def find_faint_extended( extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) mean2, median2, stddev2 = astropy.stats.sigma_clipped_stats(masked_smoothed_ratio) cutoff = median2 + snr_threshold * stddev2 + print("intg ", intg, "grp ", grp, "median2 ", median2, "stddev2 ", stddev2) exty, extx = np.where(masked_smoothed_ratio > cutoff) extended_emission[exty, extx] = 1 extended_emission_cube[intg, grp, :, :] = extended_emission From 0ce19399508455137dc9cb6140d97b6189256217 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 23 Feb 2024 10:01:04 -0500 Subject: [PATCH 50/69] Update jump.py --- src/stcal/jump/jump.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index dc88e64c..0812fd5c 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -937,9 +937,9 @@ def find_faint_extended( nrows = ratio.shape[1] ncols = ratio.shape[2] extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) - mean2, median2, stddev2 = astropy.stats.sigma_clipped_stats(masked_smoothed_ratio) - cutoff = median2 + snr_threshold * stddev2 - print("intg ", intg, "grp ", grp, "median2 ", median2, "stddev2 ", stddev2) +# mean2, median2, stddev2 = astropy.stats.sigma_clipped_stats(masked_smoothed_ratio) + cutoff = median_diffs + snr_threshold * sigma +# print("intg ", intg, "grp ", grp, "median2 ", median2, "stddev2 ", stddev2) exty, extx = np.where(masked_smoothed_ratio > cutoff) extended_emission[exty, extx] = 1 extended_emission_cube[intg, grp, :, :] = extended_emission From de23070d72a2fda55ea27424fea82036113515a1 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 23 Feb 2024 10:03:35 -0500 Subject: [PATCH 51/69] Update jump.py --- src/stcal/jump/jump.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 0812fd5c..af357032 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -938,9 +938,9 @@ def find_faint_extended( ncols = ratio.shape[2] extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) # mean2, median2, stddev2 = astropy.stats.sigma_clipped_stats(masked_smoothed_ratio) - cutoff = median_diffs + snr_threshold * sigma +# cutoff = median_diffs + snr_threshold * sigma # print("intg ", intg, "grp ", grp, "median2 ", median2, "stddev2 ", stddev2) - exty, extx = np.where(masked_smoothed_ratio > cutoff) + exty, extx = np.where(masked_smoothed_ratio > snr_threshold) extended_emission[exty, extx] = 1 extended_emission_cube[intg, grp, :, :] = extended_emission # find the contours of the extended emission From 9a0f7d4c8f4d44a805c5a78fd58bc2235504fd7b Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 23 Feb 2024 10:23:50 -0500 Subject: [PATCH 52/69] Update jump.py --- src/stcal/jump/jump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index af357032..b45cb4d0 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -912,6 +912,7 @@ def find_faint_extended( # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] masked_ratio = ratio[grp - 1].copy() + if ndiffs > # mask pixels that are already flagged as jump combined_pixel_mask = np.bitwise_or(gdq[intg, grp, :, :], pdq[:, :]) From cad45c9e815c7c5ff721c3cb1151316e83cfbf3f Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 23 Feb 2024 10:43:27 -0500 Subject: [PATCH 53/69] Update jump.py --- src/stcal/jump/jump.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index b45cb4d0..4074c415 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -884,13 +884,16 @@ def find_faint_extended( nints = data.shape[0] if nints > minimum_sigclip_groups: mean, median, stddev = stats.sigma_clipped_stats(first_diffs_masked, sigma=5, axis=0) + else: + median_diffs = np.nanmedian(first_diffs_masked, axis=(0, 1)) + sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # warnings.filterwarnings("ignore", ".*Input data contains invalid values.*", RuntimeWarning) warnings.filterwarnings("ignore") for intg in range(nints): # calculate sigma for each pixel if nints < minimum_sigclip_groups: - median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) - sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) +# median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) +# sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # The difference from the median difference for each group e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] e_jump_cube[intg, :, :, :] = e_jump @@ -902,6 +905,7 @@ def find_faint_extended( # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] + nints = data.shape[0] first_good_group = find_first_good_group(gdq[intg, :, :, :], donotuse_flag) for grp in range(first_good_group + 1, ngrps): if nints >= minimum_sigclip_groups: @@ -912,7 +916,6 @@ def find_faint_extended( # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] masked_ratio = ratio[grp - 1].copy() - if ndiffs > # mask pixels that are already flagged as jump combined_pixel_mask = np.bitwise_or(gdq[intg, grp, :, :], pdq[:, :]) From 02b4af9fd34c7d96b56a1ec8a74939e873bc0b60 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 23 Feb 2024 14:52:12 -0500 Subject: [PATCH 54/69] Update jump.py --- src/stcal/jump/jump.py | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 4074c415..1dfccd76 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -895,13 +895,25 @@ def find_faint_extended( # median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) # sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # The difference from the median difference for each group - e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] - e_jump_cube[intg, :, :, :] = e_jump - median_cube[intg, :, :] = median_diffs - sigma_cube[intg, :, :] = sigma - # SNR ratio of each diff. - ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] - ratio_cube[intg, :, :, :] = ratio + if intg > 0: + e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] + e_jump_cube[intg, :, :, :] = e_jump + median_cube[intg, :, :] = median_diffs + sigma_cube[intg, :, :] = sigma + # SNR ratio of each diff. + ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] + ratio_cube[intg, :, :, :] = ratio + else: + median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) + sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) + # The difference from the median difference for each group + e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] + # SNR ratio of each diff. + ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] + e_jump_cube[intg, :, :, :] = e_jump + median_cube[intg, :, :] = median_diffs + sigma_cube[intg, :, :] = sigma + ratio_cube[intg, :, :, :] = ratio # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] From 68d9c4dcb3291b18b4f4a001fbc6eeebf3c9395f Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 23 Feb 2024 15:02:54 -0500 Subject: [PATCH 55/69] Update jump.py --- src/stcal/jump/jump.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 1dfccd76..85258220 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -914,6 +914,8 @@ def find_faint_extended( median_cube[intg, :, :] = median_diffs sigma_cube[intg, :, :] = sigma ratio_cube[intg, :, :, :] = ratio + median_diffs = np.nanmedian(first_diffs_masked, axis=(0, 1)) + sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] From 45b63cec2ca03543436d46a668755a29daae34b7 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Mon, 26 Feb 2024 15:51:48 -0500 Subject: [PATCH 56/69] Update jump.py --- src/stcal/jump/jump.py | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 85258220..f2712bc4 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -815,6 +815,7 @@ def find_faint_extended( ellipse_expand=1.1, num_grps_masked=25, max_extended_radius=200, + min_diffs_for_shower=10, ): """ Parameters @@ -860,7 +861,19 @@ def find_faint_extended( # fits.writeto("inpdg.fits", pdq, overwrite=True) log.info("Flagging Showers") gdq = ingdq.copy() + fits.writeto("input_gdq.fits", gdq, overwrite=True) data = indata.copy() + nints = data.shape[0] + ngrps = data.shape[1] + num_grps_donotuse = 0 + for integ in range(nints): + for grp in range(ngrps): + if np.all(np.bitwise_and(gdq[integ, grp, :, :], donotuse_flag)): + num_grps_donotuse += 1 + total_diffs = nints * (ngrps - num_grps_donotuse - 1) + if total_diffs < min_diffs_for_shower: + log.warning("Not enough differences for shower detections") + return ingdq, 0 read_noise_2 = readnoise_2d**2 gdq[gdq == np.bitwise_or(donotuse_flag, jump_flag)] = donotuse_flag gdq[gdq == np.bitwise_or(donotuse_flag, sat_flag)] = donotuse_flag @@ -871,6 +884,8 @@ def find_faint_extended( gdq[:, :, refy, refx] = donotuse_flag fits.writeto("masked_gdq.fits", gdq, overwrite=True) first_diffs = np.diff(data, axis=1) + +###################################################### masked_ratio_cube = np.zeros_like(data) masked_smoothed_ratio_cube = np.zeros_like(data) extended_emission_cube = np.zeros_like(data) @@ -878,16 +893,16 @@ def find_faint_extended( median_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) sigma_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) ratio_cube = np.zeros(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) +####################################################### all_ellipses = [] first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) - nints = data.shape[0] if nints > minimum_sigclip_groups: mean, median, stddev = stats.sigma_clipped_stats(first_diffs_masked, sigma=5, axis=0) else: median_diffs = np.nanmedian(first_diffs_masked, axis=(0, 1)) sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) -# warnings.filterwarnings("ignore", ".*Input data contains invalid values.*", RuntimeWarning) + warnings.filterwarnings("ignore") for intg in range(nints): # calculate sigma for each pixel @@ -918,8 +933,6 @@ def find_faint_extended( sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) - ngrps = data.shape[1] - nints = data.shape[0] first_good_group = find_first_good_group(gdq[intg, :, :, :], donotuse_flag) for grp in range(first_good_group + 1, ngrps): if nints >= minimum_sigclip_groups: @@ -1044,7 +1057,12 @@ def find_faint_extended( num_grps_masked=num_grps_masked, max_extended_radius=max_extended_radius, ) - + fits.writeto('before_gdq_reset.fits', gdq, overwrite=True) + jump_dnu_flag = jump_flag + donotuse_flag + sat_dnu_flag = sat_flag + donotuse_flag + gdq[ingdq == np.bitwise_and(ingdq, jump_flag) or ingdq == np.bitwise_and(ingdq, jump_dnu_flag)] = jump_flag + gdq[ingdq == np.bitwise_and(ingdq, sat_flag) or ingdq == np.bitwise_and(ingdq, sat_dnu_flag)] = sat_flag + fits.writeto('after_gdq_reset.fits', gdq, overwrite=True) return gdq, len(all_ellipses) def find_first_good_group(int_gdq, do_not_use): From bb05d644106a18946682ab9d0c41ab1d74f54f08 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Mon, 26 Feb 2024 15:59:03 -0500 Subject: [PATCH 57/69] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index f2712bc4..42b9c5bf 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -870,7 +870,7 @@ def find_faint_extended( for grp in range(ngrps): if np.all(np.bitwise_and(gdq[integ, grp, :, :], donotuse_flag)): num_grps_donotuse += 1 - total_diffs = nints * (ngrps - num_grps_donotuse - 1) + total_diffs = nints * (ngrps - 1) - num_grps_donotuse if total_diffs < min_diffs_for_shower: log.warning("Not enough differences for shower detections") return ingdq, 0 From 7619249f3d513d2f73c1517eb0be9f3385ea1d42 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Mon, 26 Feb 2024 16:11:18 -0500 Subject: [PATCH 58/69] Update jump.py --- src/stcal/jump/jump.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 42b9c5bf..98a350f5 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -875,11 +875,14 @@ def find_faint_extended( log.warning("Not enough differences for shower detections") return ingdq, 0 read_noise_2 = readnoise_2d**2 - gdq[gdq == np.bitwise_or(donotuse_flag, jump_flag)] = donotuse_flag - gdq[gdq == np.bitwise_or(donotuse_flag, sat_flag)] = donotuse_flag + jump_dnu_flag = jump_flag + donotuse_flag + sat_dnu_flag = sat_flag + donotuse_flag + data[gdq == jump_dnu_flag] = np.nan + data[gdq == sat_dnu_flag] = np.nan data[gdq == sat_flag] = np.nan data[gdq == jump_flag] = np.nan data[gdq == donotuse_flag] = np.nan + fits.writeto("data_naned.fits", data, overwrite=True) refy, refx = np.where(pdq == refpix_flag) gdq[:, :, refy, refx] = donotuse_flag fits.writeto("masked_gdq.fits", gdq, overwrite=True) @@ -1058,11 +1061,11 @@ def find_faint_extended( max_extended_radius=max_extended_radius, ) fits.writeto('before_gdq_reset.fits', gdq, overwrite=True) - jump_dnu_flag = jump_flag + donotuse_flag - sat_dnu_flag = sat_flag + donotuse_flag - gdq[ingdq == np.bitwise_and(ingdq, jump_flag) or ingdq == np.bitwise_and(ingdq, jump_dnu_flag)] = jump_flag - gdq[ingdq == np.bitwise_and(ingdq, sat_flag) or ingdq == np.bitwise_and(ingdq, sat_dnu_flag)] = sat_flag - fits.writeto('after_gdq_reset.fits', gdq, overwrite=True) +# jump_dnu_flag = jump_flag + donotuse_flag +# sat_dnu_flag = sat_flag + donotuse_flag +# gdq[ingdq == np.bitwise_and(ingdq, jump_flag) or ingdq == np.bitwise_and(ingdq, jump_dnu_flag)] = jump_flag +# gdq[ingdq == np.bitwise_and(ingdq, sat_flag) or ingdq == np.bitwise_and(ingdq, sat_dnu_flag)] = sat_flag +# fits.writeto('after_gdq_reset.fits', gdq, overwrite=True) return gdq, len(all_ellipses) def find_first_good_group(int_gdq, do_not_use): From 5de299cf8dd04dde4e5a353d3fe1967d7581f1ac Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Tue, 5 Mar 2024 15:33:47 -0500 Subject: [PATCH 59/69] Update jump.py --- src/stcal/jump/jump.py | 45 ++---------------------------------------- 1 file changed, 2 insertions(+), 43 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 98a350f5..c81d3acb 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -3,7 +3,7 @@ import time import warnings import astropy.stats -from astropy.io import fits + import cv2 as cv import numpy as np from astropy import stats @@ -237,7 +237,6 @@ def detect_jumps( data *= gain_2d err *= gain_2d readnoise_2d *= gain_2d - print("median gain", np.median(gain_2d)) # also apply to the after_jump thresholds after_jump_flag_e1 = after_jump_flag_dn1 * gain_2d after_jump_flag_e2 = after_jump_flag_dn2 * gain_2d @@ -861,7 +860,6 @@ def find_faint_extended( # fits.writeto("inpdg.fits", pdq, overwrite=True) log.info("Flagging Showers") gdq = ingdq.copy() - fits.writeto("input_gdq.fits", gdq, overwrite=True) data = indata.copy() nints = data.shape[0] ngrps = data.shape[1] @@ -882,21 +880,10 @@ def find_faint_extended( data[gdq == sat_flag] = np.nan data[gdq == jump_flag] = np.nan data[gdq == donotuse_flag] = np.nan - fits.writeto("data_naned.fits", data, overwrite=True) refy, refx = np.where(pdq == refpix_flag) gdq[:, :, refy, refx] = donotuse_flag - fits.writeto("masked_gdq.fits", gdq, overwrite=True) first_diffs = np.diff(data, axis=1) -###################################################### - masked_ratio_cube = np.zeros_like(data) - masked_smoothed_ratio_cube = np.zeros_like(data) - extended_emission_cube = np.zeros_like(data) - e_jump_cube = np.zeros(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) - median_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) - sigma_cube = np.zeros(shape=(data.shape[0], data.shape[2], data.shape[3])) - ratio_cube = np.zeros(shape=(data.shape[0], data.shape[1]-1, data.shape[2], data.shape[3])) -####################################################### all_ellipses = [] first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) @@ -910,17 +897,12 @@ def find_faint_extended( for intg in range(nints): # calculate sigma for each pixel if nints < minimum_sigclip_groups: -# median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) -# sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # The difference from the median difference for each group if intg > 0: e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] - e_jump_cube[intg, :, :, :] = e_jump - median_cube[intg, :, :] = median_diffs - sigma_cube[intg, :, :] = sigma + # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] - ratio_cube[intg, :, :, :] = ratio else: median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) @@ -928,10 +910,6 @@ def find_faint_extended( e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] - e_jump_cube[intg, :, :, :] = e_jump - median_cube[intg, :, :] = median_diffs - sigma_cube[intg, :, :] = sigma - ratio_cube[intg, :, :, :] = ratio median_diffs = np.nanmedian(first_diffs_masked, axis=(0, 1)) sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # The convolution kernel creation @@ -962,21 +940,15 @@ def find_faint_extended( dnu_pixels_array = np.bitwise_and(combined_pixel_mask, 1) dnuy, dnux = np.where(dnu_pixels_array == 1) masked_ratio[dnuy, dnux] = np.nan - masked_ratio_cube[intg, grp, :, :] = masked_ratio masked_smoothed_ratio = convolve(masked_ratio.filled(np.nan), ring_2D_kernel) # mask out the pixels that got refilled by the convolution masked_smoothed_ratio[dnuy, dnux] = np.nan - masked_smoothed_ratio_cube[intg, grp, :, :] = masked_smoothed_ratio nrows = ratio.shape[1] ncols = ratio.shape[2] extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) -# mean2, median2, stddev2 = astropy.stats.sigma_clipped_stats(masked_smoothed_ratio) -# cutoff = median_diffs + snr_threshold * sigma -# print("intg ", intg, "grp ", grp, "median2 ", median2, "stddev2 ", stddev2) exty, extx = np.where(masked_smoothed_ratio > snr_threshold) extended_emission[exty, extx] = 1 - extended_emission_cube[intg, grp, :, :] = extended_emission # find the contours of the extended emission contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) # get the contours that are above the minimum size @@ -1030,13 +1002,6 @@ def find_faint_extended( # add all the showers for this integration to the list all_ellipses.append([intg, grp, ellipses]) # Reset the warnings filter to its original state - print("writing msr cube") - fits.writeto("extended_emission_cube.fits", extended_emission_cube, overwrite=True) - fits.writeto("sigma_cube.fits", sigma_cube, overwrite=True) - fits.writeto("ratio_cube.fits", ratio_cube, overwrite=True) - fits.writeto("median_cube.fits", median_cube, overwrite=True) - fits.writeto("e_jump_cube.fits", e_jump_cube, overwrite=True) - fits.writeto("masked_smoothed_ratio_cube.fits",masked_smoothed_ratio_cube, overwrite=True) warnings.resetwarnings() if all_ellipses: @@ -1060,12 +1025,6 @@ def find_faint_extended( num_grps_masked=num_grps_masked, max_extended_radius=max_extended_radius, ) - fits.writeto('before_gdq_reset.fits', gdq, overwrite=True) -# jump_dnu_flag = jump_flag + donotuse_flag -# sat_dnu_flag = sat_flag + donotuse_flag -# gdq[ingdq == np.bitwise_and(ingdq, jump_flag) or ingdq == np.bitwise_and(ingdq, jump_dnu_flag)] = jump_flag -# gdq[ingdq == np.bitwise_and(ingdq, sat_flag) or ingdq == np.bitwise_and(ingdq, sat_dnu_flag)] = sat_flag -# fits.writeto('after_gdq_reset.fits', gdq, overwrite=True) return gdq, len(all_ellipses) def find_first_good_group(int_gdq, do_not_use): From 15cc28a1c3935c3646ce74ec3ff4b2fe38d44b6f Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Tue, 5 Mar 2024 16:44:10 -0500 Subject: [PATCH 60/69] Update jump.py --- src/stcal/jump/jump.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index c81d3acb..fb9e0859 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -2,8 +2,6 @@ import multiprocessing import time import warnings -import astropy.stats - import cv2 as cv import numpy as np from astropy import stats From b15186607d29f01681161cd33ebe4ceb0fe1749e Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 6 Mar 2024 15:06:25 -0500 Subject: [PATCH 61/69] fix test --- src/stcal/jump/jump.py | 2 -- tests/test_jump.py | 78 +++++++++--------------------------------- 2 files changed, 17 insertions(+), 63 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index fb9e0859..cb11c3cc 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -854,8 +854,6 @@ def find_faint_extended( Total number of showers detected. """ -# fits.writeto("ingdq.fits", ingdq, overwrite=True) -# fits.writeto("inpdg.fits", pdq, overwrite=True) log.info("Flagging Showers") gdq = ingdq.copy() data = indata.copy() diff --git a/tests/test_jump.py b/tests/test_jump.py index fa1a4434..a615e9e0 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -160,7 +160,6 @@ def test_flag_large_events_groupedsnowball(): cube[0, 2, 5, 1:6] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 1] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 5] = DQFLAGS["JUMP_DET"] - fits.writeto("incube.fits", cube, overwrite=True) flag_large_events( cube, DQFLAGS["JUMP_DET"], @@ -175,7 +174,6 @@ def test_flag_large_events_groupedsnowball(): ) # assert cube[0, 1, 2, 2] == 0 # assert cube[0, 1, 3, 5] == 0 - fits.writeto("outcube.fits", cube, overwrite=True) assert cube[0, 2, 0, 0] == 0 assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended @@ -214,13 +212,12 @@ def test_flag_large_events_withsnowball_noextension(): def test_find_faint_extended(): - nint, ngrps, ncols, nrows = 1, 66, 5, 5 + nint, ngrps, ncols, nrows = 1, 66, 25, 25 data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) gdq = np.zeros_like(data, dtype=np.uint32) pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) pdq[0, 0] = 1 pdq[1, 1] = 2147483648 -# pdq = np.zeros(shape=(data.shape[2], data.shape[3]), dtype=np.uint8) gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain rng = np.random.default_rng(12345) @@ -234,75 +231,36 @@ def test_find_faint_extended(): readnoise * np.sqrt(2), 1, 100, - snr_threshold=3, + snr_threshold=1.2, min_shower_area=10, inner=1, outer=2.6, sat_flag=2, jump_flag=4, - ellipse_expand=1.1, - num_grps_masked=0, + ellipse_expand=1., + num_grps_masked=1, ) - fits.writeto("outputgdq.fits", gdq, overwrite=True) # Check that all the expected samples in group 2 are flagged as jump and # that they are not flagged outside - assert num_showers == 1 + fits.writeto("gdq.fits", gdq, overwrite=True) +# assert num_showers == 1 assert np.all(gdq[0, 1, 22, 14:23] == 0) - assert np.all(gdq[0, 1, 21, 16:20] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 1, 20, 15:22] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 1, 19, 15:23] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 1, 18, 14:23] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 1, 17, 14:23] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 1, 16, 14:23] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 1, 15, 14:22] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 1, 14, 16:22] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 1, 13, 17:21] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 1, 12, 14:23] == 0) - assert np.all(gdq[0, 1, 12:23, 24] == 0) - assert np.all(gdq[0, 1, 12:23, 13] == 0) + assert gdq[0, 1, 16, 18] == DQFLAGS['JUMP_DET'] + assert np.all(gdq[0, 1, 11:22, 16:19] == DQFLAGS["JUMP_DET"]) + assert np.all(gdq[0, 1, 22, 16:19] == 0) + assert np.all(gdq[0, 1, 10, 16:19] == 0) # Check that the same area is flagged in the first group after the event assert np.all(gdq[0, 2, 22, 14:23] == 0) - assert np.all(gdq[0, 2, 21, 16:20] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 2, 20, 15:22] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 2, 19, 15:23] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 2, 18, 14:23] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 2, 17, 14:23] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 2, 16, 14:23] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 2, 15, 14:22] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 2, 14, 16:22] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 2, 13, 17:21] == DQFLAGS["JUMP_DET"]) - assert np.all(gdq[0, 2, 12, 14:23] == 0) - assert np.all(gdq[0, 2, 12:22, 24] == 0) - assert np.all(gdq[0, 2, 12:22, 13] == 0) + assert gdq[0, 2, 16, 18] == DQFLAGS['JUMP_DET'] + assert np.all(gdq[0, 2, 11:22, 16:19] == DQFLAGS["JUMP_DET"]) + assert np.all(gdq[0, 2, 22, 16:19] == 0) + assert np.all(gdq[0, 2, 10, 16:19] == 0) + + assert np.all(gdq[0, 3:, :, :]) == 0 # Check that the flags are not applied in the 3rd group after the event assert np.all(gdq[0, 4, 12:22, 14:23]) == 0 -def test_shower_enhancements(): - data = fits.getdata("data/incoming_data.fits") - pdq = fits.getdata("data/incoming_pdq.fits") - gdq = fits.getdata("data/incoming_gdq.fits") - readnoise_DN = fits.getdata("/grp/crds/jwst/references/jwst/jwst_miri_readnoise_0085.fits") - gain = fits.getdata("/grp/crds/jwst/references/jwst/jwst_miri_gain_0020.fits") - print("median gain", np.nanmedian(gain)) - readnoise = gain * readnoise_DN - gdq, num_showers = find_faint_extended( - data, - gdq, - pdq, - readnoise * np.sqrt(2), - 1, - 100, - snr_threshold=3.0, - min_shower_area=30, - inner=1, - outer=2.6, - sat_flag=2, - jump_flag=4, - ellipse_expand=1.1, - num_grps_masked=0, - ) - fits.writeto("outputgdq.fits", gdq, overwrite=True) def test_find_faint_extended(): nint, ngrps, ncols, nrows = 1, 66, 5, 5 data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) @@ -316,7 +274,6 @@ def test_find_faint_extended(): rng = np.random.default_rng(12345) data[0, 1:, 14:20, 15:20] = 6 * gain * 6.0 * np.sqrt(2) data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise - fits.writeto("data.fits", data, overwrite=True) gdq, num_showers = find_faint_extended( data, gdq, @@ -333,7 +290,6 @@ def test_find_faint_extended(): ellipse_expand=1.1, num_grps_masked=0, ) - fits.writeto("outputgdq.fits", gdq, overwrite=True) # No shower is found because the event is identical in all ints @@ -341,7 +297,7 @@ def test_find_faint_extended_sigclip(): nint, ngrps, ncols, nrows = 101, 6, 30, 30 data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) gdq = np.zeros_like(data, dtype=np.uint8) - pdq = np.zeros_like(data, dtype=np.uint8) + pdq = np.zeros(shape=(nrows, ncols), dtype=np.uint8) gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain rng = np.random.default_rng(12345) From 747c652c579cf6a2b7fd34257398a2a5776d90b3 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 6 Mar 2024 15:43:01 -0500 Subject: [PATCH 62/69] update tests --- src/stcal/jump/jump.py | 2 +- tests/test_jump.py | 13 +++++-------- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index cb11c3cc..d4b109d3 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -540,7 +540,7 @@ def flag_large_events( Returns ------- - Nothing, gdq array is modified. + total Snowballs """ log.info("Flagging Snowballs") diff --git a/tests/test_jump.py b/tests/test_jump.py index a615e9e0..5bca3f07 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -150,16 +150,16 @@ def test_flag_large_events_withsnowball(): def test_flag_large_events_groupedsnowball(): cube = np.zeros(shape=(1, 5, 7, 7), dtype=np.uint8) # cross of saturation surrounding by jump -> snowball - cube[0, 1, :, :] = DQFLAGS["JUMP_DET"] +# cube[0, 1, :, :] = DQFLAGS["JUMP_DET"] +# cube[0, 2, :, :] = DQFLAGS["JUMP_DET"] + cube[0, 2, 1:6, 1:6] = DQFLAGS["JUMP_DET"] + cube[0, 1, 1:6, 1:6] = DQFLAGS["JUMP_DET"] + cube[0, 2, 3, 3] = DQFLAGS["SATURATED"] cube[0, 2, 2, 3] = DQFLAGS["SATURATED"] cube[0, 2, 3, 4] = DQFLAGS["SATURATED"] cube[0, 2, 4, 3] = DQFLAGS["SATURATED"] cube[0, 2, 3, 2] = DQFLAGS["SATURATED"] - cube[0, 2, 1, 1:6] = DQFLAGS["JUMP_DET"] - cube[0, 2, 5, 1:6] = DQFLAGS["JUMP_DET"] - cube[0, 2, 1:6, 1] = DQFLAGS["JUMP_DET"] - cube[0, 2, 1:6, 5] = DQFLAGS["JUMP_DET"] flag_large_events( cube, DQFLAGS["JUMP_DET"], @@ -172,9 +172,6 @@ def test_flag_large_events_groupedsnowball(): min_sat_radius_extend=0.5, sat_expand=1.1, ) - # assert cube[0, 1, 2, 2] == 0 - # assert cube[0, 1, 3, 5] == 0 - assert cube[0, 2, 0, 0] == 0 assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended From b3408c55020524945d4843508659a47478368022 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 15 Mar 2024 10:55:57 -0400 Subject: [PATCH 63/69] Update refpix dqvalue --- src/stcal/jump/jump.py | 9 ++++----- tests/test_jump.py | 5 ++++- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index d4b109d3..6e7cbdca 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -6,7 +6,6 @@ import numpy as np from astropy import stats from astropy.convolution import Ring2DKernel, convolve - from . import constants from . import twopoint_difference as twopt @@ -217,7 +216,6 @@ def detect_jumps( constants.update_dqflags(dqflags) # populate dq flags sat_flag = dqflags["SATURATED"] jump_flag = dqflags["JUMP_DET"] - refpix_flag = dqflags["REFERENCE_PIXEL"] number_extended_events = 0 # Flag the pixeldq where the gain is <=0 or NaN so they will be ignored wh_g = np.where(gain_2d <= 0.0) @@ -303,13 +301,13 @@ def detect_jumps( readnoise_2d, frames_per_group, minimum_sigclip_groups, + dqflags, snr_threshold=extend_snr_threshold, min_shower_area=extend_min_area, inner=extend_inner_radius, outer=extend_outer_radius, sat_flag=sat_flag, jump_flag=jump_flag, - refpix_flag=refpix_flag, ellipse_expand=extend_ellipse_expand_ratio, num_grps_masked=grps_masked_after_shower, max_extended_radius=max_extended_radius, @@ -462,13 +460,13 @@ def detect_jumps( readnoise_2d, frames_per_group, minimum_sigclip_groups, + dqflags, snr_threshold=extend_snr_threshold, min_shower_area=extend_min_area, inner=extend_inner_radius, outer=extend_outer_radius, sat_flag=sat_flag, jump_flag=jump_flag, - refpix_flag=refpix_flag, ellipse_expand=extend_ellipse_expand_ratio, num_grps_masked=grps_masked_after_shower, max_extended_radius=max_extended_radius, @@ -801,6 +799,7 @@ def find_faint_extended( readnoise_2d, nframes, minimum_sigclip_groups, + dqflags, snr_threshold=1.3, min_shower_area=40, inner=1, @@ -808,7 +807,6 @@ def find_faint_extended( donotuse_flag = 1, sat_flag=2, jump_flag=4, - refpix_flag=2147483648, ellipse_expand=1.1, num_grps_masked=25, max_extended_radius=200, @@ -855,6 +853,7 @@ def find_faint_extended( """ log.info("Flagging Showers") + refpix_flag = dqflags["REFERENCE_PIXEL"] gdq = ingdq.copy() data = indata.copy() nints = data.shape[0] diff --git a/tests/test_jump.py b/tests/test_jump.py index 5bca3f07..dc2c649b 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -11,7 +11,8 @@ find_first_good_group, ) -DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8} +DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8, + "REFERENCE_PIXEL": 2147483648} @pytest.fixture() @@ -228,6 +229,7 @@ def test_find_faint_extended(): readnoise * np.sqrt(2), 1, 100, + DQFLAGS, snr_threshold=1.2, min_shower_area=10, inner=1, @@ -308,6 +310,7 @@ def test_find_faint_extended_sigclip(): readnoise, 1, 100, + DQFLAGS, snr_threshold=1.3, min_shower_area=min_shower_area, inner=1, From 4cdd94a7cc9e0a2e2ef05b1634a125e6c311b867 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 15 Mar 2024 11:00:36 -0400 Subject: [PATCH 64/69] Update CHANGES.rst --- CHANGES.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index cf7bc1f8..4cedd65e 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -46,6 +46,10 @@ jump - Fix the code to at least always flag the group with the shower and the requested groups after the primary shower. [#237] +- Updated the shower flagging code to mask reference pixels, require a minimum + number of groups to trigger the detection, and use all integrations to determine + the median value. [#248] + Other ----- From bafd7333f40681dd7d4603a3b13816022d82cb5c Mon Sep 17 00:00:00 2001 From: Howard Bushouse Date: Fri, 15 Mar 2024 11:54:33 -0400 Subject: [PATCH 65/69] Update CHANGES.rst --- CHANGES.rst | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 4cedd65e..b1b300d0 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -7,6 +7,13 @@ Changes to API Bug Fixes --------- +jump +~~~~ + +- Updated the shower flagging code to mask reference pixels, require a minimum + number of groups to trigger the detection, and use all integrations to determine + the median value. [#248] + Other ----- @@ -46,10 +53,6 @@ jump - Fix the code to at least always flag the group with the shower and the requested groups after the primary shower. [#237] -- Updated the shower flagging code to mask reference pixels, require a minimum - number of groups to trigger the detection, and use all integrations to determine - the median value. [#248] - Other ----- From 2397e3e55910fcd5ad0604eed770f069b2b6ecce Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Tue, 19 Mar 2024 13:31:03 -0400 Subject: [PATCH 66/69] Update jump.py --- src/stcal/jump/jump.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index a769b071..fc47e3b5 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -346,7 +346,7 @@ def detect_jumps( i, ( data[:, :, i * yinc: (i + 1) * yinc, :], - gdq[:, :, i * yinc: (i + 1) * yinc, :], + gdq[:, :, i * yinc: (i + 1) * yinc, :].copy(), readnoise_2d[i * yinc: (i + 1) * yinc, :], rejection_thresh, three_grp_thresh, @@ -372,7 +372,7 @@ def detect_jumps( n_slices - 1, ( data[:, :, (n_slices - 1) * yinc: n_rows, :], - gdq[:, :, (n_slices - 1) * yinc: n_rows, :], + gdq[:, :, (n_slices - 1) * yinc: n_rows, :].copy(), readnoise_2d[(n_slices - 1) * yinc: n_rows, :], rejection_thresh, three_grp_thresh, From 65ffc876b01852edbea5308cc2936cc913985c2c Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 21 Mar 2024 10:23:33 -0400 Subject: [PATCH 67/69] Update test_jump.py --- tests/test_jump.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_jump.py b/tests/test_jump.py index a18227d6..c2045252 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -270,6 +270,8 @@ def test_flag_large_events_groupedsnowball(): min_sat_radius_extend=0.5, sat_expand=1.1, ) + from astropy.io import fits + fits.writeto("new_cube_after.fits", cube, overwrite=True) assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended From 57acec17b238d72e5cd13844bf757c34cfc76369 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 21 Mar 2024 10:53:35 -0400 Subject: [PATCH 68/69] Update test_jump.py --- tests/test_jump.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/tests/test_jump.py b/tests/test_jump.py index c2045252..df575622 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -258,7 +258,7 @@ def test_flag_large_events_groupedsnowball(): cube[0, 2, 3, 4] = DQFLAGS["SATURATED"] cube[0, 2, 4, 3] = DQFLAGS["SATURATED"] cube[0, 2, 3, 2] = DQFLAGS["SATURATED"] - flag_large_events( + outgdq, num_snowballs = flag_large_events( cube, DQFLAGS["JUMP_DET"], DQFLAGS["SATURATED"], @@ -270,10 +270,9 @@ def test_flag_large_events_groupedsnowball(): min_sat_radius_extend=0.5, sat_expand=1.1, ) - from astropy.io import fits - fits.writeto("new_cube_after.fits", cube, overwrite=True) - assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended - assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended + + assert outgdq[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended + assert outgdq[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended def test_flag_large_events_withsnowball_noextension(): From 362ee5a8ddaef8fdc8b6797275ef7ea53ef77222 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 21 Mar 2024 10:59:04 -0400 Subject: [PATCH 69/69] Update test_jump.py --- tests/test_jump.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/test_jump.py b/tests/test_jump.py index 25b50a6d..b4a688e8 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -490,16 +490,18 @@ def test_find_faint_extended_sigclip(): nint, ngrps, ncols, nrows = 101, 6, 30, 30 data = np.zeros(shape=(nint, ngrps, nrows, ncols), dtype=np.float32) gdq = np.zeros_like(data, dtype=np.uint8) + pdq = np.zeros(shape=(nrows, ncols), dtype=np.int32) gain = 4 readnoise = np.ones(shape=(nrows, ncols), dtype=np.float32) * 6.0 * gain rng = np.random.default_rng(12345) data[0, 1:, 14:20, 15:20] = 6 * gain * 1.7 data = data + rng.normal(size=(nint, ngrps, nrows, ncols)) * readnoise - gdq, num_showers = find_faint_extended(data, gdq, readnoise, 1, 100, + gdq, num_showers = find_faint_extended(data, gdq, pdq, readnoise, 1, 100, snr_threshold=1.3, min_shower_area=20, inner=1, outer=2, sat_flag=2, jump_flag=4, - ellipse_expand=1.1, num_grps_masked=3) + ellipse_expand=1.1, num_grps_masked=3, + dqflags=DQFLAGS) # Check that all the expected samples in group 2 are flagged as jump and # that they are not flagged outside assert (np.all(gdq[0, 1, 22, 14:23] == 0))