Advanced Spectrum List Patterns and Workarounds

muler is built on specutils, which has many spectrum list operations, but not everything we’d ideally want. In this notebook, we introduce a few workarounds to deal with common patterns.

%config Completer.use_jedi = False

from muler.igrins import IGRINSSpectrum, IGRINSSpectrumList
from specutils import Spectrum1D, SpectrumList
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
%matplotlib inline
%config InlineBackend.figure_format='retina'

Simply tell the IGRINS Spectrum where your file is located.

path = 'https://github.com/OttoStruve/muler_example_data/raw/main/IGRINS/01_IGRINS_test_data/'
filename1='SDCH_20201202_0059.spec_a0v.fits'
full_path = path + filename1

spec_list = IGRINSSpectrumList.read(full_path).normalize()

spec_list.remove_nans().trim_edges().normalize(order_index=15).plot(color=None, ylo=0, yhi=1.5);

Unifying H and K bands into a single IGRINSSpectrumList

filename2 = 'SDCK_20201202_0059.spec_a0v.fits'
spec_list2 = IGRINSSpectrumList.read(path+filename2).normalize()

The extend() method works in-place, overriding the first argument! We have to copy the spec_list if we want to keep it pristine:

import copy

full_H_and_K_spectrum = copy.deepcopy(spec_list)

full_H_and_K_spectrum.extend(spec_list2)

len(spec_list), len(spec_list2), len(full_H_and_K_spectrum)

(28, 26, 54)

full_H_and_K_spectrum.remove_nans().trim_edges().flatten().plot(color=None, ylo=0.8, yhi=1.2);

Workarounds

Here we introduce workarounds to three friction points.

Issue 1: Slicing a IGRINSSpectrumList

First, we want to be able to index into the IGRINSSpectrumList in the same way you would with a regular old Python List, and still get back an IGRINSSpectrumList.

type(spec_list)

muler.igrins.IGRINSSpectrumList

sublist = spec_list[4:7]

type(sublist)

list

No! We want this sublist to be an IGRINSSpectrumList

sublist_IGRINS = IGRINSSpectrumList(sublist)

type(sublist_IGRINS)

muler.igrins.IGRINSSpectrumList

Yay! This workaround achieved our goal.

Issue 2: Create an IGRINSSpectrum from a bare (flux, wavelength) array

Often we have already preprocessed a spectrum and have it in-hand as a pair of wavelength coordinates and flux values. We wish to turn this most elemental format of a spectrum into a IGRINSSpectrum.

import astropy.units as u

whole_spectrum = spec_list.remove_nans().trim_edges().normalize().stitch()

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[18], line 1
----> 1 whole_spectrum = spec_list.remove_nans().trim_edges().normalize().stitch()

File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/muler/echelle.py:881, in EchelleSpectrumList.stitch(self)
    872         unc_anc = None
    874     meta_out[ancillary_spectrum] = spec[0].__class__(
    875         spectral_axis=wls_anc,
    876         flux=fluxes_anc,
    877         uncertainty=unc_anc,
    878         meta=meta_of_meta,
    879     )
--> 881 return spec[0].__class__(
    882     spectral_axis=wls, flux=fluxes, uncertainty=unc_out, meta=meta_out, wcs=None
    883 )

File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/muler/igrins.py:203, in IGRINSSpectrum.__init__(self, file, order, sn_used, cached_hdus, wavefile, *args, **kwargs)
    193     super().__init__(
    194         spectral_axis=lamb.to(u.Angstrom),
    195         flux=flux,
   (...)
    200         **kwargs,
    201     )
    202 else:
--> 203     super().__init__(*args, **kwargs)

File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/muler/echelle.py:74, in EchelleSpectrum.__init__(self, *args, **kwargs)
     71 def __init__(self, *args, **kwargs):
     72
     73     # self.ancillary_spectra = None
---> 74     super().__init__(*args, **kwargs)

File ~/checkouts/readthedocs.org/user_builds/muler/envs/latest/lib/python3.8/site-packages/specutils/spectra/spectrum1d.py:302, in Spectrum1D.__init__(self, flux, spectral_axis, wcs, velocity_convention, rest_value, redshift, radial_velocity, bin_specification, **kwargs)
    299 # make sure that spectral axis is strictly increasing or decreasing,
    300 # raise an error if not.
    301 if not self._check_strictly_increasing_decreasing():
--> 302     raise ValueError('Spectral axis must be strictly increasing or decreasing.')
    304 if hasattr(self, 'uncertainty') and self.uncertainty is not None:
    305     if not flux.shape == self.uncertainty.array.shape:

ValueError: Spectral axis must be strictly increasing or decreasing.

len(whole_spectrum.flux.value)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[19], line 1
----> 1 len(whole_spectrum.flux.value)

NameError: name 'whole_spectrum' is not defined

wavelength = whole_spectrum.wavelength.value
flux = whole_spectrum.flux.value

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[20], line 1
----> 1 wavelength = whole_spectrum.wavelength.value
      2 flux = whole_spectrum.flux.value

NameError: name 'whole_spectrum' is not defined

Can we “round-trip” this pair of coordinates and flux values back into an IGRINSSpectrum?

roundtrip_spec = IGRINSSpectrum(spectral_axis=wavelength*u.Angstrom,
                                flux=flux*u.dimensionless_unscaled, )

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[21], line 1
----> 1 roundtrip_spec = IGRINSSpectrum(spectral_axis=wavelength*u.Angstrom,
      2                                 flux=flux*u.dimensionless_unscaled, )

NameError: name 'wavelength' is not defined

type(roundtrip_spec)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[22], line 1
----> 1 type(roundtrip_spec)

NameError: name 'roundtrip_spec' is not defined

Some methods will work out-of-the-box:

roundtrip_spec.normalize() # works fine!

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[23], line 1
----> 1 roundtrip_spec.normalize() # works fine!

NameError: name 'roundtrip_spec' is not defined

But others do not!

output = whole_spectrum.barycentric_correct() # works fine!

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[24], line 1
----> 1 output = whole_spectrum.barycentric_correct() # works fine!

NameError: name 'whole_spectrum' is not defined

Watch out! This next line will not work!

force = False

if force:
    roundtrip_spec.barycentric_correct() # Oops! Will not work!