# Licensed under a 3-clause BSD style license - see LICENSE.rst
# -*- coding: utf-8 -*-
"""
===============
prospect.viewer
===============
Run a spectral viewer (plot spectra and show widgets).
Spectra can be:
- DESI Spectra or Frames,
- `specutils`_-compatible objects (see :mod:`prospect.specutils` for objects and IO routines).
.. _`specutils`: https://specutils.readthedocs.io
"""
import os
import numpy as np
from numpy.ma.core import MaskedConstant
import bokeh.plotting as bk
_specutils_imported = True
try:
from specutils import SpectrumList
except ImportError:
_specutils_imported = False
else:
try:
from specutils import Spectrum
except ImportError:
# support specutils 1.x
from specutils import Spectrum1D as Spectrum
_desispec_imported = True
try:
import desispec.io
import desispec.spectra
import desispec.frame
from desispec.interpolation import resample_flux
except ImportError:
_desispec_imported = False
_redrock_imported = True
try:
import redrock.templates
from redrock.archetypes import All_archetypes
except ImportError:
_redrock_imported = False
from ..utilities import frames2spectra, create_zcat_from_redrock_cat, load_redrock_templates
from .cds import ViewerCDS
from .plots import ViewerPlots
from .widgets import ViewerWidgets
from .vi_widgets import ViewerVIWidgets
from .layouts import ViewerLayout, StandaloneThumbLayout
[docs]
def create_model(spectra, zcat, archetype_fit=False, archetypes_dir=None, template_dir=None):
'''
Returns model_wave[nwave], model_flux[nspec, nwave], row matched to zcat,
which can be in a different order than spectra.
- zcat must be entry-matched to spectra.
'''
assert _redrock_imported
assert _desispec_imported # for resample_flux
if np.any(zcat['TARGETID'] != spectra.fibermap['TARGETID']) :
raise ValueError('zcatalog and spectra do not match (different targetids)')
#- Get template versions from header if possible
if hasattr(zcat, 'meta') and ('TEMNAM00' in zcat.meta):
zcat_header = zcat.meta
else:
zcat_header = None
if archetype_fit:
archetypes = All_archetypes(archetypes_dir=archetypes_dir).archetypes
else:
templates = load_redrock_templates(template_dir=template_dir, zcat_header=zcat_header)
#- Empty model flux arrays per band to fill
model_flux = dict()
for band in spectra.bands:
model_flux[band] = np.zeros(spectra.flux[band].shape)
for i in range(len(zcat)):
zb = zcat[i]
if archetype_fit:
archetype = archetypes[zb['SPECTYPE']]
coeff = zb['COEFF']
for band in spectra.bands:
wave = spectra.wave[band]
wavehash = hash((len(wave), wave[0], wave[1], wave[-2], wave[-1], spectra.R[band].data.shape[0]))
dwave = {wavehash: wave}
mx = archetype.eval(zb['SUBTYPE'], dwave, coeff, wave, zb['Z']) * (1+zb['Z'])
model_flux[band][i] = spectra.R[band][i].dot(mx)
else:
if isinstance(zb['SUBTYPE'], MaskedConstant):
subtype = zcat['SUBTYPE'].data.data[i].decode('utf-8')
else:
subtype = zb['SUBTYPE']
tx = templates[(zb['SPECTYPE'], subtype)]
coeff = zb['COEFF'][0:tx.nbasis]
model = tx.flux.T.dot(coeff).T
for band in spectra.bands:
mx = resample_flux(spectra.wave[band], tx.wave*(1+zb['Z']), model)
model_flux[band][i] = spectra.R[band][i].dot(mx)
#- Combine, if needed, to a single wavelength grid across all cameras
# The following works for any given set of bands
keep = dict()
meanwaves = [ np.mean(spectra.wave[band]) for band in spectra.bands] # sort bands by increasing wave
sorted_bands = [ x for _,x in sorted(zip(meanwaves, spectra.bands)) ]
for i_band, band in enumerate(sorted_bands):
wavecut_low = 0
wavecut_up = 1.e10
if i_band>0:
band_low = sorted_bands[i_band-1]
wavecut_low = 0.5*(spectra.wave[band_low][-1] + spectra.wave[band][0])
if i_band<len(sorted_bands)-1:
band_up = sorted_bands[i_band+1]
wavecut_up = 0.5*(spectra.wave[band][-1] + spectra.wave[band_up][0])
keep[band] = (spectra.wave[band]>wavecut_low) & (spectra.wave[band]<wavecut_up)
model_wave = np.concatenate([spectra.wave[band][keep[band]] for band in sorted_bands])
mflux = np.concatenate([model_flux[band][:, keep[band]] for band in sorted_bands],
axis=1)
return model_wave, mflux
[docs]
def plotspectra(spectra, zcatalog=None, redrock_cat=None, notebook=False,
html_dir=None, outfile=None, title=None,
colors=None, with_imaging=True, with_noise=True, with_thumb_tab=True, with_vi_widgets=True,
top_metadata=None, vi_countdown=-1, with_thumb_only_page=False,
with_coaddcam=True, mask_type='DESI_TARGET',
model_from_zcat=True, model=None, with_other_model=True,
num_approx_fits=None, with_full_2ndfit=True,
template_dir=None, archetype_fit=False, archetypes_dir=None,
std_template_file=None, zmax_slider=5.0):
'''Main prospect routine. From a set of spectra, creates a bokeh document
used for VI, to be displayed as an HTML page or within a Jupyter notebook.
Parameters
----------
spectra : :class:`~desispec.spectra.Spectra` or :class:`~specutils.Spectrum` or :class:`~specutils.SpectrumList` or list of :class:`~desispec.frame.Frame`
Input spectra. :class:`~specutils.Spectrum` are assumed to be SDSS/BOSS/eBOSS.
Otherwise DESI spectra or frames is assumed.
zcatalog : :class:`~astropy.table.Table`, optional
Redshift values, matched one-to-one with the input spectra.
redrock_cat : :class:`~astropy.table.Table`, optional
Redrock output (as defined in :func:`~prospect.utilities.match_rrdetails_to_spectra`).
Entries must be matched one-by-one (in order) to spectra.
notebook : :class:`bool`, optional
If ``True``, bokeh outputs the viewer to a Jupyter notebook.
html_dir : :class:`str`, optional
Directory to store the HTML page if `notebook` is ``False``.
outfile : :class:`str`, optional
Output filename to write, including path
title : :class:`str`, optional
Title used to name the HTML page / the bokeh figure / the VI file.
colors : :class:`list`, optional
Customize the curve's colors: 3 colors should be given, associated respectively to
the coadded data, the model and the noise.
with_imaging : :class:`bool`, optional
If ``False``, don't include thumb image from https://www.legacysurvey.org/viewer.
with_noise : :class:`bool`, optional
If ``False``, don't include uncertainty for each spectrum.
with_thumb_tab : :class:`bool`, optional
If ``False``, don't include a tab with spectra thumbnails.
with_vi_widgets : :class:`bool`, optional
Include widgets used to enter VI information. Set it to ``False`` if
you do not intend to record VI files.
top_metadata : :class:`list`, optional
List of metadata to be highlighted in the top (most visible) table.
Default values ['TARGETID', 'EXPID', 'COADD_NUMEXP', 'COADD_EXPTIME']
vi_countdown : :class:`int`, optional
If ``>0``, add a countdown widget in the VI panel, with a value in minutes given
by `vi_countdown``.
with_thumb_only_page : :class:`bool`, optional
When creating a static HTML (`notebook` is ``False``), a light HTML
page including only the thumb gallery will also be produced.
with_coaddcam : :class:`bool`, optional
Include camera-coaddition, only relevant for DESI.
mask_type : :class:`str`, optional (default: DESI_TARGET)
Bitmask type to identify target categories in the spectra.
Supported types are in `..utilities.supported_desitarget_masks`
model_from_zcat : :class:`bool`, optional
If ``True``, model spectra will be computed from the input `zcatalog`.
model : :func:`tuple`, optional
If set, use this input set of model spectra instead of computing it from `zcatalog`.
model consists of (mwave, mflux); model must be entry-matched to `zcatalog`.
num_approx_fits : :class:`int`, optional
Number of best-fit models to display, if `redrock_cat` is provided.
By default, all best-fit models available in `redrock_cat` are diplayed.
with_full_2ndfit : :class:`bool`, optional
If ``True``, the second best-fit model from `redrock_cat` will be displayed
without approximation (no undersampling, full resolution).
template_dir : :class:`str`, optional
Redrock template directory.
archetype_fit : :class:`bool`, optional
If ``True``, assume `zcatalog` derived from :command:`redrock --archetypes`
and plot model accordingly.
archetypes_dir : :class:`str`, optional
Directory path for archetypes if not :envvar:`RR_ARCHETYPE_DIR`.
std_template_file : :class:`str`, optional
File containing standard templates to display in viewer.
zmax_slider : :class:`float`, optional
Maximum range of the redshift slider widget.
'''
#- Check input spectra.
#- Set masked bins to NaN for compatibility with bokeh.
if _specutils_imported and isinstance(spectra, Spectrum):
# We will assume this is from an SDSS/BOSS/eBOSS spPlate file.
survey = 'SDSS'
nspec = spectra.flux.shape[0]
# Historically, SDSS ignored any mask when marking bad pixels in plots.
bad = (spectra.uncertainty.array == 0.0)
spectra.flux[bad] = np.nan
elif _specutils_imported and isinstance(spectra, SpectrumList):
# We will assume this is from a DESI spectra or coadd file.
survey = 'DESI'
nspec = spectra[0].flux.shape[0]
for s in spectra:
# For DESI, anything that has a non-zero mask should also already
# have ivar == 0, so this may be redundant, but should also be harmless.
bad = (s.uncertainty.array == 0.0) | s.mask
s.flux[bad] = np.nan
else:
# DESI object (Spectra or list of Frame)
survey = 'DESI'
if _desispec_imported and isinstance(spectra, desispec.spectra.Spectra):
nspec = spectra.num_spectra()
elif _desispec_imported and isinstance(spectra, list) and isinstance(spectra[0], desispec.frame.Frame):
# If inputs are frames, convert to a spectra object
spectra = frames2spectra(spectra)
nspec = spectra.num_spectra()
if title is None:
title = 'Night {} ExpID {} Spectrograph {}'.format(
spectra.meta['NIGHT'], spectra.meta['EXPID'], spectra.meta['CAMERA'][1],
)
else:
raise ValueError("Unsupported type for input spectra. \n"+
" _specutils_imported = "+str(_specutils_imported)+"\n"+
" _desispec_imported = "+str(_desispec_imported))
for band in spectra.bands:
# For DESI, anything that has a non-zero mask should also already
# have ivar == 0, so this may be redundant, but should also be harmless.
bad = (spectra.ivar[band] == 0.0) | (spectra.mask[band] != 0)
spectra.flux[band][bad] = np.nan
#- No coaddition if spectra is already single-band
if len(spectra.bands)==1 : with_coaddcam = False
if title is None:
title = "prospect"
#- Input zcatalog / model
if zcatalog is not None:
if survey == 'SDSS':
if len(zcatalog) != spectra.flux.shape[0]:
raise ValueError('zcatalog and spectra do not match (different lengths)')
else:
if np.any(zcatalog['TARGETID'] != spectra.fibermap['TARGETID']) :
raise ValueError('zcatalog and spectra do not match (different targetids)')
if model is not None:
# SDSS spectra will supply the model.
assert not model_from_zcat
mwave, mflux = model
if len(mflux) != nspec:
raise ValueError("model fluxes do not match spectra (different nb of entries)")
if model_from_zcat :
# DESI spectra will obtain the model from templates.
model = create_model(spectra, zcatalog,
archetype_fit=archetype_fit,
archetypes_dir=archetypes_dir,
template_dir=template_dir)
#-----
#- Gather information into ColumnDataSource objects for Bokeh
viewer_cds = ViewerCDS()
viewer_cds.load_spectra(spectra, with_noise)
viewer_cds.compute_median_spectra(spectra)
if with_coaddcam :
viewer_cds.init_coaddcam_spec(spectra, with_noise)
if model is not None:
viewer_cds.init_model(model)
if zcatalog is not None and with_other_model:
viewer_cds.init_othermodel(zcatalog)
if with_other_model:
viewer_cds.load_std_templates(std_template_file=std_template_file)
if redrock_cat is not None :
if np.any(redrock_cat['TARGETID'] != spectra.fibermap['TARGETID']) :
raise RuntimeError('redrock_cat and spectra do not match (different targetids)')
if zcatalog is None :
raise ValueError('Redrock_cat was provided but not zcatalog.')
#- Get template versions from header if possible
if hasattr(redrock_cat, 'meta') and ('TEMNAM00' in redrock_cat.meta):
zcat_header = redrock_cat.meta
else:
zcat_header = None
if num_approx_fits!=0 :
# TODO un-hardcode nbpts_templates ?
viewer_cds.load_fit_templates(template_dir=template_dir, nbpts_templates=4000, zcat_header=zcat_header)
viewer_cds.load_rrdetails(redrock_cat)
#- define num_approx_fits, used in the "model_select" widget:
nfits_redrock_cat = viewer_cds.dict_rrdetails['Nfit']
if num_approx_fits is None : num_approx_fits = nfits_redrock_cat
if (num_approx_fits > nfits_redrock_cat):
print("Warning, num_approx_fits too large wrt redrock_cat. Changing it.")
num_approx_fits = nfits_redrock_cat
if with_full_2ndfit :
zcat_2ndfit = create_zcat_from_redrock_cat(redrock_cat, fit_num=1)
model_2ndfit = create_model(spectra, zcat_2ndfit, archetype_fit=archetype_fit,
archetypes_dir=archetypes_dir, template_dir=template_dir)
viewer_cds.init_model(model_2ndfit, second_fit=True)
viewer_cds.load_metadata(spectra, mask_type=mask_type, zcatalog=zcatalog, survey=survey)
#-------------------------
#-- Graphical objects --
#-------------------------
viewer_plots = ViewerPlots(colors=colors)
viewer_plots.create_mainfig(spectra, title, viewer_cds, survey,
with_noise=with_noise, with_coaddcam=with_coaddcam)
viewer_plots.create_zoomfig(viewer_cds,
with_noise=with_noise, with_coaddcam=with_coaddcam)
if with_imaging :
viewer_plots.create_imfig(spectra)
#-----
#- Emission and absorption lines
z = zcatalog['Z'][0] if (zcatalog is not None) else 0.0
viewer_cds.load_spectral_lines(z)
viewer_plots.add_spectral_lines(viewer_cds, figure='main')
viewer_plots.add_spectral_lines(viewer_cds, figure='zoom', label_offset_top=50)
#-------------------------
#-- Widgets and callbacks --
#-------------------------
viewer_widgets = ViewerWidgets(viewer_plots, nspec)
viewer_widgets.add_navigation(nspec)
viewer_widgets.add_resetrange(viewer_cds, viewer_plots)
viewer_widgets.add_redshift_widgets(z, viewer_cds, viewer_plots, zmax_slider)
viewer_widgets.add_oii_widgets(viewer_plots)
viewer_plots.add_imfig_callback(viewer_widgets)
if viewer_cds.cds_coaddcam_spec is not None :
viewer_widgets.add_coaddcam(viewer_plots)
if zcatalog is not None :
show_zcat = True
else : show_zcat = False
if top_metadata is None: top_metadata = ['TARGETID', 'EXPID', 'COADD_NUMEXP', 'COADD_EXPTIME']
viewer_widgets.add_metadata_tables(viewer_cds, top_metadata=top_metadata,
show_zcat=show_zcat)
viewer_widgets.add_specline_toggles(viewer_cds, viewer_plots)
if with_other_model:
viewer_widgets.add_model_select(viewer_cds, num_approx_fits)
#-----
#- VI-related widgets
## TODO if with_vi_widgets (need to adapt update_plot.js..)
viewer_vi_widgets = ViewerVIWidgets(title, viewer_cds)
viewer_vi_widgets.add_filename()
viewer_vi_widgets.add_vi_issues(viewer_cds, viewer_widgets)
viewer_vi_widgets.add_vi_z(viewer_cds, viewer_widgets)
viewer_vi_widgets.add_vi_spectype(viewer_cds, viewer_widgets)
viewer_vi_widgets.add_vi_comment(viewer_cds, viewer_widgets)
viewer_vi_widgets.add_vi_quality(viewer_cds, viewer_widgets)
viewer_vi_widgets.add_vi_scanner(viewer_cds)
viewer_vi_widgets.add_guidelines()
viewer_vi_widgets.add_vi_storage(viewer_cds, viewer_widgets)
viewer_vi_widgets.add_vi_table(viewer_cds)
if (vi_countdown > 0) :
viewer_vi_widgets.add_countdown(vi_countdown)
viewer_widgets.add_update_plot_callback(viewer_cds, viewer_plots,
viewer_vi_widgets)
#-----
#- Bokeh layout and output
bokeh_layout = ViewerLayout(viewer_plots, viewer_widgets, viewer_vi_widgets,
with_vi_widgets=with_vi_widgets)
if with_thumb_tab:
bokeh_layout.add_thumb_tab(spectra, viewer_plots, viewer_widgets, nspec)
if notebook:
bk.output_notebook()
bk.show(bokeh_layout.full_viewer)
else:
if html_dir is None and outfile is None:
raise RuntimeError("Need html_dir or outfile")
if outfile is None:
outfile = os.path.join(html_dir, title+".html")
bk.output_file(outfile, title='DESI spectral viewer')
bk.save(bokeh_layout.full_viewer)
#-----
#- "Light" Bokeh layout including only the thumbnail gallery
if with_thumb_only_page :
assert not notebook
thumb_page = os.path.join(html_dir, "thumbs_"+title+".html")
bk.output_file(thumb_page, title='DESI spectral viewer - thumbnail gallery')
thumb_grid = StandaloneThumbLayout(spectra, viewer_plots, title)
bk.save(thumb_grid.thumb_viewer)
return