Introducing spectral extractions for non-flat traces

[bmorris3]

Description

This PR implements Horne spectral extraction with non-flat traces. For a demo, see this example notebook.

specreduce Horne extractions currently support flat traces. In the flat trace case, the sum of the spectrum in the dispersion dimension produces a flux profile in the spatial dimension, which we fit with a (1D) Gaussian. The model profile is used to construct a kernel, which applies weights to each pixel in the extracted 1D spectrum.

For a non-flat trace, the profile in the will be "smeared" when summed in the dispersion dimension, which is not accurately approximated by a Gaussian. Less flat traces will produce poorer kernels, and less accurate spectral extractions.

The trick used here to support non-flat traces is to roll the pixels along the spatial dimension until the trace is flat, and then perform the extraction on the shifted image as usual. There are no actual changes to the Horne implementation, just a manipulation of the image in the case of a non-flat trace to "make it flat".

The trick within the trick is that we can use fancy indexing to flatten the trace without interpolation or resampling, and do it efficiently. In this PR, we only re-index the underlying image.

My strategy here is to (1) never make up data (no interp/sampling business), and (2) keep the diff small.

Possible Improvements

• We'll need some tests and documentation, of course. (I'm posting this draft PR to get initial feedback.)
• @kecnry pointed out a possible (small?) improvement. This PR casts floats in the trace to integers for rolling the pixels along the spatial dimension (without interpolation). That means we can expect a small spatial difference between the center of the spatial profile and the intended center of the new flat trace. We could pass along that small difference to the initial guess for the mean in the Gaussian 1D spatial profile. This may not affect the results much, because the rolls are accurate to better than a pixel, and the spatial width of the observed trace is usually a few pixels -- so the initial guess is usually good enough.

Closes #102

🐱

[codecov]

Codecov Report

Merging #148 (fc525be) into main (05a174d) will decrease coverage by 2.47%.
The diff coverage is 92.85%.

❗ Current head fc525be differs from pull request most recent head b4c5db9. Consider uploading reports for the commit b4c5db9 to get more accurate results

@@            Coverage Diff             @@
##             main     #148      +/-   ##
==========================================
- Coverage   75.46%   72.98%   -2.48%     
==========================================
  Files           9        9              
  Lines         693      659      -34     
==========================================
- Hits          523      481      -42     
- Misses        170      178       +8     

Impacted Files	Coverage Δ
specreduce/extract.py	88.66% <92.85%> (-6.64%)	⬇️
specreduce/core.py	47.05% <0.00%> (-27.31%)	⬇️
specreduce/tracing.py	92.76% <0.00%> (-0.79%)	⬇️
specreduce/background.py	89.28% <0.00%> (+0.98%)	⬆️

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[kecnry]

would just swapping this out with the "remainder" from casting trace_object.trace to their integer rolls be enough to correct the means to the actual trace position?

[bmorris3]

I think it would improve the initial guess for the mean parameter given to the Gaussian spatial profile fit by <1 pixel.

[ojustino]

This is a creative approach that seems to work upon initial investigation. I would like to test the extractions more thoroughly, but I have some comments on the code in the meantime.

[kecnry]

I think this is ready for a final code-review once test coverage is added.

[bmorris3]

The remaining failure comes from astropy convolution, see astropy/astropy#14025

[ojustino]

It looks like the big difference since my last check is the test, which seems reasonable.

I put more thought into the resampling procedure and am wondering if it's problematic in cases where it's not actually needed because the trace is truly flat. The following is the same as cell 4 of your gist's notebook except it skips the transformation step:

image_nddata = NDData(data, uncertainty=error)
flat_trace = FlatTrace(data, 30)
spec1d_init = HorneExtract(image_nddata, trace)()

fig, ax = plot_result(image_nddata, flat_trace, spec1d_init)

Is the difference in the results a cause for concern?

[bmorris3]

Hi @ojustino, thanks for your investigation and this comment. The "untransformed" image should indeed produce a distinct spectral extraction from the one I made in the example notebook which does transform image_nddata via transform_flat_to_custom_trace. That function does some interpolation and resampling of the "original" image which takes advantage of interpolation and resampling via scipy.ndimage.zoom. This method is not flux conserving and not written in a way so that the "untransformed" image matches the "transformed" one when the trace is flat. I only wrote that method to fudge non-flat traces from a real and flat example.

[bmorris3]

I've reverted the very last commit. The reverted commit pinned the setuptools version, due to an upstream astropy failure related to the convolution module. The astropy fix was merged in astropy/astropy#14035, so we no longer need the pin here.

[ojustino]

@bmorris3, that makes sense, thanks. Looks like I was fixating on the example instead of the actual code. The code looks good to me and the test makes sense, so I approve from that perspective. I would be interested in seeing feedback from others with more experience in extraction using non-flat traces, but I have no objections based on the ones I've tried.

[kecnry]

I agree - let's get this in and make any small tweaks once there is some user feedback!

[bmorris3]

Thanks both for the reviews. 🎉