Smoothing models to resolution

[mileslucas]

I have an idea for implementing a method for SourceSpectrum objects to smooth them to a given resolution. I have some working code that convolves gaussian kernels to the spectrum in logspace. I was curious if there is merit in putting the effort into creating a nice pull request with documentation and tests.

[pllim]

@mileslucas , thank you for your offer! I think I had seen feature request like this over at AstroBetter. But before you put in all that effort, how about first, a Jupyter notebook gist of what you already have, and we'll go from there?

[mileslucas]

Functions:

import pysynphot as S
from astropy.convolution import convolve_fft
from astropy.convolution import Gaussian1DKernel

def smear(sp, R, w_sample=1):
    '''
    Smears a model spectrum with a gaussian kernel to the given resolution, R.

    Parameters
    -----------

    sp: SourceSpectrum
        Pysynphot object that we willsmear

    R: int
        The resolution (dL/L) to smear to

    w_sample: int
        Oversampling factor for smoothing

    Returns
    -----------

    sp: PySynphot Source Spectrum
        The smeared spectrum
    '''

    # Save original wavelength grid and units
    w_grid = sp.wave
    w_units = sp.waveunits
    f_units = sp.fluxunits
    sp_name = sp.name

    # Generate logarithmic wavelength grid for smoothing
    w_logmin = np.log10(np.nanmin(w_grid))
    w_logmax = np.log10(np.nanmax(w_grid))
    n_w = np.size(w_grid)*w_sample
    w_log = np.logspace(w_logmin, w_logmax, num=n_w)

    # Find stddev of Gaussian kernel for smoothing
    R_grid = (w_log[1:-1]+w_log[0:-2])/(w_log[1:-1]-w_log[0:-2])/2
    sigma = np.median(R_grid)/R
    if sigma < 1:
        sigma = 1

    # Interpolate on logarithmic grid
    f_log = np.interp(w_log, w_grid, sp.flux)

    # Smooth convolving with Gaussian kernel
    gauss = Gaussian1DKernel(stddev=sigma)
    f_conv = convolve_fft(f_log, gauss)

    # Interpolate back on original wavelength grid
    f_sm = np.interp(w_grid, w_log, f_conv)

    # Write smoothed spectrum back into Spectrum object
    return S.ArraySpectrum(w_grid, f_sm, waveunits=w_units,
                            fluxunits=f_units, name=sp_name)

Some examples smoothing HiRes ACES Phoenix models:

Smoothed to R=2000

Profile Comparison at Hydrogen Pa-7

My ideal usage would be implementing the smoothing as a member function of SourceSpectrum, so example ideal usage would be

sp_smoothed = sp.smooth(R=2000, w_sample=100)

Note:
This code was originally made by my mentor. I am not as comfortable with the mathematics of the gaussian convolving. Perhaps there is a better method of smoothing?

[pllim]

@mileslucas , thank you but after thinking about this for a bit, I think this is better off as a separate tutorial-type documentation rather than being an actual function in this package. As you mentioned, different people might have their own favorite way to convolve a spectrum, hence being a tutorial would provide more flexibility.

Furthermore, since you are using Astropy for the convolution kernel, perhaps this is even better using our refactored synphot package, which is an Astropy-affiliated package and uses Astropy models. If you go that route, then such a tutorial would naturally be part of https://github.com/astropy/astropy-tutorials .