From 22878e6875937a498c4e8b264f299d6bcdd64dc3 Mon Sep 17 00:00:00 2001 From: David Law Date: Mon, 30 Sep 2024 17:06:41 -0700 Subject: [PATCH 1/3] Add flashfrac --- src/stcal/jump/jump.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 9eaecf664..5a0fe9d60 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -60,6 +60,7 @@ def detect_jumps( min_diffs_single_pass=10, mask_persist_grps_next_int=True, persist_grps_flagged=25, + mmflashfrac=1.0, ): """ This is the high-level controlling routine for the jump detection process. @@ -221,6 +222,11 @@ def detect_jumps( min_diffs_single_pass : int The minimum number of groups to switch to flagging all outliers in a single pass. + mmflashfrac: float + Fraction of the array in a given group that has to have a jump detected + in order to flag the entire array as a jump. This is designed to deal with + sporadic micrometeorite flashes. + Returns ------- gdq : int, 4D array @@ -338,6 +344,12 @@ def detect_jumps( ) log.info("Total showers= %i", num_showers) number_extended_events = num_showers + + # This is where we look for micrometeorite flashes if the triggering + # threshold is less than the entire array + if (mmflashfrac < 1): + gdq = find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac) + else: yinc = int(n_rows // n_slices) slices = [] @@ -503,6 +515,12 @@ def detect_jumps( ) log.info("Total showers= %i", num_showers) number_extended_events = num_showers + + # This is where we look for micrometeorite flashes if the triggering + # threshold is less than the entire array + if (mmflashfrac < 1): + gdq = find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac) + elapsed = time.time() - start log.info("Total elapsed time = %g sec", elapsed) @@ -514,6 +532,21 @@ def detect_jumps( # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev +# Find micrometeorite flashes that light up the entire array in a small number +# of groups. Do this by looking for cases where greater than mmflash_pct percent +# of the array is flagged as 'JUMP_DET' and flag the entire array. +def find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac): + npix = gdq.shape[2]*gdq.shape[3] # Number of pixels in array + # Loop over integrations + for ii in range(0,gdq.shape[0]): + # Loop over groups + for jj in range(0,gdq.shape[1]): + indx = np.where(gdq[ii,jj,:,:] & jump_flag != 0) + fraction = float(len(indx[0])) / npix + if (fraction > mmflashfrac): + gdq[ii,jj,:,:] = gdq[ii,jj,:,:] | jump_flag + log.info("Detected a micrometeorite flash in integ = %i, group= %i", ii, jj) + return gdq def flag_large_events( gdq, From 66f0bb691d7959629537536356b179741afaffaf Mon Sep 17 00:00:00 2001 From: David Law Date: Thu, 17 Oct 2024 15:48:41 -0400 Subject: [PATCH 2/3] Code style cleanup --- src/stcal/jump/jump.py | 175 ++++++++++++++++++----------------------- 1 file changed, 77 insertions(+), 98 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 5a0fe9d60..3a6e2b4ec 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -304,51 +304,6 @@ def detect_jumps( 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: - gdq, total_snowballs = flag_large_events( - gdq, - jump_flag, - sat_flag, - min_sat_area=min_sat_area, - min_jump_area=min_jump_area, - expand_factor=expand_factor, - sat_required_snowball=sat_required_snowball, - min_sat_radius_extend=min_sat_radius_extend, - edge_size=edge_size, - sat_expand=sat_expand, - max_extended_radius=max_extended_radius, - mask_persist_grps_next_int=mask_persist_grps_next_int, - persist_grps_flagged=persist_grps_flagged, - ) - log.info("Total snowballs = %i", total_snowballs) - number_extended_events = total_snowballs - if find_showers: - gdq, num_showers = find_faint_extended( - data, - gdq, - pdq, - 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, - ellipse_expand=extend_ellipse_expand_ratio, - num_grps_masked=grps_masked_after_shower, - max_extended_radius=max_extended_radius, - ) - log.info("Total showers= %i", num_showers) - number_extended_events = num_showers - - # This is where we look for micrometeorite flashes if the triggering - # threshold is less than the entire array - if (mmflashfrac < 1): - gdq = find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac) else: yinc = int(n_rows // n_slices) @@ -475,51 +430,53 @@ def detect_jumps( 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: - gdq, total_snowballs = flag_large_events( - gdq, - jump_flag, - sat_flag, - min_sat_area=min_sat_area, - min_jump_area=min_jump_area, - expand_factor=expand_factor, - sat_required_snowball=sat_required_snowball, - min_sat_radius_extend=min_sat_radius_extend, - edge_size=edge_size, - sat_expand=sat_expand, - max_extended_radius=max_extended_radius, - mask_persist_grps_next_int=mask_persist_grps_next_int, - persist_grps_flagged=persist_grps_flagged, - ) - log.info("Total snowballs = %i", total_snowballs) - number_extended_events = total_snowballs - if find_showers: - gdq, num_showers = find_faint_extended( - data, - gdq, - pdq, - 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, - ellipse_expand=extend_ellipse_expand_ratio, - num_grps_masked=grps_masked_after_shower, - max_extended_radius=max_extended_radius, - ) - log.info("Total showers= %i", num_showers) - number_extended_events = num_showers + # Look for snowballs in near-IR data + if expand_large_events: + gdq, total_snowballs = flag_large_events( + gdq, + jump_flag, + sat_flag, + min_sat_area=min_sat_area, + min_jump_area=min_jump_area, + expand_factor=expand_factor, + sat_required_snowball=sat_required_snowball, + min_sat_radius_extend=min_sat_radius_extend, + edge_size=edge_size, + sat_expand=sat_expand, + max_extended_radius=max_extended_radius, + mask_persist_grps_next_int=mask_persist_grps_next_int, + persist_grps_flagged=persist_grps_flagged, + ) + log.info("Total snowballs = %i", total_snowballs) + number_extended_events = total_snowballs + + # Look for showers in mid-IR data + if find_showers: + gdq, num_showers = find_faint_extended( + data, + gdq, + pdq, + 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, + ellipse_expand=extend_ellipse_expand_ratio, + num_grps_masked=grps_masked_after_shower, + max_extended_radius=max_extended_radius, + ) + log.info("Total showers= %i", num_showers) + number_extended_events = num_showers - # This is where we look for micrometeorite flashes if the triggering - # threshold is less than the entire array - if (mmflashfrac < 1): - gdq = find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac) + # Look for micrometeorite flashes if the triggering + # threshold is less than the entire array + if (mmflashfrac < 1): + gdq = find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac) elapsed = time.time() - start log.info("Total elapsed time = %g sec", elapsed) @@ -532,20 +489,42 @@ def detect_jumps( # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev -# Find micrometeorite flashes that light up the entire array in a small number -# of groups. Do this by looking for cases where greater than mmflash_pct percent -# of the array is flagged as 'JUMP_DET' and flag the entire array. +# Find micrometeorite flashes def find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac): - npix = gdq.shape[2]*gdq.shape[3] # Number of pixels in array + """ + This routine looks for micrometeorite flashes that light up the entire array in a small number + of groups. It does this by looking for cases where greater than mmflash_pct percent + of the array is flagged as 'JUMP_DET', and flags the entire array as JUMP_DET in such cases. + + Parameters + ---------- + gdq : int, 4D array + Group dq array + jump_flag : int + DQ flag for jump detection. + mmflashfrac : float + Fraction of the array that can be flagged as a jump before the entire array + will be flagged. + Returns + ------- + gdq : int, 4D array + Updated group dq array + + """ + log.info("Looking for micrometeorite flashes") + + nints, ngroups, nrows, ncols = gdq.shape + + npix = nrows * ncols # Number of pixels in array # Loop over integrations - for ii in range(0,gdq.shape[0]): + for ii in range(0, nints): # Loop over groups - for jj in range(0,gdq.shape[1]): - indx = np.where(gdq[ii,jj,:,:] & jump_flag != 0) + for jj in range(0, ngroups): + indx = np.where(gdq[ii, jj, :, :] & jump_flag != 0) fraction = float(len(indx[0])) / npix if (fraction > mmflashfrac): - gdq[ii,jj,:,:] = gdq[ii,jj,:,:] | jump_flag - log.info("Detected a micrometeorite flash in integ = %i, group= %i", ii, jj) + gdq[ii, jj, :, :] = gdq[ii, jj, :, :] | jump_flag + log.info(f"Detected a micrometeorite flash in integ = {ii}, group = {jj}") return gdq def flag_large_events( From 272dcfc0ceaf89ac0d2e6beb874950c2a0236a4a Mon Sep 17 00:00:00 2001 From: David Law Date: Mon, 21 Oct 2024 13:52:02 -0400 Subject: [PATCH 3/3] Safety catch to ensure we dont flag too much with micrometeorite routine --- src/stcal/jump/jump.py | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 3a6e2b4ec..507359f65 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -490,12 +490,16 @@ def detect_jumps( return gdq, pdq, total_primary_crs, number_extended_events, stddev # Find micrometeorite flashes -def find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac): +def find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac, mmflashnmax=5): """ This routine looks for micrometeorite flashes that light up the entire array in a small number of groups. It does this by looking for cases where greater than mmflash_pct percent of the array is flagged as 'JUMP_DET', and flags the entire array as JUMP_DET in such cases. + Since such flashes are rare (a few per instrument per year, although they can affect 2-3 groups) + this routine also applies a safety catch to ensure that not too many are flagged in any one + exposure due to unexpected detector effects. + Parameters ---------- gdq : int, 4D array @@ -505,6 +509,9 @@ def find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac): mmflashfrac : float Fraction of the array that can be flagged as a jump before the entire array will be flagged. + mmflashnmax : float + Maximum number of groups in any given exposure that can be affected by + micrometeorites before we declare a detection failure and flag nothing. Returns ------- gdq : int, 4D array @@ -516,6 +523,9 @@ def find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac): nints, ngroups, nrows, ncols = gdq.shape npix = nrows * ncols # Number of pixels in array + + # Initial loop over integrations + groups to find flashes + flash_int, flash_group = [], [] # Loop over integrations for ii in range(0, nints): # Loop over groups @@ -523,8 +533,20 @@ def find_micrometeorite_flashes(gdq, jump_flag, mmflashfrac): indx = np.where(gdq[ii, jj, :, :] & jump_flag != 0) fraction = float(len(indx[0])) / npix if (fraction > mmflashfrac): - gdq[ii, jj, :, :] = gdq[ii, jj, :, :] | jump_flag - log.info(f"Detected a micrometeorite flash in integ = {ii}, group = {jj}") + flash_int.append(ii) + flash_group.append(jj) + + # If an unrealistically large number of flashes were found, fail out + nflash=len(flash_int) + if (nflash > mmflashnmax): + log.info(f"Unreasonably large number of possible micrometeorite flashes detected ({nflash})") + log.info("Quitting micrometeorite routine and flagging nothing.") + else: + for kk in range(0, nflash): + ii, jj = flash_int[kk], flash_group[kk] + gdq[ii, jj, :, :] = gdq[ii, jj, :, :] | jump_flag + log.info(f"Flagged a micrometeorite flash in integ = {ii}, group = {jj}") + return gdq def flag_large_events(