Remove outlier data points from an EELS spectrum by doing a LOWESS fit and then keeping only spots within a specific median sigma.

EELSFuncs.py

%pylab inline
import hyperspy.api as hs
from statsmodels.nonparametric.smoothers_lowess import lowess

def DespikeEELS(spec, m=3., interp_outliers=True):
    # Make a copy of the data
    E = spec.axes_manager[0].axis.copy()
    I = spec.data.copy()
    
    # Lowess smooth.
    smooth = lowess(I, E, frac=0.02)
    
    # Using median, find inliers and outliers as # sigma away.
    dist = np.abs(I-smooth[:,1])
    d = np.abs(dist - np.median(dist))
    mdev = np.median(d)
    s = d/mdev if mdev else 0.
    
    # Which points are inliers and which are out?
    inliers = np.where(s<=m)[0]
    outliers = np.where(s>m)[0]
    
    if interp_outliers == False:
        # Keep the inliers only (drop outliers).
        E = E[inliers]
        I = I[inliers]
        spec2 = spec.deepcopy()
        spec2.axes_manager[0].axis = E
        spec2.data = I
        return spec2
    else:
        # Interpolate the outliers to the lowess curve (retains same spectrum dimensions).
        E[outliers] = smooth[outliers,0]
        I[outliers] = smooth[outliers,1]
        spec2 = spec.deepcopy()
        spec2.axes_manager[0].axis = E
        spec2.data = I
        return spec2

def ProcessEELS(CoreLoss=None, LowLoss=None, BkgRange=None, DespikeThreshhold=None, NormRange=None):
    # Load the core and low loss spectra.
    cl = hs.load(CoreLoss)
    if len(cl) > 1:
        cl = cl[0]
        
    if LowLoss is not None:
        ll = hs.load(LowLoss)
        # Align to ZLP.
        ll.align_zero_loss_peak(also_align=[cl], subpixel=True)
        # Estimate the thickness.
        ll.estimate_thickness(threshold=5.).data
    
    if DespikeThreshhold is not None:
        cl = DespikeEELS(cl, m=DespikeThreshhold)
    
    if BkgRange is not None:
        x = cl.remove_background(background_type='PowerLaw',signal_range=BkgRange, fast=False)
        cl.data = x
        
    if LowLoss is not None:
        deconv = cl.fourier_ratio_deconvolution(ll)

    if NormRange is not None:
        cl.data -= np.min(cl.data)
        cl.data /= np.mean(cl.data[(cl.axes_manager[0].axis > NormRange[0]) & (cl.axes_manager[0].axis < NormRange[1])])
        deconv.data -= np.min(deconv.data)
        deconv.data /= np.mean(deconv.data[(deconv.axes_manager[0].axis > NormRange[0]) & (deconv.axes_manager[0].axis < NormRange[1])])


    return ll, cl, deconv