Added splice module

[ladsantos]

Added the splice module, which is based on the standalone code stissplice, including a simple test suite and documentation.

[sean-lockwood]

Testing would be more robust if you upload an input file and a "truth" output file. Then you could compare all data values.

I'm concerned that testing one flux value might not be sensitive to changes in all of your code's logic.

[ladsantos]

Improved the test suite of splice in commit bece227. Now the whole spectrum is tested against a truth file.

[sean-lockwood]

Leo,

Will most users want to use stistools.splice.splice() or stistools.splice.splice_pipeline()?

I'm wondering if splice_pipeline should be renamed to splice and the original splice to something like splice_overlap (or something else descriptive). I know that if I were trying it for the first time, I would either check stistools.splice() (which is actually module-level and not how routines are stored in stistools) or stistools.splice.splice, following the pattern set by most of the other routines.

[sean-lockwood]

Docstrings should follow numpy/scipy docstring format:

Current:

    Parameters
    ----------
    unique_spectra_list (``list``):
        List of unique spectra.

    merged_pair_list (``list``):
        List of merged overlapping pair spectra.
    ...

Suggested:

    Parameters
    ----------
    unique_spectra_list : list
        List of unique spectra.

    merged_pair_list : list
        List of merged overlapping pair spectra.
    ...

Note that the colons won't currently render properly on RTD, but we're working on that in a separate PR.

[sean-lockwood]

Should probably check bitwise, i.e.:

np.where(dq_ref & adq)

and

np.where(dq_interp & adq)

In case more than one flag is present at a particular location. Perhaps this doesn't occur with the default flags, but could with other options provided by the user.

[sean-lockwood]

In fact, you can probably skip the loop over the various acceptable_dq_flags values and compare to the scalar value that is the bitwise-or or the adq flags (sum in this case).

[sean-lockwood]

Ok, so wrapping it up together:

from functools import reduce

acceptable_dq_flags = reduce(np.bitwise_or, acceptable_dq_flags)  # scalar value

dq_weights_ref[np.where(dq_ref & acceptable_dq_flags)] = 1
dq_weights_interp[np.where(dq_interp & acceptable_dq_flags)] = 1

(Should be tested.)

[ladsantos]

It seems that this code snippet catches all the good flags, except for 0. But I should be able to add a line or two that catches it.

[sean-lockwood]

Good catch! We don't typically flag acceptable pixels, so that slipped by.

Perhaps:

bad_dq_flags = ~reduce(np.bitwise_or, acceptable_dq_flags)  # scalar value

dq_weights_ref = np.ones(..., dtype=...)
dq_weights_interp = np.ones(..., dtype=...)

dq_weights_ref[dq_ref & bad_dq_flags] = 0
dq_weights_interp[dq_interp & bad_dq_flags] = 0

Looking at it now, we probably don't need the np.where either.

[sean-lockwood]

@ladsantos -
Sorry, I lost track of where we were on this. Do you think this is ready for some sanity checks and almost ready to merge?

[ladsantos]

Hi @sean-lockwood, no worries! It is good to go for sanity checks and merge.

[sean-lockwood]

@ladsantos -
We just did a quick test on oewia2010_x1d.fits and found the splice code is introducing spikes at the end of some orders. Applying DQ filtering doesn't seem to help.

Full range:

cc: @Jackie-Brown

[sean-lockwood]

I'm noticing near the spikes there are a lot of DQ = 32768 = 2**15, which isn't defined in the STIS DQ values.

Might be indicative of a problem.

[ladsantos]

It seems to me that splice is not really "introducing" these spikes, they are already there in the x1d files in the first place. It is difficult to correct them without manual inspection/manipulation. If I change the weight parameter to 'snr' instead of 'sensitivity' (which is the default), and I remove the DQ flag 2048 as an acceptable flag, the results are much better -- although not completely alleviated.

[ladsantos]

Also, the DQ flag 32768 simply means a merged pixel where there was an overlap.

[ladsantos]

Here's a closer look at the spikes at longer wavelengths (spliced spectrum is offset for visualization purposes). It seems to me that the issue lies in the x1d orders not overlapping, so there is no merging to be done to correct for the edge effects. In such cases, there is not much that the splice code can do besides just concatenating the orders.

[ladsantos]

Found a solution to the spikes near the order edges. The new version of the code simply adds a DQ flag of 4096 at a user-specified edge truncation number. Note that it does not fix spikes that are present in the x1d data when there is no overlap.