JP-3737: Handle micrometeorite flashes

src/stcal/jump/jump.py

@@ -60,6 +60,7 @@
     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 @@
     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
@@ -298,46 +304,7 @@
             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
+
     else:
         yinc = int(n_rows // n_slices)
         slices = []
@@ -463,46 +430,54 @@
         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
+
+    # 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 +489,65 @@
     # Return the updated data quality arrays
     return gdq, pdq, total_primary_crs, number_extended_events, stddev
 
+# Find micrometeorite flashes
+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
+        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.
+    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
+        Updated group dq array
+
+    """
+    log.info("Looking for micrometeorite flashes")
+
+    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
+        for jj in range(0, ngroups):
+            indx = np.where(gdq[ii, jj, :, :] & jump_flag != 0)
+            fraction = float(len(indx[0])) / npix
+            if (fraction > mmflashfrac):
+                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(
     gdq,