Preprocessing

Scikit-learn transformers for preprocessing soil spectra

Each transformer follows the scikit-learn estimator API, with fit, transform, fit_transform, get_params and set_params. You can chain transformers in a Pipeline and tune their parameters with GridSearchCV. Every transformer takes a 2-D array X with one row per spectrum and one column per wavenumber, and returns a 2-D array.

from sklearn.pipeline import Pipeline
from scipy.stats import iqr
from soilspectfm.visualization import plot_spectra
from matplotlib import pyplot as plt
from soilspectfm.visualization import plot_spectra_comparison
from soilspectfm.datasets import load_toy_mir, load_toy_noisy_mir

The examples use the toy MIR dataset from load_toy_mir: 50 spectra measured from 600 to 4000 cm⁻¹, every 2 cm⁻¹.

X, ws = load_toy_mir()

Baseline corrections

source

SNV

def SNV(
    center_func:Callable=mean, # Function to center the data
    scale_func:Callable=std, # Function to scale the data
    eps:float=1e-10, # Small value to avoid division by zero
):

Center and scale each spectrum independently (Standard Normal Variate)

Exported source 
class SNV(BaseEstimator, TransformerMixin):
    "Center and scale each spectrum independently (Standard Normal Variate)"
    def __init__(self,
        center_func: Callable=np.mean, # Function to center the data
        scale_func: Callable=np.std, # Function to scale the data
        eps: float=1e-10 # Small value to avoid division by zero
    ):
        store_attr()
    def fit(self, X, y=None): return self
    def __sklearn_is_fitted__(self): return True # `fit` learns nothing
    def transform(self,
        X: np.ndarray # Spectral data to be transformed
    ) -> np.ndarray: # Transformed spectra
        center = self.center_func(X, axis=1, keepdims=True)
        scale = self.scale_func(X - center, axis=1, keepdims=True) + self.eps
        return (X - center) / scale

SNV subtracts a center from each spectrum and divides by a spread, both computed across its wavenumbers. This removes baseline offsets and multiplicative scatter effects that differ between samples. SNV needs no fitted statistics: fit does nothing, and each spectrum is corrected on its own.

center_func and scale_func accept any function that takes numpy’s axis and keepdims arguments. Common choices:

Role Function Behavior
Center np.mean (default) Standard choice, sensitive to outliers
Center np.median Less affected by outliers, slower
Center np.min Makes all values non-negative, sensitive to noise
Center lambda x, **kw: 0 No centering, keeps absolute values
Scale np.std (default) Standard choice, assumes a normal distribution
Scale lambda x, **kw: np.sqrt(np.mean(x**2, **kw)) Root mean square, suits baseline variations
Scale scipy.stats.iqr Less affected by outliers and extreme peaks
Scale lambda x, **kw: np.max(x, **kw) - np.min(x, **kw) Range, keeps relative peak heights
Scale lambda x, **kw: np.median(np.abs(x - np.median(x, **kw)), **kw) Median absolute deviation, least affected by outliers, slower

With the defaults, every corrected spectrum has mean 0 and standard deviation 1:

X_tfm = SNV().fit_transform(X)
test_close(X_tfm.mean(axis=1), 0)
test_close(X_tfm.std(axis=1), 1)

A median center with an interquartile-range scale reduces the influence of strong peaks. Each corrected spectrum then has a median of 0:

X_rob = SNV(center_func=np.median, scale_func=iqr).fit_transform(X)
test_close(np.median(X_rob, axis=1), 0)
plot_spectra_comparison(X, X_tfm, ws, raw_title='Raw Spectra', transformed_title='SNV Transformed Spectra');

source

MSC

def MSC(
    reference_method:Union[str, numpy.ndarray]='mean', # 'mean', 'median', or a reference spectrum
):

Multiplicative Scatter Correction against a reference spectrum

Exported source 
class MSC(BaseEstimator, TransformerMixin):
    "Multiplicative Scatter Correction against a reference spectrum"
    def __init__(self,
        reference_method: Union[str, np.ndarray] = 'mean', # 'mean', 'median', or a reference spectrum
    ):
        store_attr()
        
    def _compute_reference(self, x: np.ndarray):
        "Compute reference spectrum from array using specified method"
        if isinstance(self.reference_method, str):
            if self.reference_method not in ('mean', 'median'): raise ValueError("reference_method must be 'mean', 'median' or a spectrum")
            return np.mean(x, axis=0) if self.reference_method == 'mean' else np.median(x, axis=0)
        return np.array(self.reference_method)
    
    def fit(self, X: np.ndarray, y=None):
        "Compute the reference spectrum"
        self.reference_ = self._compute_reference(X)
        return self
    
    def transform(self,
        X: np.ndarray # Spectral data to be transformed
    ) -> np.ndarray: # Transformed spectra
        "Apply MSC to the spectra"
        check_is_fitted(self)
        X = np.asarray(X)
        slope, intercept = np.polyfit(self.reference_, X.T, deg=1)
        return (X - intercept[:, None]) / slope[:, None]

MSC models each spectrum x as a linear function of a reference spectrum, x ≈ a + b·reference. It fits a and b by least squares, then returns (x - a) / b. This removes additive and multiplicative scatter effects relative to the reference.

fit sets reference_ to the mean or median of the training spectra, or to the spectrum you pass as reference_method. transform corrects new spectra against that stored reference, so test spectra are corrected against the training reference. After correction, regressing each spectrum on the reference gives a slope of 1 and an intercept of 0:

msc = MSC().fit(X)
X_msc = msc.transform(X)
slope, intercept = np.polyfit(msc.reference_, X_msc.T, deg=1)
test_close(slope, 1)
test_close(intercept, 0)

reference_method='median' uses the median spectrum as the reference, which is less sensitive to outlying samples:

X_tfm = MSC(reference_method='median').fit_transform(X)
plot_spectra_comparison(X, X_msc, ws, raw_title='Raw Spectra', transformed_title='MSC Transformed Spectra');

Asymmetric least squares (ALS) baseline correction is planned. References:

  • https://eigenvector.com/wp-content/uploads/2020/01/RobustFittingtoBasisFunctionsIII.pdf
  • https://nirpyresearch.com/two-methods-baseline-correction-spectral-data/
  • https://diposit.ub.edu/dspace/bitstream/2445/188026/1/2014_IEEE_Adaptive_MarcoS_postprint.pdf

Savitzky-Golay filter

source

SavitzkyGolay

def SavitzkyGolay(
    window_length:int=15, # Number of points in each filter window, odd
    polyorder:int=3, # Degree of the polynomial fitted in each window
    deriv:int=0, # Derivative order, 0 to smooth only
):

Savitzky-Golay smoothing or derivative of each spectrum

Exported source 
class SavitzkyGolay(BaseEstimator, TransformerMixin):
    "Savitzky-Golay smoothing or derivative of each spectrum"
    def __init__(self,
        window_length:int=15, # Number of points in each filter window, odd
        polyorder:int=3, # Degree of the polynomial fitted in each window
        deriv:int=0 # Derivative order, 0 to smooth only
    ): store_attr()

    def fit(self,
        X:np.ndarray, # Spectra, one per row
        y:np.ndarray=None # Ignored
    ):
        "Check that the filter parameters are consistent"
        if self.window_length % 2 == 0: raise ValueError("window_length must be odd")
        if self.window_length <= self.polyorder: raise ValueError("window_length must be greater than polyorder")
        if self.deriv > self.polyorder: raise ValueError("deriv must be <= polyorder")
        return self
    def __sklearn_is_fitted__(self): return True # `fit` learns nothing

    def transform(self,
        X:np.ndarray # Spectra, one per row
    ) -> np.ndarray: # Filtered spectra
        "Filter each spectrum"
        return savgol_filter(np.asarray(X), self.window_length, self.polyorder, self.deriv)

SavitzkyGolay fits a polynomial of degree polyorder within each window of window_length points. With the default deriv=0, it replaces each point by the fitted value, which smooths the spectrum. With deriv=1 or deriv=2, it returns the derivative of the fitted polynomial instead. Derivatives are taken per wavenumber step, not per cm⁻¹. fit checks that window_length is odd and greater than polyorder, and that deriv is at most polyorder.

Smoothing leaves a polynomial of degree up to polyorder unchanged:

poly = ((ws / 1000) ** 3)[None]
test_close(SavitzkyGolay().fit_transform(poly), poly)

A first derivative removes constant baseline offsets, and a second derivative also removes linear baseline slopes. Adding a constant offset to every spectrum leaves the first derivative unchanged:

sg1 = SavitzkyGolay(deriv=1)
test_close(sg1.fit_transform(X + 0.5), sg1.fit_transform(X))

Derivatives usually follow a scatter correction. In a Pipeline, SNV runs first, then SavitzkyGolay. A pipeline fitted on some spectra transforms others, such as spectra from another instrument. Both transformers work on each spectrum on its own. The result does not depend on the spectra used to fit:

pipe = Pipeline([('snv', SNV()), ('deriv', SavitzkyGolay(deriv=1))])
pipe.fit(X[:40])
X_tfm = pipe.transform(X)
test_close(X_tfm[40:], pipe.fit_transform(X[40:]))
plot_spectra_comparison(X, X_tfm, ws, raw_title='Raw Spectra', transformed_title='SNV + Derivative (1st order) Transformed Spectra');

Smoothing

source

WaveletDenoise

def WaveletDenoise(
    wavelet:str='db6', # Wavelet to use for decomposition
    level:Optional[int]=None, # Decomposition level. If None, maximum level is used
    threshold_mode:str='soft', # Thresholding mode, 'soft' or 'hard'
):

Wavelet denoising of each spectrum

Exported source 
class WaveletDenoise(BaseEstimator, TransformerMixin):
    "Wavelet denoising of each spectrum"
    def __init__(self,
        wavelet:str='db6', # Wavelet to use for decomposition
        level:Optional[int]=None, # Decomposition level. If None, maximum level is used
        threshold_mode:str='soft' # Thresholding mode, 'soft' or 'hard'
    ):
        store_attr()
        
    def _denoise_single(self, spectrum):
        "Denoise a single spectrum"
        coeffs = pywt.wavedec(spectrum, self.wavelet, level=self.level_)
        
        detail_coeffs = np.concatenate(coeffs[1:])
        sigma = np.median(np.abs(detail_coeffs)) / 0.6745
        threshold = sigma * np.sqrt(2 * np.log(len(spectrum)))
        
        new_coeffs = list(coeffs)
        for i in range(1, len(coeffs)):
            new_coeffs[i] = pywt.threshold(coeffs[i], threshold * (1 / 2 ** ((self.level_ - i) / 2)), mode=self.threshold_mode)

        denoised = pywt.waverec(new_coeffs, self.wavelet)
        return denoised[:len(spectrum)]
    
    def fit(self, X, y=None):
        "Set `level_`: `level`, or the deepest level the spectrum length allows"
        self.level_ = ifnone(self.level, pywt.dwt_max_level(np.shape(X)[1], pywt.Wavelet(self.wavelet).dec_len))
        return self
    
    def transform(self, X):
        "Denoise each spectrum"
        X = np.asarray(X)
        X_denoised = np.zeros_like(X)
        for i in range(X.shape[0]): X_denoised[i] = self._denoise_single(X[i])
        return X_denoised

WaveletDenoise decomposes each spectrum with a discrete wavelet transform, shrinks the detail coefficients, and reconstructs the spectrum. The base threshold is the universal threshold σ·√(2·ln n), where n is the spectrum length. σ is the median absolute value of all detail coefficients divided by 0.6745. For the i-th detail band, counted from the coarsest (i = 1) to the finest (i = level_), the threshold is multiplied by 2^(-(level_ - i)/2), so coarser bands get smaller thresholds. threshold_mode selects soft or hard thresholding, as in pywt.threshold.

The smoothing examples use load_toy_noisy_mir: 48 noisy MIR spectra from 599.9 to 3998.4 cm⁻¹, about every 1 cm⁻¹. From here on, X and ws refer to this dataset:

X, ws, _ = load_toy_noisy_mir()

The difference between the raw and denoised spectra shows the noise that WaveletDenoise removes:

X_denoised = WaveletDenoise(wavelet='db6', level=5, threshold_mode='soft').fit_transform(X)
plot_spectra(X - X_denoised, ws, alpha=0.1, title='Raw Spectra - Wavelet Denoised');

With the default level=None, fit sets level_ to the deepest decomposition level that the spectrum length allows:

wd = WaveletDenoise().fit(X)
test_eq(wd.level_, pywt.dwt_max_level(X.shape[1], 'db6'))
wd.level_
8

SavitzkyGolay with its defaults also smooths the noisy spectra:

X_smoothed = SavitzkyGolay().fit_transform(X)
plot_spectra_comparison(X, X_smoothed, ws, raw_title='Raw Spectra', transformed_title='Savitzky-Golay Smoothed Spectra');

Other transformations

source

ToAbsorbance

def ToAbsorbance(
    eps:float=1e-05, # Small value to avoid log(0)
):

Convert reflectance to absorbance, -log10(R)

Exported source 
class ToAbsorbance(BaseEstimator, TransformerMixin):
    "Convert reflectance to absorbance, `-log10(R)`"
    def __init__(self,
        eps: float=1e-5 # Small value to avoid log(0)
    ): self.eps = eps
    def fit(self, X, y=None): return self
    def __sklearn_is_fitted__(self): return True # `fit` learns nothing
    def transform(self, X, y=None): return -np.log10(np.clip(X, self.eps, 1))

ToAbsorbance converts reflectance R to absorbance A = -log10(R). It first clips R to the range [eps, 1], so zero or negative values don’t produce infinities, and values above 1 give an absorbance of 0:

test_close(ToAbsorbance().fit_transform(np.array([[1, 0.1, 0.01]])), [[0, 1, 2]])

source

Resample

def Resample(
    source_x:numpy.ndarray, # Source x-axis points (wavenumbers or wavelengths)
    target_x:numpy.ndarray, # Target x-axis points (wavenumbers or wavelengths) for resampling
    interpolation_kind:str='cubic', # Spline kind: 'linear', 'quadratic' or 'cubic'
):

Resample spectra onto new x-axis points by spline interpolation

Exported source 
class Resample(BaseEstimator, TransformerMixin):
    "Resample spectra onto new x-axis points by spline interpolation"
    def __init__(self,
        source_x: np.ndarray, # Source x-axis points (wavenumbers or wavelengths)
        target_x: np.ndarray, # Target x-axis points (wavenumbers or wavelengths) for resampling
        interpolation_kind: str='cubic' # Spline kind: 'linear', 'quadratic' or 'cubic'
    ):
        store_attr()
        
    def fit(self,
        X: np.ndarray, # Spectral data to be resampled
        y: np.ndarray=None # Ignored
    ):
        "Check `X` against `source_x`, check `target_x` is in range, and sort `source_x`"
        n_cols, n_src = np.shape(X)[1], len(self.source_x)
        if n_cols != n_src: raise ValueError(f"X has {n_cols} columns but source_x has {n_src} points")
        self.order_ = np.argsort(self.source_x)
        lo, hi = np.asarray(self.source_x)[self.order_[[0, -1]]]
        if np.min(self.target_x) < lo or np.max(self.target_x) > hi:
            raise ValueError(f"target_x must lie within the source_x range [{lo}, {hi}]")
        return self
    
    def transform(self,
        X: np.ndarray # Spectral data to be resampled
    ):
        "Interpolate each spectrum at `target_x`"
        k = dict(linear=1, quadratic=2, cubic=3)[self.interpolation_kind]
        xs = np.asarray(self.source_x)[self.order_]
        return make_interp_spline(xs, np.asarray(X)[:, self.order_], k=k, axis=1)(self.target_x)

Resample interpolates each spectrum from source_x onto target_x, for example to put spectra from different instruments on one wavenumber grid. Here the noisy spectra, sampled about every 1 cm⁻¹, are resampled every 2 cm⁻¹:

new_ws = np.arange(600, 3998, 2)
X_resampled = Resample(source_x=ws, target_x=new_ws).fit_transform(X)
test_eq(X_resampled.shape, (len(X), len(new_ws)))
plot_spectra(X, ws, ascending=False, title='Original Spectra', alpha=0.2);
plot_spectra(X_resampled, new_ws, ascending=False, title='Resampled Spectra', alpha=0.2, color='steelblue');

interpolation_kind sets the spline degree: 'linear', 'quadratic' or 'cubic'. Linear resampling gives the same result as np.interp:

X_lin = Resample(ws, new_ws, interpolation_kind='linear').fit_transform(X)
test_close(X_lin, np.array([np.interp(new_ws, ws, x) for x in X]))

fit raises a ValueError when target_x goes beyond the range of source_x, instead of extrapolating. The noisy spectra end at 3998.4 cm⁻¹, so a grid up to 4000 cm⁻¹ fails:

test_fail(lambda: Resample(ws, np.arange(600, 4001, 2)).fit(X), contains='source_x range')

source_x can be in decreasing order, as some instruments export spectra from high to low wavenumbers. fit sorts it, and reorders the columns of X to match:

test_close(Resample(ws[::-1], new_ws).fit_transform(X[:, ::-1]), X_resampled)

source

Trim

def Trim(
    ws:numpy.ndarray, # Wavenumbers or wavelengths
    w_min:Optional[float]=None, # Minimum wavenumber or wavelength
    w_max:Optional[float]=None, # Maximum wavenumber or wavelength
):

Keep the wavenumbers (or wavelengths) between w_min and w_max

Exported source 
class Trim(BaseEstimator, TransformerMixin):
    "Keep the wavenumbers (or wavelengths) between `w_min` and `w_max`"
    def __init__(self,
        ws: np.ndarray, # Wavenumbers or wavelengths
        w_min: Optional[float]=None, # Minimum wavenumber or wavelength
        w_max: Optional[float]=None # Maximum wavenumber or wavelength
    ):
        store_attr()
        
    def fit(self,
        X: np.ndarray, # Spectra to be trimmed
        y: Optional[np.ndarray]=None # Ignored
    ):
        "Compute `mask_`, the columns whose wavenumber is inside the range"
        self.mask_ = (self.ws >= (-np.inf if self.w_min is None else self.w_min)) & \
                     (self.ws <= (np.inf if self.w_max is None else self.w_max))
        return self
    
    def transform(self, X: np.ndarray) -> np.ndarray:
        "Keep the columns inside the range"
        return X[:, self.mask_]
    
    def get_wavenumbers(self) -> np.ndarray:
        "Return the trimmed wavenumbers"
        return self.ws[self.mask_]

Trim keeps the wavenumbers between w_min and w_max, inclusive. A bound left as None leaves that side open. get_wavenumbers returns the kept wavenumbers, to use as the x-axis of the trimmed spectra:

trimmer = Trim(ws, w_min=650, w_max=4000)
X_trimmed = trimmer.fit_transform(X)
trimmed_ws = trimmer.get_wavenumbers()
assert trimmed_ws.min() >= 650 and trimmed_ws.max() <= 4000
plot_spectra(X_trimmed, trimmed_ws, title='Trimmed Spectra', ascending=False);