Skip to content

Added splice module#152

Open
ladsantos wants to merge 20 commits into
spacetelescope:masterfrom
ladsantos:dev_ldsantos
Open

Added splice module#152
ladsantos wants to merge 20 commits into
spacetelescope:masterfrom
ladsantos:dev_ldsantos

Conversation

@ladsantos

Copy link
Copy Markdown

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

@sean-lockwood

Copy link
Copy Markdown
Member

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

ladsantos commented Mar 2, 2023

Copy link
Copy Markdown
Author

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

@sean-lockwood

sean-lockwood commented Mar 2, 2023

Copy link
Copy Markdown
Member

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

sean-lockwood commented Mar 14, 2023

Copy link
Copy Markdown
Member

Docstrings should follow numpy/scipy docstring format:

image

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.

Comment thread stistools/splice.py Outdated
dq_weights_interp = np.zeros_like(dq_interp)
# And then for each acceptable dq, if the element of the dq array is one
# of the acceptable flags, we set its dq weight to one
for adq in acceptable_dq_flags:

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 sean-lockwood Mar 31, 2023

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 ladsantos Apr 5, 2023

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

@ladsantos ladsantos requested a review from a team as a code owner April 25, 2023 15:08
@sean-lockwood

Copy link
Copy Markdown
Member

@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

Copy link
Copy Markdown
Author

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

@sean-lockwood

Copy link
Copy Markdown
Member

@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.

image

Full range:
splice_bug_full

cc: @Jackie-Brown

@sean-lockwood

Copy link
Copy Markdown
Member

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

image

Might be indicative of a problem.

@ladsantos

Copy link
Copy Markdown
Author

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.

Screenshot 2024-05-23 at 3 27 48 PM

@ladsantos

Copy link
Copy Markdown
Author

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

@ladsantos

ladsantos commented May 23, 2024

Copy link
Copy Markdown
Author

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.

Screenshot 2024-05-23 at 3 41 24 PM

@ladsantos

Copy link
Copy Markdown
Author

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.

Screenshot 2024-05-24 at 11 15 14 AM

@sean-lockwood

sean-lockwood commented Aug 14, 2024

Copy link
Copy Markdown
Member

@ladsantos -

The spikes on the right have DQ = (512 | 2048), whereas the ones on the left are 2^15.

I think making the current DQ 2^15 regions be (2^15 | 2048) makes the most sense, as we'd be able to filter all the spikes using the prior definition of SDQFLAGS (preserving behavior in legacy code) and your distinction using 2^15 would still be present.

(Pending discussion with the team.)

splice_dq_states_zoomed

@sean-lockwood sean-lockwood left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm still looking at the DQ behavior in more detail, but here are a few issues I noticed today.

Comment thread stistools/splice.py Outdated
# We interpolate the lower-SNR spectra to the wavelength bins of the higher
# SNR spectrum.
max_sens_idx = np.where(avg_sensitivity == np.nanmax(avg_sensitivity))[0][0]
overlap_ref = overlap_sections.pop(max_sens_idx)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This has the side effect of modifying the input variable. You should probably make a copy inside the function instead.

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Corrected in commit a2f5407.

Comment thread stistools/splice.py Outdated
if truncate_edge_right is not None:
for sk in spectrum:
sk['data_quality'][-truncate_edge_right:] = 4096
elif truncate_edge_left is not None:

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

elif -> if allows both edges to be truncated with default values. Is that the desired behavior?

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Corrected in commit a2f5407.

Comment thread stistools/splice.py Outdated

# Merge the overlapping spectral sections
merged_pairs = [
merge_overlap(overlap_pair_sections[k], acceptable_dq_flags)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should the weight parameter be passed here too?

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Corrected in commit a2f5407.

Comment thread stistools/splice.py
overlap_s_u = np.copy(overlap_sections[i]['uncertainty'])
overlap_s_n = np.copy(overlap_sections[i]['net'])
overlap_s_dq = overlap_sections[i]['data_quality']
where_dq_bad = np.where(overlap_s_dq & bitwise_or_acceptable_dq_flags)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If bitwise_or_acceptable_dq_flags are the flags that we're ok ignoring, shouldn't we use the inverse here to find bad locations?

unacceptable_dq_flags = int(~np.uint16(bitwise_or_acceptable_dq_flags))
...

where_dq_bad = (overlap_s_dq & unacceptable_dq_flags) != 0

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried this, and it does not remove the spikes on the right. I think there's something we are missing with the logic here (see my comment below).

@ladsantos

Copy link
Copy Markdown
Author

@ladsantos -

The spikes on the right have DQ = (512 | 2048), whereas the ones on the left are 2^15.

I think making the current DQ 2^15 regions be (2^15 | 2048) makes the most sense, as we'd be able to filter all the spikes using the prior definition of SDQFLAGS (preserving behavior in legacy code) and your distinction using 2^15 would still be present.

(Pending discussion with the team.)

The 2^15 DQ flag represents pixels that were co-added from good pixels (those with acceptable flags). This is odd, because the 2048 and 512 are not an acceptable DQ flag, so it shouldn't be co-adding the DQ = (512 | 2048) pixels. Or am I getting this logic wrong?

@sean-lockwood

Copy link
Copy Markdown
Member

I've been working on splice and uploaded a draft (still in progress) to this branch:
https://github.com/spacetelescope/stistools/tree/sl_splice

Among other changes, this fixes DQ flags and implements sensitivity weighting derived from reference files. Algorithm improvements and tests are still in progress.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants