deep copy of numpy array cutouts / masking applied in method

AegeanTools/source_finder.py

@@ -124,7 +124,7 @@ def find_islands(im, bkg, rms, seed_clip=5.0, flood_clip=4.0, log=log):
         # obey seed clip constraint
         if np.any(snr[xmin:xmax, ymin:ymax] > seed_clip):
 
-            data_box = copy.copy(
+            data_box = copy.deepcopy(
                 im[xmin:xmax, ymin:ymax]
             )  # copy so that we don't blank the master data
 
@@ -197,9 +197,10 @@ def estimate_parinfo_image(islands, im, rms, wcshelper,
         # set flags to be empty
         is_flag = 0x0
         [rmin, rmax], [cmin, cmax] = island.bounding_box
-        i_data = im[rmin:rmax, cmin:cmax]
-        i_rms = rms[rmin:rmax, cmin:cmax]
+        i_data = copy.deepcopy(im[rmin:rmax, cmin:cmax])
+        i_rms = copy.deepcopy(rms[rmin:rmax, cmin:cmax])
 
+        # Mask out the bad pixels
         island_mask = island.mask 
         i_data[island_mask] = np.nan
         i_rms[island_mask] = np.nan
@@ -737,20 +738,23 @@ def _gen_flood_wrap(self, data, rmsimg, innerclip,
             ):  # obey inner clip constraint
                 # self.log.info("{1} Island {0} is above the inner clip limit"
                 #               .format(i, data.shape))
-                data_box = copy.copy(
+
+                # Flag pixel that are either below the flood level, or belong to other
+                # islands that happen to be withing the bounding box
+                island_mask = (snr[xmin:xmax, ymin:ymax] < outerclip) | \
+                              (l[xmin:xmax, ymin:ymax] != i + 1)
+                data_box = copy.deepcopy(
                     data[xmin:xmax, ymin:ymax]
                 )  # copy so that we don't blank the master data
-                data_box[
-                    np.where(snr[xmin:xmax, ymin:ymax] < outerclip)
-                ] = np.nan  # blank pixels that are outside the outerclip
-                data_box[
-                    np.where(l[xmin:xmax, ymin:ymax] != i + 1)
-                ] = np.nan  # blank out other summits
+
+                # blank out the bad pixels
+                data_box[island_mask] = np.nan
                 # check if there are any pixels left unmasked
                 if not np.any(np.isfinite(data_box)):
                     # self.log.info("{1} Island {0} has no non-masked pixels"
                     #               .format(i,data.shape))
                     continue
+
                 if domask and (self.global_data.region is not None):
                     y, x = np.where(snr[xmin:xmax, ymin:ymax] >= outerclip)
                     # convert indices of this sub region to indices in
@@ -827,7 +831,7 @@ def estimate_lmfit_parinfo(
 
         # check to see if this island is a negative peak since we need to
         # treat such cases slightly differently
-        isnegative = max(data[np.where(np.isfinite(data))]) < 0
+        isnegative = np.nanmax(data[np.isfinite(data)]) < 0
         if isnegative:
             self.log.debug("[is a negative island]")
 
@@ -2484,10 +2488,14 @@ def find_sources_in_image(
 
         for island in islands:
             [[xmin, xmax], [ymin, ymax]] = island.bounding_box
-            i = global_data.data_pix[xmin:xmax, ymin:ymax]
+            island_mask = island.mask
+
+            i = copy.deepcopy(global_data.data_pix[xmin:xmax, ymin:ymax])
+            i[island_mask] = np.nan
+
             # ignore empty islands
             # This should now be impossible to trigger
-            if np.size(i) < 1:
+            if not np.any(np.isfinite(i)):
                 self.log.warning(
                     "Empty island detected, this should be imposisble.")
                 continue