mangadap.util.resolution module¶
Provides a set of functions to handle and alter the instrumental resolution.
Revision history¶
27 May 2015: Original implementation by K. Westfall (KBW) based on downgrader_MANGA.f provided by D. Thomas, O. Steele, D. Wilkinson, D. Goddard.13 Jun 2015: D.Wilkinson edit to not calculate unimportant convolution terms > runs 5x faster.03 Feb 2016: (KBW) Generalized M. Cappellari’slog_rebin()
toresample_vector()
to allow for linear interpolation of flux density across the pixel.04 Feb 2016: (KBW) Further fixes toresample_vector()
. Addedspectral_coordinate_step()
, and propagated change tospectrum_velocity_scale()
21 Apr 2016: (KBW) It’s the Queen’s 90th birthday! Removed log10 keyword fromspectrum_velocity_scale()
.20 May 2016: (KBW) Corrected match between number of pixels and output range computed inresample_vector_npix()
; now returns an adjusted range to make sure that the sampling and range results in an exact integer number of pixels.05 Jul 2016: (KBW) To avoid confusion, commented out log_rebin.25 Oct 2016: (KBW) Modifiedspectral_coordinate_step()
to be a mean over the full spectrum to avoid numerical precision errors.06 Apr 2017: (KBW) Addangstroms_per_pixel()
.30 Aug 2017: (KBW) Addresample1d()
;resample_vector()
should be deprecated.27 Sep 2017: (KBW) Added integral keyword tomatch_spectral_resolution()
so that it can be passed toconvolution_variable_sigma()
.19 Jul 2018: (KBW) Fixed a bug inVariableGaussianKernel
when constructing the kernel based on the error function.30 Aug 2018: (KBW) Removed sampling code and moved it to sampling.py. Renamed this file from instrument.py to resolution.py.
Copyright © 2019, SDSSIV/MaNGA Pipeline Group

class
mangadap.util.resolution.
SpectralResolution
(wave, sres, log10=False, interp_ext='extrapolate')[source]¶ Bases:
object
Container class for the resolution, \(R = \lambda/\Delta\lambda\), of a spectrum. The primary functionality is to determine the parameters necessary to match the resolution of one spectrum to another. It can also be used as a function to interpolate the spectral resolution at a given wavelenth.
Parameters:  wave (numpy.ndarray) – A 1D vector with the wavelength in angstroms. The sampling may be either in linear steps of wavelength or \(\log_{10}\) steps.
 sres (numpy.ndarray) – A 1D vector with the spectral resolution, \(R\), sampled at the positions of the provided wavelength vector.
 log10 (
bool
, optional) – Flag that the spectrum has been binned logarithmically (base 10) in wavelength  interp_ext (
int
,str
, optional) – The value to pass as ext to the interpolator, which defines its behavior when attempting to sample the spectral resolution beyond where it is defined. See scipy.interpolate.interp1d. Default is to extrapolate.
Raises: ValueError
– Raised if wave is not a 1D vector or if wave and sres do not have the same shape.
interpolator
¶ An object used to interpolate the spectral resolution at any given wavelength. The interpolation is hardwired to be linear and its extrapolation behavior is defined by
interp_ext
. The wavelength and resolution vectors are held by this object for later reference if needed.Type: scipy.interpolate.interp1d

log10
¶ Flag that the spectrum has been binned logarithmically (base 10) in wavelength
Type: bool

c
¶ The speed of light; provided by astropy.constants.
Type: float

dv
¶ The velocity step per pixel in km/s. Defined using
spectrum_velocity_scale()
iflog10
is True; otherwise set to None.Type: float

dw
¶ The wavelength step per pixel in angstroms. Defined as the wavelength step between the first two pixels if
log10
is False; otherwise set to None.Type: float

min_sig
¶ Minimum standard deviation allowed for the kernel used to match two spectral resolutions. See
GaussianKernelDifference()
.Type: float

sig_pd
¶ The standard deviation in pixels required to match the spectral resolution of this object to the resolution defined by a different
SpectralResolution
object. SeeGaussianKernelDifference()
.Type: numpy.ndarray

sig_mask
¶ A uint vector used to identify measurements of
sig_pd
that should not be used to match the spectral resolution of this object to the resolution defined by a differentSpectralResolution
object. SeeGaussianKernelDifference()
.Type: numpy.ndarray

sig_vo
¶ A constant offset of the kernal standard deviation in km/s that has been applied to
sig_pd
. SeeGaussianKernelDifference()
.Type: float
Todo
 Allow it to determine if the binning is linear or geometric, then use the log10 keyword to distinguish between natural log and \(log_{10}\) binning.
 Allow for more than one type of linespread function?
Warning
The default behavior of the interpolator is to extrapolate the input spectral resolution vector when trying to sample from regions beyond where it is sampled. Use
interp_ext
change this; see scipy.interpolate.interp1d.
GaussianKernelDifference
(new_sres, no_offset=True, min_sig_pix=0.0)[source]¶ Determine the parameters of a Gaussian kernel required to convert the resolution of this object to the resolution of a different the
SpectralResolution
object, new_sres.The spectral resolution is defined as \(R = \lambda/\Delta\lambda\), where \(\Delta\lambda\) is the FWHM of the spectral resolution element. The standard deviation of the resolution element in angstroms is then
\[\sigma_\lambda = \frac{\lambda}{f R}, \ \ {\rm where} \ \ f = \frac{{\rm FWHM_\lambda}}{\sigma_\lambda} = \sqrt{8\ln 2}.\]Assuming a Gaussian (in angstroms) linespread function:
\[\sigma^2_{\lambda,2} = \sigma^2_{\lambda,1} + \sigma^2_{\lambda,d}\]such that
\[\sigma^2_{\lambda,d} = \left(\frac{\lambda}{f}\right)^2 (R^{2}_2  R^{2}_1)\]is the defining parameter of the Gaussian kernel needed to take a spectrum of resolution \(R_1\) to one with a resolution of \(R_2\).
For input to
convolution_variable_sigma()
, the wavelengthdependent parameter of the Gaussian kernel is converted to pixels. This function allows for the input spectra to be linearly sampled in angstroms or log10(angstroms). For the former (log10=False),\[\sigma^2_{p,d} = \left(\frac{\lambda}{f\ \delta\lambda}\right)^2 (R^{2}_2  R^{2}_1)\]where \(\delta\lambda\) is the size of the pixel in angstroms. If the units are log10(angstrom) (log10=True), we approximate the velocity width of each pixel to be \(\delta v = c \ln(10.0) (\log\lambda[1]\log\lambda[0])\), such that
\[\begin{split}\sigma^2_{p,d} &= \left(\frac{c}{ \delta v \lambda}\right)^2 \sigma^2_{\lambda,d} \\ &= \left(\frac{c}{ f\ \delta v}\right)^2 (R^{2}_2  R^{2}_1)\ ;\end{split}\]\(c\) is the speed of light in km/s.
The nominal use of this algorithm assumes \(R_1 \geq R_2\). However, in practice,
convolution_variable_sigma()
only uses a Gaussian kernel up to some minimum value of \(\epsilon_\sigma\); below this, the kernel is assumed to be a Delta function. Therefore, as long as\[\sigma_{p,d} \equiv \sigma^2_{p,d}/\sqrt{\sigma^2_{p,d}} \geq \epsilon_\sigma\ ,\]the behavior of
convolution_variable_sigma()
should not be affected. However, in regions with\[\sigma_{p,d} \equiv \sigma^2_{p,d}/\sqrt{\sigma^2_{p,d}} \leq \epsilon_\sigma\ ,\]the behavior of
convolution_variable_sigma()
does not produce an accurate convolution!To deal with spectral regions that do not have \(\sigma_{p,d} \geq \epsilon_\sigma\), there are three choices:
(Option 1) trim the spectral range to only those spectral regions where the existing resolution is better than the target resolution up to this limit,
(Option 2) match the existing resolution to the target resolution up to some constant offset that must be accounted for in subsequent analyses, or
(Option 3) allow for a wavelength dependent difference in the spectral resolution that must be accounted for in subsequent analyses.
The choice of either Option 1 or 2 is selected by setting no_offset to, respectively, True or False; Option 1 is the default behavior. Currently, Option 3 is not allowed.
For Option 1, pixels with \(\sigma_{p,d} < \epsilon_\sigma\) are masked (sigma_mask = 1); however, the returned values of \(\sigma_{p,d}\) are left unchanged.
For Option 2, we define
\[\sigma^2_{v,o} =  {\rm min}(0.0, {\rm min}(\sigma^2_{v,d})  {\rm max}(\epsilon_\sigma \delta v)^2) \]where \(\delta v\) is constant for the logarithmically binned spectrum and is wavelength dependent for the linearly binned spectra; in the latter case, the velocity step is determined for each pixel:
_wave = self.wave() dv = self.c * (2.0*(_wave[1:]  _wave[0:1]) / (_wave[1:] + _wave[0:1]))
If \(\sigma^2_{v,o} > 0.0\), it must be that \({\rm min}(\sigma^2_{v,d}) < {\rm max}(\epsilon_\sigma \delta v)^2\), such that an offset should be applied. In that case, the returned kernel parameters are
\[\sigma^\prime_{v,d} = \sqrt{\sigma^2_{v,d} + \sigma^2_{v,o}}\ .\]with the units converted to pixels using the equations above, no pixels are masked, and \(\sqrt{\sigma^2_{v,o}}\) is returned for the offset. Otherwise, the offset is set to 0.
Parameters:  new_sres (
SpectralResolution
) – Spectral resolution to match to.  no_offset (
bool
, optional) – Force \(\sigma^2_{v,o} = 0\) by masking regions with \(\sigma_{p,d} < \epsilon_\sigma\); i.e., the value of this arguments selects Option 1 (True) or Option 2 (False).  min_sig_pix (
float
, optional) – Minimum value of the standard deviation allowed before assuming the kernel is a Delta function.
 new_sres (

_convert_pd2vd
(sig2_pd)[source]¶ Convert from \(\sigma^2_{p,d}\) to \(\sigma^2_{v,d}\). See
GaussianKernelDifference()
.

_convert_vd2pd
(sig2_vd)[source]¶ Convert from \(\sigma^2_{v,d}\) to \(\sigma^2_{p,d}\). See
GaussianKernelDifference()
.

_finalize_GaussianKernelDifference
(sig2_pd)[source]¶ Given the calculated \(\sigma^2_{p,d}\), calculate and save the attributes
sig_pd
andsig_mask
. SeeGaussianKernelDifference()
.

adjusted_resolution
(indx=None)[source]¶ Return the resolution that should result from applying
convolution_variable_sigma()
to the spectrum associated with this spectral resolution object usingsigma_pd
. I.e., calculate:\[R_2 = \left[ \left(\frac{f}{c}\right)^2 \sigma^2_{v,d} + R^{2}_1\right]^{1/2}\ .\]Parameters: indx (tuple) – (Optional) Selection tuple used to return a subset of the full resolution vector. Returns: The (full or selected) vector with the adjusted resolution. Return type: numpy.ndarray Todo
Allow to reset the resolution of this object to the adjusted resolution and reset the kernel variables to None.

instrumental_dispersion
(w=None)[source]¶ Return the instrumental dispersion by converting from \(R\) to \(\sigma_{v,inst}\) according to:
\[\sigma_{v,inst} = \frac{c}{\sqrt{8\ln 2}\ R}.\]If w is None, just convert the internal interpolator values. Otherwise, return the values sampled at w.

match
(new_sres, no_offset=True, min_sig_pix=0.0)[source]¶ Currently only an alias for
GaussianKernelDifference()
.

offset_GaussianKernelDifference
(offset)[source]¶ If the properties required to match the resolution of one spectrum to another has already been calculated (see
GaussianKernelDifference()
), this allows for one to apply an additional offset. The additional offset must be in km/s (not pixels).The offset is applied in quadrature; however, the offset can be negative such that one can reduce
sig_pd
. Once converted to km/s, the offset is applied by calculating:\[\sigma^{\prime\ 2}_{v,d} = \sigma^{2}_{v,d} + \sigma_{off}\sigma_{off}\ .\]sig_vo
is adjusted in the same way, and the change in \(\sigma^{\prime\ 2}_{v,d}\) is then propagated tosig_pd
andsig_mask
.Parameters: offset (float) – Value of the standard deviation in km/s to add in quadrature to a previously calculated sig_pd
.Raises: ValueError
– Raised if the kernel properties have not yet been defined.

sres
(copy=False)[source]¶ Return the resolution vector; held by
interpolator
.

wave
(copy=False)[source]¶ Return the wavelength vector; held by
interpolator
.

class
mangadap.util.resolution.
VariableGaussianKernel
(sigma, minsig=0.01, nsig=3.0, integral=False)[source]¶ Bases:
object
Support class for variable sigma convolution. See
convolution_variable_sigma()
.Modest modification of M. Cappellari’s
ppxf.ppxf_util.gaussian_filter1d()
function.Parameters:  sigma (numpy.ndarray) – Coordinatedependent standard deviation of the Gaussian kernel in pixels.
 minsig (
float
, optional) – The minimum sigma value allowed.  nsig (
float
, optional) – The extent of the calculation of the kernel for the convolution. By default, the kernel is not calculated beyond 3*sigma.  integral (
bool
, optional) – Construct the kernel based on the integral of the Gaussian over the pixel width using the error function, instead of just sampling the Gaussian function directly.
Raises: ValueError
– Raised if the input sigma has no length.
sigma
¶ Coordinatedependent standard deviation of the Gaussian kernel.
Type: numpy.ndarray

p
¶ The maximum number of pixels on either side of the kernel center needed for the convolution.
Type: int

m
¶ Number of pixels in the kernel (\(m = 2p + 1\)).
Type: int

kernel
¶ The calculated kernel to use for the convolution.
Type: numpy.ndarray

convolve
(y, ye=None)[source]¶ Convolve the input vector with the variablewidth Gaussian.
The input y must be a 1D vector with the same length as
sigma
. If errors are provided, they must also adhere to these shape limitations and are propagated to the convolved vector.Parameters:  y (numpy.ndarray) – Vector to convolve
 ye (numpy.ndarray, optional) – Error in the vector to convolve
Raises: ValueError
– Raised if y is not a 1D vector, or if the shape of y and sigma (and ye if provided) are different.Returns: The convovled vector. If the errors are provided, the function returns two vectors, the convolved vector and its error.
Return type: numpy.ndarray

mangadap.util.resolution.
convolution_variable_sigma
(y, sigma, ye=None, integral=False)[source]¶ Convolve a discretely sampled function \(y(x)\) with a Gaussian kernel, \(g\), where the standard deviation of the kernel is a function of \(x\), \(\sigma(x)\). Nominal calculations can be performed to propagate the error in the result; these calculations do not include the covariance between the pixels, which will mean that the calculations likely have significant error!
The convolution is defined as:
\[\begin{split}(y\ast g)(x) &= \int_{\infty}^{\infty} y(X)\ g(\sigma,xX)\ dX \\ &= \int_{\infty}^{\infty} \frac{y(X)}{\sqrt{2\pi}\ \sigma(X)}\ \exp\left(\frac{(xX)^2}{2\ \sigma(X)^2}\right) dX .\end{split}\]To minimize edge effects and account for the censoring of the data (finite range in \(x\)), the convolution is actually calculated as a definite integral and normalized as follows:
\[(y\ast g)(x) \sim\frac{ \int_{xn_\sigma*\sigma(x)}^{x+n_\sigma*\sigma(x)} y(X)\ g(\sigma,xX)\ dX}{ \int_{xn_\sigma*\sigma(x)}^{x+n_\sigma*\sigma(x)} g(\sigma,xX)\ dX} .\]The above is identical to getting the weighted mean of \(y\) at each position \(x\), where the weights are defined by a Gaussian kernel centered at \(x\) with a variable dispersion.
 Use of this function requires:
 y and sigma must be 1D vectors
 y and sigma must be uniformly sampled on the same grid
 sigma must be in pixel units.
Todo
This should be deprecated and/or merged with
VariableGaussianKernel
.Parameters:  y (numpy.ndarray) – A uniformly sampled function to convolve.
 sigma (numpy.ndarray) – The standard deviation of the Gaussian kernel sampled at the same positions as y. The units of sigma must be in pixels.
 ye (numpy.ndarray) – (Optional) Errors in the function \(y(x)\).
Returns: Arrays with the convolved function \((y\ast g)(x)\) sampled at the same positions as the input \(x\) vector and its error. The second array will be returned as None if the error vector is not provided.
Return type: numpy.ndarray

mangadap.util.resolution.
match_spectral_resolution
(wave, flux, sres, new_sres_wave, new_sres, ivar=None, mask=None, min_sig_pix=0.0, no_offset=True, variable_offset=False, log10=False, new_log10=False, integral=True, quiet=False)[source]¶ Adjust the existing spectral resolution of a spectrum to a lower resolution as best as possible. The primary functionality is in
SpectralResolution
, which determines the Gaussian kernel parameters needed to match the resolution, andconvolve_variable_sigma()
, which actually performs the convolution to match the resolution.In particular, see
SpectralResolution.GaussianKernelDifference()
for a description of how the kernel parameters are determined and how regions where the target resolution is higher than the existing resolution. In this case, one of the options is to adopt an offset of the resolution (in km/s) that could be corrected for in subsequent analysis. In this case, setting variable_offset to True allows the offset to be different for all input spectra. If one expects to combine the spectra, the default behavior should be used, forcing all the spectra to have a constant offset.Parameters:  wave (numpy.ndarray) – A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the wavelength in angstroms for a set of spectra. The sampling may be either in linear steps of wavelength or \(\log_{10}\) steps, as set using log10.
 flux (numpy.ndarray) – A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the flux sampled at the provided wavelengths.
 sres (numpy.ndarray) – A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the spectral resolution, \(R\), at the provided wavelengths.
 new_sres_wave (numpy.ndarray) – A 1D vector with the wavelength in angstroms at which the new resolution of the input spectra has been sampled. The sampling may be either in linear steps of wavelength or \(\log_{10}\) steps, as set using new_log10.
 new_sres (numpy.ndarray) – A 1D vector with the new resolution for the input spectra.
 ivar (numpy.ndarray) –
(Optional) A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the inverse variance of the flux sampled at the provided wavelengths. This vector is used to estimate the noise in the resolutionmatched spectra.
Warning
The accuracy of the errors still remain untested!
 mask (numpy.ndarray) – (Optional) A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with a uint mask for the flux sampled at the provided wavelengths.
 no_offset (bool) – (Optional) Force \(\sigma^2_{v,o} =
0\) by masking regions with \(\sigma_{p,d} <
\epsilon_\sigma\); i.e., the value of this arguments selects
Option 1 (True) or Option 2 (False). See
SpectralResolution.GaussianKernelDifference()
.  min_sig_pix (float) – (Optional) Minimum value of the standard deviation in pixels allowed before assuming the kernel is a Delta function.
 variable_offset (bool) – (Optional) Flag to allow the offset applied to each spectrum (when the input contains more than one spectraum) to be tailored to each spectrum. Otherwise (variable_offset=False) the offset is forced to be the same for all spectra.
 log10 (bool) – (Optional) Flag that the spectrum has been binned logarithmically (base 10) in wavelength
 new_log10 (bool) – (Optional) Flag that the coordinates of the new spectral resolution are spectrum as been binned logarithmically (base 10) in wavelength.
Returns: Five objects are returned:
 A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the resolutionmatched flux sampled at the input wavelengths.
 A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the spectral resolution, \(R\), of the resolutionmatched spectra at the provided wavelengths.
 A 1D vector with any constant offset in resolution in
km/s between the targetted value and the result. See
SpectralResolution.GaussianKernelDifference()
.  A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with a uint mask for the resolutionmatched flux sampled at the input wavelengths. This is returned regardless of whether an input mask was provided. Any pixel that had a resolution that was lower than the target resolution (up to some tolerance defined by min_sig_pix) is returned as masked.
 A 1D or 2D (\(N_{\rm spec}\times N_{\rm pix}\)) array with the inverse variance of the resolutionmatched flux sampled at the input wavelengths. If ivar is not provided, a ‘None’ returned as the last element
Return type: numpy.ndarray
Raises: ValueError
– Raised if: the input wave array is 2D and the sres array is not; a 1D wavelength array is allowed for a 2D sres array but not vice versa
 the number of spectral pixels in wave, flux, and sres is not the same
 the shape of the flux, mask (if provided), and ivar (if provided) are not the same
 the shape of the new_sres_wave and new_sres arrays are not the same and/or not 1D
Todo
 Add interp_ext != ‘extrapolate’ option?
 Better way to use warnings?