ENH: Automatically reject bad frames due to cosmic rays

amical/data_processing.py

@@ -117,13 +117,22 @@ def select_data(cube, clip_fact=0.5, clip=False, verbose=True, display=True):
         ind_clip = []
         cube_cleaned_checked = np.array(good_fram)
 
+    cube_cleaned_checked = np.array(cube_cleaned_checked)
+
     ind_clip2 = np.where(fluxes <= limit_flux)[0]
     if ((worst_fr in ind_clip2) and clip) or (worst_fr in flag_fram):
         ext = "(rejected)"
     else:
         ext = ""
 
-    diffmm = 100 * abs(np.max(fluxes) - np.min(fluxes)) / med_flux
+    fluxes_check = np.array([x.sum() for x in cube_cleaned_checked])
+
+    best_fr = np.argmax(fluxes_check)
+    worst_fr = np.argmin(fluxes_check)
+
+    med_flux = np.median(fluxes_check)
+    std_flux = np.std(fluxes_check)
+    diffmm = 100 * abs(np.max(fluxes_check) - np.min(fluxes_check)) / med_flux
     if display:
         import matplotlib.pyplot as plt
         from matplotlib.colors import PowerNorm
@@ -429,7 +438,7 @@ def show_clean_params(
 
     noBadPixel = False
     bad_pix_x, bad_pix_y = [], []
-    if np.any(bmap0):
+    if np.any(bmap0) or np.any(ab0):
         if len(ab0) != 0:
             for j in range(len(ab0)):
                 bmap0[ab0[j][1], ab0[j][0]] = 1
@@ -544,6 +553,39 @@ def _remove_dark(img1, darkfile=None, ihdu=0, verbose=False):
     return img1
 
 
+def _cosmic_bad_frames_finder(data, f_kernel):
+    """
+    Summary
+    -------------
+    Review the datacube and attempt centering. If centering fails, mark the frame as
+    a "bad frame" and exclude it from further cleaning steps.
+
+    Parameters
+    ----------
+    `data` : {numpy.array}
+       datacube containing the NRM data.\n
+
+    Returns:
+    --------
+    `bad_frames` : {list}
+       List of bad frames index.\n
+    """
+    n_im = data.shape[0]
+    bad_frames = []
+    for i in range(n_im):
+        filtmed = f_kernel is not None
+        try:
+            crop_max(data[i], 64, iframe=i, filtmed=filtmed, f=f_kernel)[0]
+        except ValueError:
+            bad_frames.append(i)
+    print(
+        "[AMICAL] %i unusable frames have been identified in the data cube,"
+        % (len(bad_frames))
+        + " primarily due to cosmic rays or persistent bad pixels."
+    )
+    return bad_frames
+
+
 def clean_data(
     data,
     isz=None,
@@ -585,6 +627,8 @@ def clean_data(
 
     bad_map, add_bad = _get_3d_bad_pixels(bad_map, add_bad, data)
 
+    bad_frames = _cosmic_bad_frames_finder(data, f_kernel)
+
     for i in track(range(n_im), description="Cleaning"):
         img0 = data[i]
         img0 = _apply_edge_correction(img0, edge=edge)
@@ -617,34 +661,42 @@ def clean_data(
             img_biased = img1.copy()
         img_biased[img_biased < 0] = 0  # Remove negative pixels
 
-        if isz is not None:
-            # Get expected center for sky correction
-            filtmed = f_kernel is not None
-            im_rec_max = crop_max(
-                img_biased, isz, offx=offx, offy=offy, filtmed=filtmed, f=f_kernel
-            )[0]
-        else:
-            im_rec_max = img_biased.copy()
-
-        if (
-            (im_rec_max.shape[0] != im_rec_max.shape[1])
-            or (isz is not None and im_rec_max.shape[0] != isz)
-            or (isz is None and im_rec_max.shape[0] != img0.shape[0])
-        ):
-            l_bad_frame.append(i)
-        else:
-            if apod and window is not None:
-                img = apply_windowing(im_rec_max, window=window)
-            elif apod:
-                warnings.warn(
-                    "apod is set to True, but window is None. Skipping apodisation",
-                    RuntimeWarning,
-                    stacklevel=2,
-                )
-                img = im_rec_max.copy()
+        if i not in bad_frames:
+            if isz is not None:
+                # Get expected center for sky correction
+                filtmed = f_kernel is not None
+                im_rec_max = crop_max(
+                    img_biased,
+                    isz,
+                    iframe=i,
+                    offx=offx,
+                    offy=offy,
+                    filtmed=filtmed,
+                    f=f_kernel,
+                )[0]
+            else:
+                im_rec_max = img_biased.copy()
+
+            if (
+                (im_rec_max.shape[0] != im_rec_max.shape[1])
+                or (isz is not None and im_rec_max.shape[0] != isz)
+                or (isz is None and im_rec_max.shape[0] != img0.shape[0])
+            ):
+                l_bad_frame.append(i)
             else:
-                img = im_rec_max.copy()
+                if apod and window is not None:
+                    img = apply_windowing(im_rec_max, window=window)
+                elif apod:
+                    warnings.warn(
+                        "apod is set to True, but window is None. Skipping apodisation",
+                        RuntimeWarning,
+                        stacklevel=2,
+                    )
+                    img = im_rec_max.copy()
+                else:
+                    img = im_rec_max.copy()
             cube_cleaned.append(img)
+
     if verbose:
         print("Bad centering frame number:", l_bad_frame)
     cube_cleaned = np.array(cube_cleaned)

amical/tools.py

@@ -79,7 +79,7 @@ def find_max(img, filtmed=True, f=3):
     return X, Y
 
 
-def crop_max(img, dim, offx=0, offy=0, filtmed=True, f=3):
+def crop_max(img, dim, iframe=0, offx=0, offy=0, filtmed=True, f=3):
     """
     Summary
     -------------
@@ -112,7 +112,7 @@ def crop_max(img, dim, offx=0, offy=0, filtmed=True, f=3):
     if isz_max < dim:
         size_msg = (
             f"The specified cropped image size, {dim}, is greater than the distance to"
-            " the PSF center in at least one dimension. The max size for this image is"
+            f" the PSF center in at least one dimension (frame {iframe}). The max size for this image is"
             f" {isz_max}"
         )
         raise ValueError(size_msg)