stissplice documentation

stissplice is the splicer of Echelle Spectra from Hubble Space Telescope. This code splices Echelle spectra obtained with the Space Telescope Imaging Spectrograph (STIS) instrument. It can be adapted to work with spectra obtained with other instruments as well.

This code is currently in development and accepting suggestions of features to implement, as well as contributions.

Installation

Option 1: Using pip (stable version)

Simply run the following command:

pip install stissplice

Option 2: Compile from source (development version)

First, clone the repository and then navigate to it:

git clone https://github.com/ladsantos/stissplice
cd stissplice

And then compile it from source:

python setup.py install

You can test the installation from source with pytest (you may need to install pytest first):

pytest tests

stissplice API

splicer

splicer.read_spectrum(filename, prefix, truncate_edge_left=None, truncate_edge_right=None)

This is a fairly straightforward function to read the spectrum from a x1d FITS file.

Parameters:

filename (``str``):

Name of the *_x1d.fits file containing the spectrum.

prefix (``str``):

Path to the *_x1d.fits file containing the spectrum.

truncate_edge_left (``int``, optional):

Set the number of low-resolution pixels at the left edge of the detector where the spectra should be truncated. If None, then no truncation is applied. Default is None.

truncate_edge_right (``int``, optional):

Set the number of low-resolution pixels at the right edge of the detector where the spectra should be truncated. If None, then no truncation is applied. Default is None.

Returns:

spectrum (list):

List of all the orders contained in the Echelle spectrum and their respective fluxes.

splicer.find_overlap(spectrum)

Find and return the overlapping sections of the Echelle spectrum.

Parameters:

spectrum (``list``):

List of dictionaries containing the orders of the Echelle spectrum. It should resemble the output of read_spectrum().

Returns:

unique_sections (list):

List containing the unique sections of the spectrum.

overlap_pair_sections (list):

List containing the overlapping pairs of the spectrum.

overlap_trio_sections (list):

List containing the overlapping trios of the spectrum.

splicer.merge_overlap(overlap_sections, acceptable_dq_flags=(0, 64, 128, 1024, 2048), weight='sensitivity')

Merges overlapping spectral regions. The basic workflow of this function is to interpolate the sections into a common wavelength table and calculate the weighted mean flux for each wavelength bin. If the fluxes are inconsistent between each other, the code can use the flux with higher SNR instead of the mean. If there are still outlier fluxes (compared to neighboring pixels), the code uses the flux from the lower SNR section instead. Co-added (merged) pixels will have their DQ flag set to 32768 if they are the result of combining good pixels (according to the list of acceptable flags). Their DQ flag will be set to 65536 if the combined pixels do not have an acceptable DQ flag.

Parameters:

overlap_sections (``list``):

List of dictionaries containing the overlapping spectra of neighboring orders.

acceptable_dq_flags (array-like, optional):

Data-quality flags that are acceptable when co-adding overlapping spectra. The default values are (0, 64, 128, 1024, 2048), which correspond to: 0 = regular pixel, 64 = vignetted pixel, 128 = pixel in overscan region, 1024 = small blemish, 2048 = more than 30% of background pixels rejected by sigma-clipping in the data reduction.

weightstr, optional

Defines how to merge the overlapping sections. The options currently implemented are 'sensitivity' and 'snr' (inverse square of the uncertainties). Default is 'sensitivity'.

Returns:

overlap_merged (dict):

Dictionary containing the merged overlapping spectrum.

splicer.splice(unique_spectra_list, merged_pair_list, merged_trio_list)

Concatenate the unique and the (merged) overlapping spectra.

Parameters:

unique_spectra_list (``list``):

List of unique spectra.

merged_pair_list (``list``):

List of merged overlapping pair spectra.

merged_trio_list (``list``):

List of merged overlapping trio spectra.

Returns:

spliced_wavelength (numpy.ndarray):

Array containing the wavelengths in the entire spectrum.

spliced_flux (numpy.ndarray):

Array containing the fluxes in the entire spectrum.

spliced_uncertainty (numpy.ndarray):

Array containing the flux uncertainties in the entire spectrum.

splicer.splice_pipeline(dataset, prefix='./', update_fits=False, output_file=None, weight='sensitivity', acceptable_dq_flags=(0, 64, 128, 1024, 2048))

The main workhorse of the package. This pipeline performs all the steps necessary to merge overlapping spectral sections and splice them with the unique sections.

Parameters:

dataset (``str``):

String that identifies the dataset (example: 'oblh01040').

prefix (``str``):

Path to the *_x1d.fits file containing the spectrum.

update_fits (``bool``, optional):

Use carefully, since it can modify fits files permanently. Parameter that decides whether to update the *_x1d.fits file with a new extension containing the spliced spectrum.

output_file (``str`` or ``None``, optional):

String containing the location to save the output spectrum as an ascii file. If None, no output file is saved and the code returns an Astropy Table instead. Default is None.

weightstr, optional

Defines how to merge the overlapping sections. The options currently implemented are 'sensitivity' and 'snr' (inverse square of the uncertainties). Default is 'sensitivity'.

acceptable_dq_flags (array-like, optional):

Data-quality flags that are acceptable when co-adding overlapping spectra. The default values are (0, 64, 128, 1024, 2048), which correspond to: 0 = regular pixel, 64 = vignetted pixel, 128 = pixel in overscan region, 1024 = small blemish, 2048 = more than 30% of background pixels rejected by sigma-clipping in the data reduction.

Returns:

spliced_spectrum_table (astropy.Table object):

Astropy Table containing the spliced spectrum. Only returned if output_file is None.

tools

tools.nearest_index(array, target_value)

Finds the index of a value in array that is closest to target_value.

Parameters:

array (``numpy.array``):

Target array.

target_value (``float``):

Target value.

Returns:

index (int):

Index of the value in array that is closest to target_value.

Usage example

In this example notebook we learn the basics of how to use stissplice to splice STIS echèlle spectra measured with HST. We start by importing the necessary packages.

import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
from stissplice import splicer, tools
import numpy as np

# Uncomment the next line if you have a MacBook with retina screen
# %config InlineBackend.figure_format = 'retina'
pylab.rcParams['figure.figsize'] = 9.0,6.5
pylab.rcParams['font.size'] = 18

STIS echèlle spectra are composed of many orders that overlap in their edges. Here’s how it looks when we read an x1d file:

dataset = 'oblh01040'
prefix = '../../data/'

spectrum = splicer.read_spectrum(dataset, prefix)

for s in spectrum:
    plt.plot(s['wavelength'], s['flux'], alpha=0.3)
_ = plt.xlabel(r'Wavelength (${\rm \AA}$)')
_ = plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$)')

If we zoom in, we can see the order overlap in more detail.

for s in spectrum:
    plt.plot(s['wavelength'], s['flux'], alpha=0.5)
_ = plt.xlabel(r'Wavelength (${\rm \AA}$)')
_ = plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$)')
_ = plt.xlim(2200, 2300)
_ = plt.ylim(-.5E-11, 1.5E-11)

The objective of the stissplice package is precisely to deal with the order overlap and co-add them to produce a one-dimensional spectrum. We do that using the splice_pipeline() function.

spliced_spectrum = splicer.splice_pipeline(dataset, prefix, weight='sensitivity')

And now we plot it in black and compare it to the original spectra in colors.

plt.plot(spliced_spectrum['WAVELENGTH'], spliced_spectrum['FLUX'], color='k')

for s in spectrum:
    plt.plot(s['wavelength'], s['flux'], alpha=0.3)
_ = plt.xlabel(r'Wavelength (${\rm \AA}$)')
_ = plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$)')
_ = plt.ylim(-.5E-11, 3E-11)

And here’s how the zoomed-in part looks like:

plt.plot(spliced_spectrum['WAVELENGTH'], spliced_spectrum['FLUX'], color='k')

for s in spectrum:
    plt.plot(s['wavelength'], s['flux'], alpha=0.5)
_ = plt.xlabel(r'Wavelength (${\rm \AA}$)')
_ = plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$)')
_ = plt.xlim(2200, 2300)
_ = plt.ylim(-.5E-11, 1.5E-11)

How it works

In this notebook we learn how the stissplice code works in more detail. We start by importing the necessary packages.

import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
from stissplice import splicer, tools
import numpy as np

# Uncomment the next line if you have a MacBook with retina screen
# %config InlineBackend.figure_format = 'retina'
pylab.rcParams['figure.figsize'] = 13.0,6.5
pylab.rcParams['font.size'] = 18

The workhorse of stissplice is the splice_pipeline() function. Here’s how it works, step by step.

1. Read an _x1d spectrum using splicer.read_spectrum(): This code uses astropy.io.fits.open() to read a .fits file and extract the wavelength, flux and error arrays.

dataset = 'oblh01040'
prefix = '../../data/'

spectrum = splicer.read_spectrum(dataset, prefix)

2. Find the overlap between the different spectral orders using splicer.find_overlap(): This code identifies the sections of the orders that are either unique (no overlap with any order) or that overlap once or twice with other orders. This function returns the variables unique_sections, overlap_pair_sections, overlap_trio_sections.

In the plot below, black represents the unique sections. Red represents sections that are double overlaps (one of them is shifted in the y-axis for clarity purposes). And green represents the sections in triple overlaps (two of them shifted for clarity).

unique_sections, overlap_pair_sections, overlap_trio_sections = splicer.find_overlap(spectrum)

for us in unique_sections:
    plt.plot(us['wavelength'], us['flux'], color='k')

for ops in overlap_pair_sections:
    plt.plot(ops[0]['wavelength'], ops[0]['flux'], color='C3', alpha=0.5)
    plt.plot(ops[1]['wavelength'], ops[1]['flux'] + 1.0E-11, color='C3', alpha=0.5)

for ots in overlap_trio_sections:
    plt.plot(ots[0]['wavelength'], ots[0]['flux'], color='C2', alpha=0.5)
    plt.plot(ots[1]['wavelength'], ots[1]['flux'] - 0.7E-11, color='C2', alpha=0.5)
    plt.plot(ots[2]['wavelength'], ots[2]['flux'] - 1.6E-11, color='C2', alpha=0.5)

_ = plt.xlabel(r'Wavelength (${\rm \AA}$)')
_ = plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$)')

3. Merge the overlapping sections using splicer.merge_overlaps(): This function takes either the double or triple overlap sections calculated above, and co-add them. The co-adding process begins by finding which section has the higher SNR, and then interpolates the other sibling section(s) into the same wavelength grid as the section with the higher SNR. The sibling sections are then merged by taking the weighted-average flux and error for each pixel, where the weights are the inverse-square of the errors. The merged pixels will then be assigned a data quality flag with the value 32768, which correspond to a co-added pixel. This function returns the co-added section

When co-adding pixels, the code identifies the data quality flag (DQ) of each one of them, and only co-adds pixels that have acceptable DQ flags. The user can specify which flags are acceptable. The default acceptable flags are: 0 = regular pixel, 64 = vignetted pixel, 128 = pixel in overscan region, 1024 = small blemish, 2048 = more than 30% of background pixels rejected by sigma-clipping in the data reduction.

merged_pairs = [
    splicer.merge_overlap(overlap_pair_sections[k])
    for k in range(len(overlap_pair_sections))
]

merged_trios = [
    splicer.merge_overlap(overlap_trio_sections[k])
    for k in range(len(overlap_trio_sections))
]

By now we have three lists: unique_sections, merged_pairs and merged_trios. The next step is to concatenate everything in the correct order using the splicer.splice() function.

wavelength, flux, uncertainty, dq = splicer.splice(
    unique_sections, merged_pairs, merged_trios)

# Plot the spliced spectrum
plt.plot(wavelength, flux)
_ = plt.xlabel(r'Wavelength (${\rm \AA}$)')
_ = plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$)')

As shown in the usage example, this whole process can be executed in one line using the splicer.splice_pipeline() function.

spliced_spectrum = splicer.splice_pipeline(dataset, prefix, )
plt.plot(spliced_spectrum['WAVELENGTH'], spliced_spectrum['FLUX'])
_ = plt.xlabel(r'Wavelength (${\rm \AA}$)')
_ = plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$)')