Merge pull request #1014 from ibusko/specviz-load-and-combine

Specviz: load a SpectrumList and combine everything

CHANGES.rst

@@ -11,6 +11,8 @@ New Features
 - The line analysis plugin now includes logic to account for the background
   continuum. [#1060]
 
+- Specviz can load a ``SpectrumList`` and combine all its elements into a single spectrum. [#1014]
+
 Cubeviz
 ^^^^^^^
 

docs/specviz/import_data.rst

@@ -46,6 +46,11 @@ Below is an example of importing the :class:`~jdaviz.configs.specviz.helper.Spec
     >>> specviz = Specviz()  # doctest: +SKIP
     >>> specviz.load_spectrum(spec1d)  # doctest: +SKIP
 
+You can also pass the path to a file that `~specutils.Spectrum1D` understands directly to the
+:meth:`jdaviz.configs.specviz.helper.Specviz.load_spectrum` method::
+
+    >>> specviz.load_spectrum("path/to/data/file") #doctest: +SKIP
+
 This method works well for data files that ``specutils`` understands.  However, if you are using your own data file or in-memory data, you can instead create a :class:`~specutils.Spectrum1D` object directly. In this example that is done using randomly generated data, and then that :class:`~specutils.Spectrum1D` object is loaded into the Specviz application::
 
     >>> from jdaviz import Specviz
@@ -60,3 +65,35 @@ This method works well for data files that ``specutils`` understands.  However,
 
 For more information about using the Specutils package, please see the
 `Specutils documentation <https://specutils.readthedocs.io>`_.
+
+Loading multiple spectra via the API
+------------------------------------
+In addition to loading single spectra as above, in some cases it may be useful
+to load multiple related spectra at once into the Jdaviz application. The most common
+such case is when you have spectra of the same object covering multiple wavelength
+ranges and want to look at and analyze the entire range of spectral coverage
+simultaneously. The :meth:`jdaviz.configs.specviz.helper.Specviz.load_spectrum` accepts
+a `~specutils.SpectrumList` object, in which case it will both load the
+individual `~specutils.Spectrum1D` objects in the list and additionally attempt
+to stitch together the spectra into a single data object so that
+they can be manipulated and analyzed in the application as a single entity::
+
+    >>> from specutils import SpectrumList
+    >>> spec_list = SpectrumList([spec1d_1, spec1d_2]) #doctest: +SKIP
+    >>> specviz.load_spectrum(spec_list) #doctest: +SKIP
+
+In the screenshot below, the combined spectrum is plotted in gray, and one of
+the single component spectra are also selected and plotted in red. Note that the
+"stitching" algorithm to combine the spectra is a simple concatenation of data,
+so in areas where the wavelength ranges of component spectra overlap you may see
+the line plot jumping between points of the two spectra, as at the beginning and
+end of the red region in the screenshot below:
+
+.. image:: img/spectrumlist_combined.png
+
+This functionality is also available in limited instances by providing a directory path
+to the :meth:`jdaviz.configs.specviz.helper.Specviz.load_spectrum` method. Note
+that the ``read`` method of :class:`~specutils.SpectrumList` is only set up to handle
+directory input in limited cases, for example JWST MIRI MRS data, and will throw an error
+in other cases. In cases that it does work, only files in the directory level specified
+will be read, with no recursion into deeper folders.

jdaviz/configs/specviz/plugins/parsers.py

@@ -2,7 +2,11 @@
 import pathlib
 import uuid
 
+import numpy as np
+
 from astropy.io.registry import IORegistryError
+from astropy.nddata import StdDevUncertainty
+
 from specutils import Spectrum1D, SpectrumList, SpectrumCollection
 
 from jdaviz.core.registries import data_parser_registry
@@ -38,8 +42,11 @@ def specviz_spectrum1d_parser(app, data, data_label=None, format=None, show_in_v
     elif isinstance(data, Spectrum1D):
         data = [data]
         data_label = [data_label]
+    # No special processing is needed in this case, but we include it for completeness
     elif isinstance(data, SpectrumList):
         pass
+    elif isinstance(data, list):
+        data = SpectrumList.read(data, format=format)
     else:
         path = pathlib.Path(data)
 
@@ -50,6 +57,11 @@ def specviz_spectrum1d_parser(app, data, data_label=None, format=None, show_in_v
             except IORegistryError:
                 # Multi-extension files may throw a registry error
                 data = SpectrumList.read(str(path), format=format)
+        elif path.is_dir():
+            data = SpectrumList.read(str(path), format=format)
+            if data == []:
+                raise ValueError(f"`specutils.SpectrumList.read('{str(path)}')` "
+                                 "returned an empty list")
         else:
             raise FileNotFoundError("No such file: " + str(path))
 
@@ -70,14 +82,58 @@ def specviz_spectrum1d_parser(app, data, data_label=None, format=None, show_in_v
         spec_key = list(current_spec.keys())[0]
         current_unit = current_spec[spec_key].spectral_axis.unit
     with app.data_collection.delay_link_manager_update():
-        for i in range(len(data)):
-            spec = data[i]
+
+        # these are used to build a combined spectrum with all
+        # input spectra included (taken from https://github.com/spacetelescope/
+        # dat_pyinthesky/blob/main/jdat_notebooks/MRS_Mstar_analysis/
+        # JWST_Mstar_dataAnalysis_analysis.ipynb)
+        wlallorig = []
+        fnuallorig = []
+        dfnuallorig = []
+
+        for i, spec in enumerate(data):
+
+            wave_units = spec.spectral_axis.unit
+            flux_units = spec.flux.unit
+
             if current_unit is not None and spec.spectral_axis.unit != current_unit:
                 spec = Spectrum1D(flux=spec.flux,
                                   spectral_axis=spec.spectral_axis.to(current_unit))
 
             app.add_data(spec, data_label[i])
 
-            # Only auto-show the first spectrum in a list
-            if i == 0 and show_in_viewer:
-                app.add_data_to_viewer("spectrum-viewer", data_label[i])
+            # handle display, with the SpectrumList special case in mind.
+            if show_in_viewer:
+                if isinstance(data, SpectrumList):
+
+                    # add spectrum to combined result
+                    for wlind in range(len(spec.spectral_axis)):
+                        wlallorig.append(spec.spectral_axis[wlind].value)
+                        fnuallorig.append(spec.flux[wlind].value)
+                        dfnuallorig.append(spec.uncertainty[wlind].array)
+
+                elif i == 0:
+                    app.add_data_to_viewer("spectrum-viewer", data_label[i])
+
+        # reset display ranges, or build combined spectrum, when input is a SpectrumList instance
+        if isinstance(data, SpectrumList):
+
+            # build combined spectrum
+            wlallarr = np.array(wlallorig)
+            fnuallarr = np.array(fnuallorig)
+            dfnuallarr = np.array(dfnuallorig)
+            srtind = np.argsort(wlallarr)
+            wlall = wlallarr[srtind]
+            fnuall = fnuallarr[srtind]
+            fnuallerr = dfnuallarr[srtind]
+
+            # units are not being handled properly yet.
+            unc = StdDevUncertainty(fnuallerr * flux_units)
+            spec = Spectrum1D(flux=fnuall * flux_units, spectral_axis=wlall * wave_units,
+                              uncertainty=unc)
+
+            # needs perhaps a better way to label the combined spectrum
+            label = "Combined " + data_label[0]
+            app.add_data(spec, label)
+            if show_in_viewer:
+                app.add_data_to_viewer("spectrum-viewer", label)

jdaviz/configs/specviz/tests/test_helper.py

@@ -1,13 +1,16 @@
+from zipfile import ZipFile
+
 import pytest
 from astropy import units as u
 from astropy.tests.helper import assert_quantity_allclose
 from glue.core.roi import XRangeROI
 from glue.core.edit_subset_mode import OrMode
 from specutils import Spectrum1D, SpectrumList, SpectrumCollection
+from astropy.utils.data import download_file
 
 from jdaviz.app import Application
 from jdaviz.configs.specviz.plugins.unit_conversion import unit_conversion as uc
-from jdaviz.configs.specviz.helper import Specviz
+from jdaviz import Specviz
 
 
 class TestSpecvizHelper:
@@ -31,14 +34,14 @@ def test_load_spectrum1d(self):
 
     def test_load_spectrum_list_no_labels(self):
         self.spec_app.load_spectrum(self.spec_list)
-        assert len(self.spec_app.app.data_collection) == 4
+        assert len(self.spec_app.app.data_collection) == 5
         for i in (1, 2, 3):
             assert "specviz_data" in self.spec_app.app.data_collection[i].label
 
     def test_load_spectrum_list_with_labels(self):
         labels = ["List test 1", "List test 2", "List test 3"]
         self.spec_app.load_spectrum(self.spec_list, data_label=labels)
-        assert len(self.spec_app.app.data_collection) == 4
+        assert len(self.spec_app.app.data_collection) == 5
 
     def test_mismatched_label_length(self):
         with pytest.raises(ValueError, match='Length'):
@@ -224,3 +227,19 @@ def test_app_links(specviz_helper, spectrum1d):
     sv = specviz_helper.app.get_viewer('spectrum-viewer')
     assert isinstance(sv.jdaviz_app, Application)
     assert isinstance(sv.jdaviz_helper, Specviz)
+
+
+@pytest.mark.remote_data
+def test_load_spectrum_list_directory(tmpdir, specviz_helper):
+    test_data = 'https://stsci.box.com/shared/static/l2azhcqd3tvzhybdlpx2j2qlutkaro3z.zip'
+    fn = download_file(test_data, cache=True, timeout=30)
+    with ZipFile(fn, 'r') as sample_data_zip:
+        sample_data_zip.extractall(tmpdir.strpath)
+    data_path = str(tmpdir.join('NIRISS_for_parser_p0171'))
+
+    with pytest.warns(UserWarning, match='SRCTYPE is missing or UNKNOWN in JWST x1d loader'):
+        specviz_helper.load_spectrum(data_path)
+    assert len(specviz_helper.app.data_collection) == 3
+    for element in specviz_helper.app.data_collection:
+        assert element.data.main_components[0] in ['flux']
+        assert element.data.main_components[1] in ['uncertainty']

notebooks/concepts/Specviz_from_List.ipynb

@@ -0,0 +1,116 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "17d43695",
+   "metadata": {},
+   "source": [
+    "# SpectrumList reading by Specviz"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7d96dae3",
+   "metadata": {},
+   "source": [
+    "This demonstrates the reading of a set of x1d files by specviz, using internally an instance of SpectrumList.\n",
+    "\n",
+    "The desired result is the combination (splicing) of all spectra in the list."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "96971b54",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "import glob\n",
+    "import tempfile\n",
+    "import urllib\n",
+    "import tarfile\n",
+    "from zipfile import ZipFile\n",
+    "\n",
+    "from astropy.utils.data import download_file\n",
+    "\n",
+    "from specutils import SpectrumList\n",
+    "\n",
+    "from jdaviz import Specviz"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5dbebe90",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Unzip files if they haven't already been unzipped\n",
+    "if os.path.exists(\"reduced/\"):\n",
+    "    print(\"Pipeline 3 Data Exists\")\n",
+    "else:\n",
+    "    url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/reduced.tar.gz'\n",
+    "    urllib.request.urlretrieve(url, './reduced.tar.gz')\n",
+    "\n",
+    "    tar = tarfile.open('./reduced.tar.gz', \"r:gz\")\n",
+    "    tar.extractall()\n",
+    "    tar.close()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f5657b71",
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "specviz = Specviz()\n",
+    "specviz.app"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4b3a8944",
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "specviz.load_spectrum('reduced/')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e03a3364",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.9.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}