Open skendrew opened 3 years ago
This is an interesting use case! Currently, specutils doesn't handle using individual points to define a continuum particularly well, but it is possible with a little finagling:
from specutils import Spectrum1D
import numpy as np
import astropy.units as u
from specutils.fitting import fit_continuum
import matplotlib.pyplot as plt
spec = Spectrum1D(flux=np.random.sample(100) * u.Jy, spectral_axis=np.linspace(0, 50, 100) * u.AA)
# The points on the spectral axis corresponding to the continuum
cont_values = [7, 16, 24, 29.5, 34.5, 35.5, 36.5, 43, 45.5] * u.AA
# Create an array of indices for the values in `spec` that are closest to the
# continuum values. This just ensures that the exact values are used.
cont_idx = np.array([np.argmin(np.abs(spec.spectral_axis - x)) for x in cont_values])
# Create a "continuum spectrum" with the data found above, which represents
# only the continuum points we care about
cont_spec = Spectrum1D(flux=spec.flux[cont_idx], spectral_axis=spec.spectral_axis[cont_idx])
# Get the continuum model
cont_mod = fit_continuum(cont_spec)
I think we should certainly look into incorporating this feature.
Thanks Nick - yes this is the right idea I think. Here's a couple of plots to illustrate further the specific use case. The spectrum of this target looks like it's noisy, but every "bump" is actually a weak emission line (usually even a blend of several). It's hard to identify "regions" of continuum, and the automated fitter interprets the lines as noise. To get a good continuum fit we have to be able to pick out individual points (and you can see what the resulting spectrum should look like in the 2nd plot).
Some of our JWST mid-IR spectra with low resolution have dense forests of emission lines (H recombination lines). The automated continuum fit tends to interpret as noise what is actually all emission lines, and there is often only a single point between the lines. To get a good continuum fit we really need to pass individual points rather than regions. Would be good if that were possible but I couldn't see it in specutils.