revert to old masking

deprecated/arc_old.py

@@ -855,3 +855,53 @@ def saturation_mask(a, satlevel):
 
     return mask.astype(int)
 
+def mask_around_peaks(spec, inbpm):
+    """
+    Find peaks in the input spectrum and mask pixels around them.
+
+    All pixels to the left and right of a peak is masked until
+    a pixel has a lower value than the adjacent pixel. At this
+    point, we assume that spec has reached the noise level.
+
+    Parameters
+    ----------
+    spec: `numpy.ndarray`_
+        Spectrum (1D array) in counts
+    inbpm: `numpy.ndarray`_
+        Input bad pixel mask
+
+    Returns
+    -------
+    outbpm: `numpy.ndarray`_
+        Bad pixel mask with pixels around peaks masked
+    """
+    # Find the peak locations
+    pks = detect_peaks(spec)
+
+    # Initialise some useful variables and the output bpm
+    xarray = np.arange(spec.size)
+    specdiff = np.append(np.diff(spec), 0.0)
+    outbpm = inbpm.copy()
+
+    # Loop over the peaks and mask pixels around them
+    for i in range(len(pks)):
+        # Find all pixels to the left of the peak that are above the noise level
+        wl = np.where((xarray <= pks[i]) & (specdiff > 0.0))[0]
+        ww = (pks[i]-wl)[::-1]
+        # Find the first pixel to the left of the peak that is below the noise level
+        nmask = np.where(np.diff(ww) > 1)[0]
+        if nmask.size != 0 and nmask[0] > 5:
+            # Mask all pixels to the left of the peak
+            mini = max(0,wl.size-nmask[0]-1)
+            outbpm[wl[mini]:pks[i]] = True
+        # Find all pixels to the right of the peak that are above the noise level
+        ww = np.where((xarray >= pks[i]) & (specdiff < 0.0))[0]
+        # Find the first pixel to the right of the peak that is below the noise level
+        nmask = np.where(np.diff(ww) > 1)[0]
+        if nmask.size != 0 and nmask[0] > 5:
+            # Mask all pixels to the right of the peak
+            maxi = min(nmask[0], ww.size)
+            outbpm[pks[i]:ww[maxi]+2] = True
+    # Return the output bpm
+    return outbpm
+

pypeit/core/arc.py

@@ -1070,57 +1070,6 @@ def detect_lines(censpec, sigdetect=5.0, fwhm=4.0, fit_frac_fwhm=1.25, input_thr
     return tampl_true, tampl, tcent, twid, centerr, ww, arc, nsig
 
 
-def mask_around_peaks(spec, inbpm):
-    """
-    Find peaks in the input spectrum and mask pixels around them.
-
-    All pixels to the left and right of a peak is masked until
-    a pixel has a lower value than the adjacent pixel. At this
-    point, we assume that spec has reached the noise level.
-
-    Parameters
-    ----------
-    spec: `numpy.ndarray`_
-        Spectrum (1D array) in counts
-    inbpm: `numpy.ndarray`_
-        Input bad pixel mask
-
-    Returns
-    -------
-    outbpm: `numpy.ndarray`_
-        Bad pixel mask with pixels around peaks masked
-    """
-    # Find the peak locations
-    pks = detect_peaks(spec)
-
-    # Initialise some useful variables and the output bpm
-    xarray = np.arange(spec.size)
-    specdiff = np.append(np.diff(spec), 0.0)
-    outbpm = inbpm.copy()
-
-    # Loop over the peaks and mask pixels around them
-    for i in range(len(pks)):
-        # Find all pixels to the left of the peak that are above the noise level
-        wl = np.where((xarray <= pks[i]) & (specdiff > 0.0))[0]
-        ww = (pks[i]-wl)[::-1]
-        # Find the first pixel to the left of the peak that is below the noise level
-        nmask = np.where(np.diff(ww) > 1)[0]
-        if nmask.size != 0 and nmask[0] > 5:
-            # Mask all pixels to the left of the peak
-            mini = max(0,wl.size-nmask[0]-1)
-            outbpm[wl[mini]:pks[i]] = True
-        # Find all pixels to the right of the peak that are above the noise level
-        ww = np.where((xarray >= pks[i]) & (specdiff < 0.0))[0]
-        # Find the first pixel to the right of the peak that is below the noise level
-        nmask = np.where(np.diff(ww) > 1)[0]
-        if nmask.size != 0 and nmask[0] > 5:
-            # Mask all pixels to the right of the peak
-            maxi = min(nmask[0], ww.size)
-            outbpm[pks[i]:ww[maxi]+2] = True
-    # Return the output bpm
-    return outbpm
-
-
 def fit_arcspec(xarray, yarray, pixt, fitp):
     """
     Fit a series of pre-identified arc spectrum lines.

pypeit/core/findobj_skymask.py

@@ -1763,16 +1763,8 @@ def objs_in_slit(image, ivar, thismask, slit_left, slit_righ,
     gpm_smash = npix_smash > 0.3*nsmash
     flux_sum_smash = np.sum((image_rect*gpm_sigclip)[find_min_max_out[0]:find_min_max_out[1]], axis=0)
     flux_smash = flux_sum_smash*gpm_smash/(npix_smash + (npix_smash == 0.0))
-    # Mask around the peaks
-    peak_bpm = arc.mask_around_peaks(flux_smash, np.logical_not(gpm_smash))
-    # Check that the peak mask is not empty
-    if np.sum(np.logical_not(peak_bpm))/peak_bpm.size < 0.05:  # 95% masked!
-        msgs.warn('Peak masking was too aggressive. Attempting to sigma clip.')
-        smash_mask = np.logical_not(gpm_smash)
-    else:
-        smash_mask = np.logical_not(gpm_smash) | peak_bpm
     flux_smash_mean, flux_smash_med, flux_smash_std = astropy.stats.sigma_clipped_stats(
-        flux_smash, mask=smash_mask, sigma_lower=3.0, sigma_upper=3.0)
+        flux_smash, mask=np.logical_not(gpm_smash), sigma_lower=3.0, sigma_upper=3.0)
     flux_smash_recen = flux_smash - flux_smash_med
 
     # Return if none found and no hand extraction

pypeit/tests/test_arc.py

@@ -15,45 +15,5 @@ def test_detect_lines():
             = arc.detect_lines(arx_sky.flux.value)
     assert (len(arx_w) > 3275)
 
-
-def test_mask_around_peaks():
-    # Generate a fake spectrum
-    np.random.seed(0)
-    npks = 5   # Number of peaks to generate
-    npix = 20  # Number of pixels in each segment
-    sigma = 2.0  # Fluctuations
-    linewid = 2.0  # Width of the peaks
-    noiselev = [np.random.normal(0.0, sigma, npix) for i in range(npks+1)]
-    ampl = np.random.randint(100*sigma, 200*sigma, npks)
-    spec = np.array([])
-    outbpm = np.array([], dtype=bool)
-    for i in range(npks):
-        spec = np.append(spec, noiselev[i])
-        outbpm = np.append(outbpm, np.zeros(npix, dtype=bool))
-        gauss = ampl[i]*np.exp(-0.5*(np.arange(npix)-npix/2)**2/linewid**2)
-        spec = np.append(spec, np.random.normal(gauss, sigma))
-        if i == 0 or i == 1:
-            tmpbpm = np.array(2*[False] + 16*[True] + 2*[False])
-        elif i == 2:
-            tmpbpm = np.array(1*[False] + 18*[True] + 1*[False])
-        elif i == 3:
-            tmpbpm = np.array(3*[False] + 15*[True] + 2*[False])
-        elif i == 4:
-            tmpbpm = np.array(20*[True])
-        else:
-            assert (False, "Input test data have changed")
-        outbpm = np.append(outbpm, tmpbpm.copy())
-    outbpm = np.append(outbpm, np.zeros(npix, dtype=bool))
-    spec = np.append(spec, noiselev[npks])
-    inbpm = np.zeros(len(spec), dtype=bool)
-    # This is the expected bpm
-    bpm = arc.mask_around_peaks(spec, inbpm=inbpm)
-    # Check that all peaks were identified and masked
-    ind = arc.detect_peaks(bpm)
-    assert (len(ind) == npks)
-    # Check the BPM matches exactly
-    assert (np.array_equal(bpm, outbpm))
-
-
 # Many more functions in pypeit.core.arc that need tests!