Post-process using Specutils

Find peaks or uncertainties using the specutils library. A Radis Spectrum object can easily be converted to a specutils specutils.spectra.spectrum.Spectrum using to_specutils().

Below, we create a noisy spectrum based on a synthetic CO spectrum, we convert it to specutils, add uncertainties by targeting a noisy region, then determine the lines using find_lines_threshold() :

import astropy.units as u
import numpy as np

from radis import spectrum_test

""" We create a synthetic CO spectrum"""

s = (
    spectrum_test(molecule="CO", wavenum_min=2000, wavenum_max=2030)
    .apply_slit(1.5, "nm")
    .take("radiance")
)
s.trim()  # removes nans created by the slit convolution boundary effects
noise = np.random.normal(0.0, s.max().value * 0.03, len(s))
s_exp = s + noise

s_exp.plot()

--------------------------------------------------------------------------------
CO - HITRAN - Downloading database
--------------------------------------------------------------------------------

Download:
- All files already downloaded.

Caching to HDF5/H5 format:
- All files already cached.
0.02s - Loaded database
Calculating Equilibrium Spectrum
Physical Conditions
----------------------------------------
   Tgas                 700 K
   isotope              1,2,3
   medium               air
   mole_fraction        0.1
   path_length          1 cm
   pressure             1.01325 bar
   self_absorption      True
   species              CO
   state                X
   wavenum_max          2030.0000 cm-1
   wavenum_min          2000.0000 cm-1
Computation Parameters
----------------------------------------
   Tref                 296 K
   add_at_used          numpy
   broadening_method    voigt_poly
   cutoff               1e-27 cm-1/(#.cm-2)
   dbformat             hitran
   dbpath               /home/docs/.radisdb/hitran/CO.hdf5
   diluent              air
   folding_thresh       1e-06
   include_neighbouring_lines  True
   isatom               False
   isneutral            None
   lbfunc               None
   memory_mapping_engine  auto
   neighbour_lines      0 cm-1
   optimization         simple
   parsum_mode          full summation
   pfsource             default
   potential_lowering   None
   pseudo_continuum_threshold  0
   sparse_ldm           True
   truncation           50 cm-1
   waveunit             cm-1
   wstep                0.01 cm-1
   zero_padding         3001
----------------------------------------
0.02s - Spectrum calculated

Determine the noise level by selecting a noisy region from the graph above :

spectrum = s_exp.to_specutils()

from specutils import SpectralRegion
from specutils.manipulation import noise_region_uncertainty

noise_region = SpectralRegion(2010.5 / u.cm, 2009.5 / u.cm)
spectrum = noise_region_uncertainty(spectrum, noise_region)

Find lines :

from specutils.fitting import find_lines_threshold

lines = find_lines_threshold(spectrum, noise_factor=2)

print(lines)


s_exp.plot(lw=2, show_ruler=True)
import matplotlib.pyplot as plt

for line in lines.to_pandas().line_center.values:
    plt.axvline(line, color="r", zorder=-1)
s.plot(nfig="same")
plt.axvspan(noise_region.lower.value, noise_region.upper.value, color="b", alpha=0.1)

/home/docs/checkouts/readthedocs.org/user_builds/radis/envs/latest/lib/python3.12/site-packages/specutils/analysis/flux.py:289: AstropyUserWarning: Spectrum is not below the threshold signal-to-noise 0.01. This may indicate you have not continuum subtracted this spectrum (or that you have but it has high SNR features).

If you want to suppress this warning either type 'specutils.conf.do_continuum_function_check = False' or see http://docs.astropy.org/en/stable/config/#adding-new-configuration-items for other ways to configure the warning.
  warnings.warn(message, AstropyUserWarning)
   line_center     line_type  line_center_index
      1 / cm
------------------ ---------- -----------------
2000.6899999999994   emission                 8
2000.7699999999993   emission                16
2000.9399999999991   emission                33
 2001.119999999999   emission                51
2001.7499999999984   emission               114
2001.7799999999984   emission               117
2001.8699999999983   emission               126
2002.0299999999982   emission               142
2002.0699999999981   emission               146
               ...        ...               ...
2015.1899999999862 absorption              1458
2015.6199999999858 absorption              1501
2017.1299999999844 absorption              1652
2018.9099999999828 absorption              1830
2019.2799999999825 absorption              1867
 2019.809999999982 absorption              1920
2024.3799999999778 absorption              2377
2025.9499999999764 absorption              2534
 2026.399999999976 absorption              2579
 2028.649999999974 absorption              2804
Length = 251 rows
/home/docs/checkouts/readthedocs.org/user_builds/radis/checkouts/latest/radis/tools/plot_tools.py:615: UserWarning: Couldn't add Ruler tool (still an experimental feature in RADIS : please report the error !)
  warn(

<matplotlib.patches.Rectangle object at 0x7dd6bdeea600>

Note: we can also create a RADIS spectrum object from Specutils specutils.spectra.spectrum.Spectrum :

from radis import Spectrum

s2 = Spectrum.from_specutils(spectrum)

s2.plot(Iunit="mW/cm2/sr/nm", wunit="nm")
s_exp.plot(Iunit="mW/cm2/sr/nm", wunit="nm", nfig="same")
assert s_exp == s2