Repository: radis/radis (language: python)

Issue:
Add bridges to/from Specutils
### 🎯 Describe what you'd like
<!--- A clear and concise description of what you want to happen. -->

As mentioned in #35, PySpeckit and Specutils are two spectroscopy-related librairies, with post-processing and fitting abilities. Following a recent call in the Astropy community to gather spectroscopy-related efforts, the goal is to build bridges to/from these librairies. 

Radis could benefit from using the post-processing & fitting functions developed and maintained there ; and Radis-calculated spectra (HITRAN/HITEMP/ExoMol) could be used as PySpeckit models or to compare with experiments processed with Specutils. 

### 💡 Possible implementations
<!--- Not obligatory, but suggest an idea for implementing addition or change -->

[specutil](https://github.com/astropy/specutils) is an Astropy affiliated package for Spectral analysis, maintained by @nmearl @crawfordsm @keflavich @eteq.  The post-processing tools are based on the [Spectrum1D](https://specutils.readthedocs.io/en/latest/api/specutils.Spectrum1D.html#specutils.Spectrum1D) class 

[pyspeckit](https://github.com/pyspeckit/pyspeckit) is the Python Spectroscopic Toolkit, "designed to allow for analysis of spectroscopic data from a wide variety of astronomical instruments. ", maintained by @keflavich and @mirochaj .  It is based on the [Spectrum](https://pyspeckit.readthedocs.io/en/latest/spectrum.html#pyspeckit.spectrum.classes.Spectrum) class 

(note that unlike these Spectrum and Spectrum1D classes above, Radis's Spectrum is more like a "column of emitting/absorbing gas", with both absorption & emission feature, and the ability to be rescaled etc. 

The goal would be to build 2 functions : 

- [x] Spectrum.to_specutils() function
- [x] Spectrum.from_specutils() function
- Spectrum.to_pyspeckit() function   --> Not needed as Pyspeckit integrates with specutils (see comment below)
- Spectrum.from_pyspeckit() function --> idem

For instance, Spectrum.to_specutils could look like this : 
```python
def to_specutils(s):
    """Convert a ``radis`` :py:class:`~radis.spectrum.spectrum.Spectrum` object to ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D`"""
    import astropy.units as u
    from specutils.spectra import Spectrum1D
    w, I = s.get("radiance", wunit=s.get_waveunit(), Iunit=s.units["radiance"])
    return Spectrum1D(flux=I*u.Unit(s.units["radiance"]), spectral_axis=w*u.Unit(s.get_waveunit()), meta=s.conditions)

```

- [x] To keep RADIS package as small as possible, support should be optional, i.e. imports of specutils/pyspeckit should only happen within function calls. 





---
<sub>*Tip : you can vote for [New Features](https://github.com/radis/radis/issues/79) on https://feathub.com/radis/radis*</sub>



Discussion / hints:
I'd suggest just integrating with specutils; pyspeckit should be able to work with any specutils objects, so you don't need a separate layer for it.

Fix patch:
diff --git a/docs/spectrum/spectrum.rst b/docs/spectrum/spectrum.rst
index 959cdfbd7..b15398ad4 100644
--- a/docs/spectrum/spectrum.rst
+++ b/docs/spectrum/spectrum.rst
@@ -103,8 +103,14 @@ convenience function::
     s = Spectrum.from_array(w, T, 'transmittance_noslit',
                                wunit='nm', unit='') # adimensioned
 
+Dimensionned arrays can also be used directly ::
 
-Convenience functions have been added to handle the usual cases:
+    import astropy.units as u
+    w = np.linspace(200, 300) * u.nm
+    I = np.random.rand(len(w)) * u.mW/u.cm**2/u.sr/u.nm
+    s = Spectrum.from_array(w, I, 'radiance_noslit')
+
+Other convenience functions have been added to handle the usual cases:
 :func:`~radis.spectrum.models.calculated_spectrum`,
 :func:`~radis.spectrum.models.transmittance_spectrum` and
 :func:`~radis.spectrum.models.experimental_spectrum`::
@@ -116,6 +122,16 @@ Convenience functions have been added to handle the usual cases:
     s3 = experimental_spectrum(w, I, wunit='nm', Iunit='W/cm2/sr/nm')   # creates 'radiance'
 
 
+Initialize from Specutils
+-------------------------
+
+Use :py:meth:`~radis.spectrum.spectrum.Spectrum.from_specutils` to convert
+from a ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D` object ::
+
+    from radis import Spectrum
+    Spectrum.from_specutils(spectrum)
+
+
 Initialize from a text file
 ---------------------------
 
@@ -156,9 +172,21 @@ Spectrum conditions and get Spectrum that suits specific parameters
 .. minigallery:: radis.load_spec
 
 
-Want a quick spectrum
+Load from a HDF5 file
 ---------------------
 
+This is the fastest way to read a Spectrum object from disk. It keeps metadata
+and units, and you can also load only a part of a very large spectrum.
+Use :py:meth:`~radis.spectrum.spectrum.Spectrum.from_hdf5`  ::
+
+    Spectrum.from_hdf5("rad_hdf.h5", wmin=2100, wmax=2200, columns=['abscoeff', 'emisscoeff'])
+
+
+
+
+Calculate a test spectrum
+-------------------------
+
 You need a spectrum immediatly, to run some tests ? Use :py:func:`~radis.test.utils.test_spectrum` ::
 
     s = radis.test_spectrum()
@@ -446,6 +474,14 @@ automatically retrieves a spectrum from a database if you calculated it already,
 or calculates it and stores it if you didn't. Very useful for spectra that
 requires long computational times!
 
+Export to hdf5
+--------------
+
+This is the fastest way to dump a Spectrum object on disk (and also, it keeps metadata
+and therefore units !). Use :py:meth:`~radis.spectrum.spectrum.Spectrum.to_hdf5`  ::
+
+    s.to_hdf5('spec01.h5')
+
 
 Export to txt
 -------------
@@ -460,6 +496,28 @@ will export a given spectral arrays::
 
     s.savetxt('radiance_W_cm2_sr_um.csv', 'radiance_noslit', wunit='nm', Iunit='W/cm2/sr/µm')
 
+Export to Pandas
+----------------
+
+Use :py:meth:`~radis.spectrum.spectrum.Spectrum.to_pandas` ::
+
+    s.to_pandas()
+
+This will return a DataFrame with all spectral arrays as columns.
+
+
+Export to Specutils
+-------------------
+
+Use :py:meth:`~radis.spectrum.spectrum.Spectrum.to_specutils` to convert
+to a to ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D` object ::
+
+    s.to_specutils()
+
+
+.. minigallery:: radis.spectrum.spectrum.Spectrum.to_specutils
+
+
 
 How to modify a Spectrum object?
 ================================
@@ -605,13 +663,15 @@ For instance, remove a constant baseline in a given unit::
 
     s -= 0.1 * u.Unit('W/cm2/sr/nm')
 
-Here, se normalize a spectrum manually, assuming the spectrum units is in "mW/cm2/sr/nm" ::
+The :py:meth:`~radis.spectrum.spectrum.Spectrum.max` function returns a dimensionned
+value, therefore it can be used to normalize a spectrum directly :  ::
 
-    s /= s.max() * u.Unit("mW/cm2/sr/nm")
+    s /= s.max()
 
-Or below, we calibrate a Spectrum, assuming the spectrum units is in "count" ::
+Or below, we calibrate a Spectrum, assuming the spectrum units is in "count",
+and that our calibration show we have 94 :math:`mW/cm2/sr/nm` per count.  ::
 
-    s *= 100 * u.Unit("mW/cm2/sr/nm/count")
+    s *= 94 * u.Unit("mW/cm2/sr/nm") / u.Unit("count")
 
 
 .. _label_spectrum_offset_crop:
@@ -666,7 +726,7 @@ Or::
 If you want to create a new spectrum, don't forget to set ``inplace=False``
 for the first command that allows it. i.e ::
 
-    s2 = s.crop(4120, 4220, 'nm', inplace=False).apply_slit(3, 'nm').take('radiance')
+    s2 = s.crop(4120, 4220, 'nm', inplace=False).apply_slit(3, 'nm').offset(1.5, 'nm')
 
 
 .. _label_spectrum_remove_baseline:
@@ -684,6 +744,9 @@ Or remove a linear baseline with:
 - :py:func:`~radis.spectrum.operations.get_baseline`
 - :py:func:`~radis.spectrum.operations.sub_baseline`
 
+You could also use the functions available in :py:mod:`pyspecutils`,
+see :py:meth:`~radis.spectrum.spectrum.Spectrum.to_specutils`.
+
 
 Calculate transmittance from radiance with Kirchoff's law
 ---------------------------------------------------------
diff --git a/environment.yml b/environment.yml
index 501b617ed..de3bd871e 100644
--- a/environment.yml
+++ b/environment.yml
@@ -6,25 +6,25 @@ channels:
 - cantera
 dependencies:
 - python=3.8
-- numpy
-- scipy>=1.4.0
-- matplotlib
-- seaborn   # other matplotlib themes
-- cantera   # for chemical equilibrium computations
-- cython
-- pandas>=1.0.5
-- pytables  # for pandas to HDF5 export
-- plotly    # used in line-survey
 - astropy  # Unit aware calculations
-- plotly>=2.5.1  # for line survey HTML output
-- termcolor  # terminal colors
-- configparser
 - astroquery>=0.3.9  # to fetch HITRAN databases
-- pytest  # to run test suite
+- cantera   # for chemical equilibrium computations
+- configparser
+- cython
+- h5py   # load HDF5
 - joblib  # for parallel loading of SpecDatabase
+- matplotlib
+- numpy
 - numba  # just-in-time compiler
+- pandas>=1.0.5
+- plotly>=2.5.1  # for line survey HTML output
 - psutil # for getting user RAM
-- h5py   # load HDF5
+- pytables  # for pandas to HDF5 export
+- pytest  # to run test suite
+- scipy>=1.4.0
+- seaborn   # other matplotlib themes
+- termcolor  # terminal colors
+- specutils
 - pip
 - pip:
   - beautifulsoup4 # parse ExoMol website
diff --git a/examples/plot_specutils_processing.py b/examples/plot_specutils_processing.py
new file mode 100644
index 000000000..358b4a6ab
--- /dev/null
+++ b/examples/plot_specutils_processing.py
@@ -0,0 +1,72 @@
+# -*- coding: utf-8 -*-
+"""
+============================
+Post-process using Specutils
+============================
+
+Find peaks or uncertainties using the :py:mod:`specutils` library. A Radis Spectrum
+object can easily be converted to a ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D`
+using :py:meth:`~radis.spectrum.spectrum.Spectrum.to_specutils`.
+
+Below, we create a noisy spectrum based on a synthetic CO spectrum,
+we convert it to :py:mod:`specutils`, add uncertainties by targeting a
+noisy region, then determine the lines using :py:func:`~specutils.fitting.find_lines_threshold` :
+
+"""
+
+import astropy.units as u
+import numpy as np
+
+from radis import test_spectrum
+
+""" We create a synthetic CO spectrum"""
+
+s = (
+    test_spectrum(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() * 0.03, len(s))
+s_exp = s + noise
+
+s_exp.plot()
+
+#%% 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)
+
+
+#%% Note: we can also create a RADIS spectrum object from the Specutils Spectrum1D :
+
+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
diff --git a/radis/io/hitran.py b/radis/io/hitran.py
index 1d128f449..cb012c9be 100644
--- a/radis/io/hitran.py
+++ b/radis/io/hitran.py
@@ -1042,6 +1042,6 @@ def parse_global_quanta(df, mol, verbose=True):
 
 if __name__ == "__main__":
 
-    from radis.test.test_io import _run_testcases
+    from radis.test.io.test_hitran_cdsd import _run_testcases
 
     print("Testing HITRAN parsing: ", _run_testcases())
diff --git a/radis/lbl/calc.py b/radis/lbl/calc.py
index 605d3aea7..ebad7593e 100644
--- a/radis/lbl/calc.py
+++ b/radis/lbl/calc.py
@@ -61,6 +61,7 @@ def calc_spectrum(
     mode="cpu",
     export_lines=False,
     verbose=True,
+    return_factory=False,
     **kwargs
 ) -> Spectrum:
     r"""Calculate a :py:class:`~radis.spectrum.spectrum.Spectrum`.
@@ -244,9 +245,18 @@ def calc_spectrum(
         Note that ``mode='gpu'`` requires CUDA compatible hardware to execute.
         For more information on how to setup your system to run GPU-accelerated
         methods using CUDA and Cython, check `GPU Spectrum Calculation on RADIS <https://radis.readthedocs.io/en/latest/lbl/gpu.html>`__
+    return_factory: bool
+        if ``True``, return the :py:class:`~radis.lbl.factory.SpectrumFactory` that
+        computes the spectrum. Useful to access computational parameters, the line database,
+        or to start batch-computations from a first spectrum calculation. Ex::
+
+                s, sf = calc_spectrum(..., return_factory=True)
+                sf.df1  # see the lines calculated
+                sf.eq_spectrum(...)  #  new calculation without reloading the database
+
     **kwargs: other inputs forwarded to SpectrumFactory
         For instance: ``warnings``.
-        See :class:`~radis.lbl.factory.SpectrumFactory` documentation for more
+        See :py:class:`~radis.lbl.factory.SpectrumFactory` documentation for more
         details on input.
 
     Returns
@@ -415,6 +425,7 @@ def _check_molecules_are_consistent(
 
     # ... Initialize and fill the master-dictionary
     molecule_dict = {}
+    factory_dict = {}
 
     for molecule in molecule_reference_set:
         molecule_dict[molecule] = {}
@@ -502,8 +513,12 @@ def _check_molecules_are_consistent(
             verbose=verbose,
             mode=mode,
             export_lines=export_lines,
+            return_factory=return_factory,
             **kwargs_molecule
         )
+        if return_factory:
+            factory_dict[molecule] = generated_spectrum[1]
+            generated_spectrum = generated_spectrum[0]  # the spectrum
         s_list.append(generated_spectrum)
 
         if condition_multiple_wstep:
@@ -519,7 +534,10 @@ def _check_molecules_are_consistent(
     if save_to:
         s.store(path=save_to, if_exists_then="replace", verbose=verbose)
 
-    return s
+    if return_factory:
+        return s, factory_dict
+    else:
+        return s
 
 
 def _calc_spectrum_one_molecule(
@@ -548,8 +566,9 @@ def _calc_spectrum_one_molecule(
     verbose,
     mode,
     export_lines,
+    return_factory=False,
     **kwargs
-):
+) -> Spectrum:
     """See :py:func:`~radis.lbl.calc.calc_spectrum`"""
 
     # Check inputs
@@ -747,7 +766,10 @@ def _is_at_equilibrium():
             name=name,
         )
 
-    return s
+    if return_factory:
+        return s, sf
+    else:
+        return s
 
 
 # --------------------------
diff --git a/radis/spectrum/spectrum.py b/radis/spectrum/spectrum.py
index a48ff6553..676fb3996 100644
--- a/radis/spectrum/spectrum.py
+++ b/radis/spectrum/spectrum.py
@@ -519,7 +519,18 @@ def __init__(
     # %% Constructors
 
     @classmethod
-    def from_array(self, w, I, quantity, wunit, unit, waveunit=None, *args, **kwargs):
+    def from_array(
+        self,
+        w,
+        I,
+        quantity,
+        wunit=None,
+        Iunit=None,
+        waveunit=None,
+        unit=None,
+        *args,
+        **kwargs,
+    ):
         """Construct Spectrum from 2 arrays.
 
         Parameters
@@ -531,8 +542,10 @@ def from_array(self, w, I, quantity, wunit, unit, waveunit=None, *args, **kwargs
         wunit: ``'nm'``, ``'cm-1'``, ``'nm_vac'``
             unit of waverange:         wavelength in air (``'nm'``), wavenumber
             (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``).
-        unit: str
-            spectral quantity unit (arbitrary). Ex: 'mW/cm2/sr/nm' for radiance_noslit
+            If ``None``, then ``w`` must be a dimensionned array.
+        Iunit: str
+            spectral quantity unit (arbitrary). Ex: ``'mW/cm2/sr/nm'`` for radiance_noslit
+            If ``None``, then ``I`` must be a dimensionned array.
         *args, **kwargs
             see :class:`~radis.spectrum.spectrum.Spectrum` doc
 
@@ -561,19 +574,24 @@ def from_array(self, w, I, quantity, wunit, unit, waveunit=None, *args, **kwargs
 
         Returns
         -------
-
-        s: Spectrum
+        Spectrum
             creates a :class:`~radis.spectrum.spectrum.Spectrum` object
 
         Examples
         --------
-
         Create a spectrum::
 
             from radis import Spectrum
             s = Spectrum.from_array(w, I, 'radiance_noslit',
                                    wunit='nm', unit='mW/cm2/sr/nm')
 
+        Dimensionned arrays can also be used directly ::
+
+            import astropy.units as u
+            w = np.linspace(200, 300) * u.nm
+            I = np.random.rand(len(w)) * u.mW/u.cm**2/u.sr/u.nm
+            s = Spectrum.from_array(w, I, 'radiance_noslit')
+
         To create a spectrum with absorption and emission components
         (e.g: ``radiance_noslit`` and ``transmittance_noslit``, or ``emisscoeff``
         and ``abscoeff``) call the :class:`~radis.spectrum.spectrum.Spectrum`
@@ -604,9 +622,56 @@ class directly. Ex::
                 DeprecationWarning,
             )
             wunit = waveunit
+        if unit is not None:
+            warn(
+                "`unit=` parameter in from_array is now named `Iunit=`",
+                DeprecationWarning,
+            )
+            Iunit = unit
+
+        # Check if dimensionned arrays
+        try:
+            I.unit
+        except AttributeError:
+            # dimensionless
+            if Iunit is None:
+                raise ValueError(
+                    "I must be a dimensionned array, or ``Iunit=`` must be given."
+                )
+        else:
+            if Iunit is not None:
+                raise ValueError(
+                    f"I is a dimensionned array (in {Iunit.unit.to_string()}), therefore ``Iunit`` cannot be given too (got {Iunit}). Set ``Iunit=None``"
+                )
+            Iunit = I.unit.to_string()
+            I = I.value
+        try:
+            w.unit
+        except AttributeError:
+            # dimensionless
+            if wunit is None:
+                raise ValueError(
+                    "``w`` must be a dimensionned array, or ``Iunit=`` must be given."
+                )
+        else:
+            if wunit is not None:
+                raise ValueError(
+                    f"``w`` is a dimensionned array (in {wunit.unit.to_string()}), therefore ``wunit`` cannot be given too (got {wunit}). Set ``wunit=None``"
+                )
+            wunit = w.unit.to_string()
+            if wunit not in WAVELEN_UNITS + WAVENUM_UNITS:
+                # Convert to something we know, for instance 'nm' :
+                if w.unit.is_equivalent("nm"):
+                    w = w.to("nm")
+                    wunit = "nm"
+                elif w.unit.is_equivalent("1 / cm"):
+                    w = w.to("1 / cm")
+                    wunit = "cm-1"
+                # else, an error will be raised anyway on Spectrum creation.
+            w = w.value
 
         quantities = {quantity: (w, I)}
-        units = {quantity: unit}
+        units = {quantity: Iunit}
 
         # Update Spectrum conditions
         conditions = kwargs.pop("conditions", {})
@@ -767,7 +832,12 @@ def from_hdf5(
         Examples
         --------
         ::
+
             Spectrum.from_hdf5("rad_hdf.h5", wmin=2100, wmax=2200, columns=['abscoeff', 'emisscoeff'])
+
+        See Also
+        --------
+        :py:func:`~radis.spectrum.spectrum.Spectrum.to_hdf5`
         """
         from radis.io.spec_hdf import hdf2spec
 
@@ -786,6 +856,64 @@ def from_spec(self, file):
 
         return load_spec(file)
 
+    @classmethod
+    def from_specutils(self, spectrum, var="radiance"):
+        """Convert a ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D`
+        to a ``radis`` :py:class:`~radis.spectrum.spectrum.Spectrum` object.
+
+        Examples
+        --------
+
+        Taken from the Specutils website (https://specutils.readthedocs.io/en/stable/#getting-started-with-specutils)
+        ::
+
+            from astropy.io import fits
+            from astropy import units as u
+            from specutils import Spectrum1D
+
+            f = fits.open('https://data.sdss.org/sas/dr16/sdss/spectro/redux/26/spectra/1323/spec-1323-52797-0012.fits')
+            # The spectrum is in the second HDU of this file.
+            specdata = f[1].data
+
+            lamb = 10**specdata['loglam'] * u.AA
+            flux = specdata['flux'] * 10**-17 * u.Unit('erg cm-2 s-1 AA-1')
+            spec = Spectrum1D(spectral_axis=lamb, flux=flux)
+
+
+            from radis import Spectrum
+            s = Spectrum.from_specutils(spec)
+            s.plot(wunit='nm')
+
+
+        .. minigallery:: radis.spectrum.spectrum.Spectrum.from_specutils
+
+        See Also
+        --------
+        :py:func:`~radis.spectrum.spectrum.Spectrum.to_specutils`
+        """
+        try:
+            import specutils
+        except ModuleNotFoundError as err:
+            raise ModuleNotFoundError(
+                "Specutils is required to use this function."
+                "Install it with:\n"
+                "   conda install -c conda-forge specutils"
+                "\nor\n"
+                "   pip install specutils"
+            ) from err
+
+        assert isinstance(spectrum, specutils.Spectrum1D)
+
+        # Convert to radis Spectrum
+        # note : in Specutils wavelengths are given as wavelength in vacuum (1e7/X conversion
+        # to wavenumbers), so here we give the frequency unit.
+        return Spectrum.from_array(
+            spectrum.wavelength.to("1 / cm"),
+            spectrum.flux,
+            quantity=var,
+            conditions=spectrum.meta,
+        )
+
     # Public functions
     # %% ======================================================================
     # ----------------
@@ -2491,7 +2619,7 @@ def apply_slit(
                 w_nm, wslit0_nm, slit_dispersion, slit_dispersion_threshold
             )
         else:
-            slice_windows = [np.ones_like(w, dtype=np.bool)]
+            slice_windows = [np.ones_like(w, dtype=bool)]
 
         # Create dictionary to store convolved
         I_conv_slices = {}
@@ -3105,7 +3233,18 @@ def to_json(self, *args, **kwargs):
         return self.store(compress=False, *args, **kwargs)
 
     def to_hdf5(self, file, engine="pytables"):
-        """Stores the Spectrum under HDF5 format. Uses :py:func:`~radis.io.spec_hdf.spec2hdf`"""
+        """Stores the Spectrum under HDF5 format. Uses :py:func:`~radis.io.spec_hdf.spec2hdf`
+
+        Examples
+        --------
+        ::
+
+            s.to_hdf5('spec01.h5')
+
+        See Also
+        --------
+        :py:func:`~radis.spectrum.spectrum.Spectrum.from_hdf5`
+        """
         from radis.io.spec_hdf import spec2hdf
 
         return spec2hdf(self, file, engine=engine)
@@ -3133,6 +3272,91 @@ def to_pandas(s, copy=True):
 
         return df
 
+    def to_specutils(self, var=None, wunit="default", Iunit="default"):
+        """Convert a ``radis`` :py:class:`~radis.spectrum.spectrum.Spectrum`
+        object to ``specutils`` :py:class:`specutils.spectra.spectrum1d.Spectrum1D`
+
+        Parameters
+        ----------
+        var: 'radiance', 'radiance_noslit', etc.
+            which spectral array to convert. If ``None`` and only one spectral
+            arry is defined, use it
+        wunit: ``'nm'``, ``'cm'``, ``'nm_vac'``.
+            wavespace unit: wavelength in air (``'nm'``), wavenumber
+            (``'cm-1'``), or wavelength in vacuum (``'nm_vac'``).
+            if ``"default"``, default unit for waveunit is used. See
+            :py:meth:`~radis.spectrum.spectrum.Spectrum.get_waveunit`.
+        Iunit: unit for variable ``var``
+            if ``"default"``, default unit for quantity `var` is used. See the
+            :py:attr:`~radis.spectrum.spectrum.Spectrum.units` attribute.
+            For ``var="radiance"``, one can use per wavelength (~ 'W/m2/sr/nm')
+            or per wavenumber (~ 'W/m2/sr/cm-1') units
+
+        .. note::
+            `nan`, that may have been added on the wings of the spectra if a
+            slit has been applied, are removed using ``trim_nan`` parameter of
+            :py:meth:`~radis.spectrum.spectrum.Spectrum.get` . The waverange
+            in the 1DSpectrum object may therefore be cropped.
+
+        Examples
+        --------
+        ::
+
+            spectrum = s.to_specutils()
+
+        Add uncertainties by reading a spectrum noise region with
+        :py:class:`~specutils.SpectralRegion` and
+        :py:func:`~specutils.manipulation.noise_region_uncertainty` ::
+
+            from specutils import SpectralRegion
+            from specutils.manipulation import noise_region_uncertainty
+            noise_region = SpectralRegion(2012 / u.cm, 2009 / u.cm)
+            spectrum = noise_region_uncertainty(spectrum, noise_region)
+
+        Find lines out of the noise using :py:func:`~specutils.fitting.find_lines_threshold` ::
+
+            from specutils.fitting import find_lines_threshold
+            lines = find_lines_threshold(spectrum)
+
+        .. minigallery:: radis.spectrum.spectrum.Spectrum.to_specutils
+            :add-heading:
+
+        See Also
+        --------
+        :py:func:`~radis.spectrum.spectrum.Spectrum.from_specutils`
+
+        """
+        import astropy.units as u
+
+        try:
+            from specutils.spectra import Spectrum1D
+        except ModuleNotFoundError as err:
+            raise ModuleNotFoundError(
+                "Specutils is required to use this function."
+                "Install it with:\n"
+                "   conda install -c conda-forge specutils"
+                "\nor\n"
+                "   pip install specutils"
+            ) from err
+
+        if var is None:
+            var = self._get_unique_var(operation_name="to_specutils")
+        if wunit == "default":
+            wunit = self.get_waveunit()
+        if Iunit == "default":
+            Iunit = self.units[var]
+
+        if wunit in WAVELEN_UNITS:
+            wunit = "nm_vac"  # AFAIK, specutil's Spectrum1D will only handle wavelengths as seen in vacuum
+            # so below we request wavelengths in vac :
+
+        w, I = self.get(var, wunit=wunit, Iunit=Iunit, copy=True, trim_nan=True)
+        return Spectrum1D(
+            flux=I * u.Unit(Iunit),
+            spectral_axis=w * u.Unit(wunit),
+            meta=self.get_conditions(),
+        )
+
     def resample(
         self,
         w_new,
@@ -4250,7 +4474,7 @@ def __add__(self, other):
         elif isinstance(other, np.ndarray):
             from radis.spectrum.operations import add_array
 
-            return add_array(self, other, inplace=True)
+            return add_array(self, other, inplace=False)
         elif isinstance(other, Spectrum):
             from radis.spectrum.operations import add_spectra
 
diff --git a/radis/tools/line_survey.py b/radis/tools/line_survey.py
index c6bdbb3ea..5c0ef99fc 100644
--- a/radis/tools/line_survey.py
+++ b/radis/tools/line_survey.py
@@ -466,7 +466,7 @@ def get_label_nist(row, attrs):
     def get_label_all(row):
         """ print all lines details """
         # label = row.__repr__()
-        label = "\n".join(
+        label = "<br>".join(
             [f"{k}: {v}" for k, v in row.items() if k not in ["wav", "shiftwav"]]
         )
         return label


Test patch:
diff --git a/radis/test/files/CO2_measured_spectrum_4-5um.spec b/radis/test/files/CO2_measured_spectrum_4-5um.spec
index eb36394c6..1dde49e69 100644
Binary files a/radis/test/files/CO2_measured_spectrum_4-5um.spec and b/radis/test/files/CO2_measured_spectrum_4-5um.spec differ
diff --git a/radis/test/files/CO_Tgas1500K_mole_fraction0.01.spec b/radis/test/files/CO_Tgas1500K_mole_fraction0.01.spec
index 6b9033c97..34b8dd3c0 100644
Binary files a/radis/test/files/CO_Tgas1500K_mole_fraction0.01.spec and b/radis/test/files/CO_Tgas1500K_mole_fraction0.01.spec differ
diff --git a/radis/test/files/CO_Tgas1500K_mole_fraction0.5.spec b/radis/test/files/CO_Tgas1500K_mole_fraction0.5.spec
index 3483b8504..49748b0cc 100644
Binary files a/radis/test/files/CO_Tgas1500K_mole_fraction0.5.spec and b/radis/test/files/CO_Tgas1500K_mole_fraction0.5.spec differ
diff --git a/radis/test/files/N2C_specair_380nm.spec b/radis/test/files/N2C_specair_380nm.spec
index fdfdeb2e4..a78380b64 100644
--- a/radis/test/files/N2C_specair_380nm.spec
+++ b/radis/test/files/N2C_specair_380nm.spec
@@ -7525,7 +7525,6 @@
             ]
         }
     },
-    "_q_conv": {},
     "units": {
         "radiance_noslit": "mW/cm2/sr/um",
         "transmittance_noslit": "1"
@@ -7561,7 +7560,6 @@
         "wo2": 0.0,
         "nprt1": 2,
         "nms": 9,
-        "medium": "air",
         "narray": 2500,
         "wh2o": 0.0,
         "bandoriginsaver": "0",
@@ -7590,6 +7588,11 @@
         "wavelength_min": "nm"
     },
     "populations": {},
+    "references": {
+        "10.1016/j.jqsrt.2018.09.027": [
+            "post-processing"
+        ]
+    },
     "lines": null,
     "name": "N2C_specair_380nm",
     "_slit": {},