Source code for mcalf.utils.spec

import numpy as np
import scipy.interpolate

from mcalf.profiles.gaussian import single_gaussian

__all__ = ['reinterpolate_spectrum', 'normalise_spectrum', 'generate_sigma']

[docs]def reinterpolate_spectrum(spectrum, original_wavelengths, constant_wavelengths): """Reinterpolate the spectrum. Reinterpolates the spectrum such that intensities at `original_wavelengths` are transformed into intensities at `constant_wavelengths`. Uses :class:`scipy.interpolate.InterpolatedUnivariateSpline` to interpolate. Parameters ---------- spectrum : numpy.ndarray, ndim=1 Spectrum to reinterpolate. original_wavelengths : numpy.ndarray, ndim=1, length=length of `spectrum` Wavelengths of `spectrum`. constant_wavelengths : numpy.ndarray, ndim=1 Wavelengths to cast `spectrum` into. Returns ------- spectrum : numpy.ndarray, length=length of `constant_wavelengths` Reinterpolated spectrum. """ s = scipy.interpolate.InterpolatedUnivariateSpline(original_wavelengths, spectrum) # Fit spline return s(constant_wavelengths) # Evaluate spline at new wavelengths
[docs]def normalise_spectrum(spectrum, original_wavelengths=None, constant_wavelengths=None, prefilter_response=None, model=None): """Normalise an individual spectrum to have intensities in range [0, 1]. .. warning:: Not recommended for normalising many spectra in a loop. Parameters ---------- spectrum : numpy.ndarray, ndim=1 Spectrum to reinterpolate and normalise. original_wavelengths : numpy.ndarray, ndim=1, length=length of `spectrum`, optional Wavelengths of `spectrum`. constant_wavelengths : numpy.ndarray, ndim=1, optional Wavelengths to cast `spectrum` into. prefilter_response : numpy.ndarray, ndim=1, length=length of `constant_wavelengths`, optional Prefilter response to divide spectrum by. model : child class of mcalf.models.ModelBase, optional Model to extract the above parameters from. Returns ------- spectrum : numpy.ndarray, ndim-1, length=length of `constant_wavelengths` The normalised spectrum. """ from mcalf.models import ModelBase if issubclass(model.__class__, ModelBase): if original_wavelengths is None: original_wavelengths = model.original_wavelengths if constant_wavelengths is None: constant_wavelengths = model.constant_wavelengths if prefilter_response is None: prefilter_response = model.prefilter_response if original_wavelengths is not None or constant_wavelengths is not None: if original_wavelengths is None or constant_wavelengths is None: raise ValueError("original_wavelengths and constant_wavelengths must be provided") spectrum = reinterpolate_spectrum(spectrum, original_wavelengths, constant_wavelengths) if prefilter_response is not None: spectrum /= prefilter_response spectrum -= min(spectrum) spectrum /= max(spectrum) return spectrum
[docs]def generate_sigma(sigma_type, wavelengths, line_core, a=-0.95, c=0.04, d=1, centre_rad=7, a_peak=0.4): """Generate the default sigma profiles. Parameters ---------- sigma_type : int Type of profile to generate. Should be either `1` or `2`. wavelengths : array_like Wavelengths to use for sigma profile. line_core : float Line core to use as centre of Gaussian sigma profile. a : float, optional, default=-0.95 Amplitude of Gaussian sigma profile. c : float, optional, default=0.04 Sigma of Gaussian sigma profile. d : float, optional, default=1 Background of Gaussian sigma profile. centre_rad : int, optional, default=7 Width of central flattened region. a_peak : float, optional, default=0.4 Amplitude of central 7 pixel section, if `sigma_type` is 2. Returns ------- sigma : numpy.ndarray, length=n_wavelengths The generated sigma profile. """ sigma = single_gaussian(wavelengths, a, line_core, c, d) # Define initial sigma profile centre_i = np.argmin(sigma) # Get central index of sigma sigma[:centre_i - centre_rad] = sigma[centre_rad:centre_i] # Shift left half left sigma[centre_i + centre_rad:] = sigma[centre_i:-centre_rad] # Shift right half right sigma[centre_i - centre_rad:centre_i + centre_rad] = sigma[centre_i] # Connect the two halves with a line if sigma_type == 2: sigma[centre_i - 3:centre_i + 3 + 1] = a_peak # Lessen importance of peak detail return sigma