Open mileslucas opened 6 years ago
@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?
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?
@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 .
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.